check_estimation_result

% check estimation result, PCV Sect. 4.6.8

 check_estimation_result(R,x_true,S0,s0q_sample,x_sample,S,text)

 R          = redundancy (d.o.f. of variance factor)
 x_true     = Ux1 vector of parameters
 C0         = UxU reference covariance matrix of estimated parameters
 s0q_sample = Kx1 vector of variance factors
 x_sampl    = KxU matrix of estimated parameters
 S          = significance level of tests
 text       = string characterizing problem

 provides three tests
 - variance factor (if not only ones, eg for algebraic estimation)
 - covariance matrix
 - bias

 Wolfgang Förstner 2016-09-28
 wfoerstn@uni-bonn.de

0001 %% check estimation result, PCV Sect. 4.6.8
0002 %
0003 % check_estimation_result(R,x_true,S0,s0q_sample,x_sample,S,text)
0004 %
0005 % R          = redundancy (d.o.f. of variance factor)
0006 % x_true     = Ux1 vector of parameters
0007 % C0         = UxU reference covariance matrix of estimated parameters
0008 % s0q_sample = Kx1 vector of variance factors
0009 % x_sampl    = KxU matrix of estimated parameters
0010 % S          = significance level of tests
0011 % text       = string characterizing problem
0012 %
0013 % provides three tests
0014 % - variance factor (if not only ones, eg for algebraic estimation)
0015 % - covariance matrix
0016 % - bias
0017 %
0018 % Wolfgang Förstner 2016-09-28
0019 % wfoerstn@uni-bonn.de
0020 
0021 function f = check_estimation_result(R, x_true, C0, s0q_sample, x_sample, S, text)
0022 
0023 % number of parameters
0024 U = length(x_true);
0025 % sample size
0026 K = length(find(s0q_sample(:) > 0));
0027 
0028 disp('++++++++++++++++++++++++++++++++++++++++++++++++')
0029 disp(strcat('Checks for ', text))
0030 disp(['Number U of unknown parameters = ', num2str(U)]);
0031 disp(['Redundancy R                   = ', num2str(R)]);
0032 disp(['Number K of samples            = ', num2str(K)]);
0033 disp('-------------------------------------------------')
0034 
0035 %% mean and covariance matrix of sample
0036 me = mean(x_sample);
0037 dx = me'-x_true;
0038 Ce = cov(x_sample);
0039 
0040 %% check of sigma_0 (two-sided), PCV (4.354) ff
0041 if sum(s0q_sample) ~= K
0042     mean_s02 = mean(s0q_sample);
0043     F_tol_o = finv(1 - (1 - S) / 2, K * R, 100000);
0044     F_tol_u = finv((1 - S) / 2, K * R, 100000);
0045     
0046     if mean_s02 > F_tol_u && mean_s02 < F_tol_o
0047         disp(['variance factor s_0^2      ok: mean(s_0^2) = ', ...
0048             num2str(mean_s02, '% 12.4f'), '     in [', ...
0049             num2str(F_tol_u, '% 12.4f'), ',', num2str(F_tol_o, '% 12.4f'), ']']);
0050     else
0051         disp(['variance factor  s_0^2 not ok: mean(s_0^2) = ', ...
0052             num2str(mean_s02, '% 12.4f'), ' not in [', ...
0053             num2str(F_tol_u, '% 12.4f'), ',', num2str(F_tol_o, '% 12.4f'), '] *****']);
0054     end
0055     
0056     % plot histogram of variance factors with expected density (F-distibution)
0057     screensize = plot_init;
0058     f = figure('Color', 'w', 'Position', [20 20 screensize / 2]);
0059     hold on
0060     % number of bins = sqrt(K)
0061     N_bin = ceil(sqrt(K));
0062     % plot histogram
0063     hist(s0q_sample, N_bin);
0064     % store bin-centres r
0065     [~, r] = hist(s0q_sample, N_bin);
0066     % determine scale
0067     range = abs(r(N_bin) - r(1)) * N_bin / (N_bin - 1);
0068     % plot density with correct scale
0069     plot(r, N_bin * range * fpdf(r, R, 10000), '-r', 'LineWidth', 2);
0070     title(['Histogram of ', num2str(K), ' variance factors with reference density (', text, ')']);
0071 end
0072 
0073 %% Check of covariance matrix (two-sided), PCV (4.356) ff
0074 % Koch Parametersch Schätzung und Hypothesentests in linearen Modellen,
0075 % 2004, (287.5)
0076 %Ce,C0,K,S
0077 [lambda, T] = sugr_check_CovM(Ce, C0, K, S);
0078 
0079 if lambda > T(2) && lambda < T(1)
0080     disp(['covariance matrix C_xx     ok: lambda      = ', ...
0081         num2str(lambda, '% 12.4f'), '     in [', num2str(T(2), '% 12.4f'), ...
0082         ',', num2str(T(1), '% 12.4f'), ']']);
0083 else
0084     disp(['covariance matrix C_xx not ok: lambda      = ', ...
0085         num2str(lambda, '% 12.4f'), ' not in [', num2str(T(2), '% 12.4f'), ...
0086         ',', num2str(T(1), '% 12.4f'), '] *****']);
0087 end
0088 
0089 %% Check of mean parameters (two-sided), PCV (4.359) ff
0090 
0091 X_mean = K * dx' * inv(Ce) * dx;                                           %#ok<MINV>
0092 X_mean_tol_h = chi2inv(1 - (1 - S) / 2, U);
0093 X_mean_tol_l = chi2inv((1 - S) / 2, U);
0094 
0095 if X_mean < X_mean_tol_h && X_mean > X_mean_tol_l
0096     disp(['mean of parameters x       ok: mean(dx)    = ', ...
0097         num2str(X_mean, '% 12.4f'), '     in [', ...
0098         num2str(X_mean_tol_l, '% 12.4f'), ',', num2str(X_mean_tol_h, '% 12.4f'), ']']);
0099 else
0100     disp(['mean of parameters x   not ok: mean(dx)    = ', ...
0101         num2str(X_mean, '% 12.4f'), ' not in [', ...
0102         num2str(X_mean_tol_l, '% 12.4f'), ',', num2str(X_mean_tol_h, '% 12.4f'), ']', ' *****']);
0103 end
0104