% Jyvaskylä Summer School 2014 % Aku Seppanen % University of Eastern Finland, Kuopio % Department of Applied Physics %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example 3.1: line-fitting % %%%%%%%%%%%%%%%%%%%%%%%%%%%%% a_true = 1; b_true = 2; var_n = .05; % variance of noise n x = [0 1 2 3 4]'; n = sqrt(var_n)*randn(length(x),1); g = a_true*x + b_true + n; figure(1), clf plot(x,g,'k*') K = [x ones(size(x))]; f_LS = inv(K'*K)*K'*g; % LS-solution (could be written also as f_LS = K\g; a_LS = f_LS(1) b_LS = f_LS(2) g_LS = a_LS*x + b_LS; figure(1), hold on, plot(x,g_LS,'-') title(['a_{true}=',num2str(a_true),', b_{true}=',num2str(b_true),... ', a_{LS}=',num2str(a_LS),', b_{LS}=',num2str(b_LS)]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example 3.2: Fitting of a 2nd order polynomial % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a_true = 2; b_true = -5; c_true = 5 var_n = .05; % variance of noise n x = [0 1 2 3 4]'; n = sqrt(var_n)*randn(length(x),1); g = a_true*x.^2 + b_true*x + c_true + n; figure(1), clf plot(x,g,'k*') K = [x.^2 x ones(size(x))]; f_LS = inv(K'*K)*K'*g; % LS-solution (could be written also as f_LS = K\g; a_LS = f_LS(1) b_LS = f_LS(2) c_LS = f_LS(3) x_denser = linspace(min(x),max(x),100)'; g_LS = a_LS*x_denser.^2 + b_LS*x_denser + c_LS; figure(1), hold on, plot(x_denser,g_LS,'-') title(['a_{true}=',num2str(a_true),', b_{true}=',num2str(b_true),', c_{true}=',num2str(c_true),... ', a_{LS}=',num2str(a_LS),', b_{LS}=',num2str(b_LS), ', c_{LS}=',num2str(c_LS)]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example 3.3: line-fitting; only a single observation % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a_true = 1; b_true = 2; var_n = .0005; % variance of noise n x = [0]'; n = sqrt(var_n)*randn(length(x),1); g = a_true*x + b_true + n; K = [x ones(size(x))]; f_LS = inv(K'*K)*K'*g; % Matlab is not happy about this... rank(K) % rank of K is 1; underdetermined system of equations f_MN = pinv(K)*g; %minimum norm solution a_MN = f_MN(1) b_MN = f_MN(2) x_denser = [-1,0,1]'; g_MN = a_MN*x_denser + b_MN; figure(1), clf plot(x,g,'k*'), hold on, plot(x_denser,a_true*x_denser+b_true,'g-') plot(x_denser,g_MN,'b--') hold on title(['a_{true}=',num2str(a_true),', b_{true}=',num2str(b_true),... ', a_{MN}=',num2str(a_MN),', b_{MN}=',num2str(b_MN)]); legend('Observation','y_{true}','y_{MN}') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example 3.4: Gauss-Markov estimate; line-fitting % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a_true = 1; b_true = 2; var_n = [.05 .05 100 .05 .05]'; % third noise component has a very large variance x = [0 1 2 3 4]'; n = sqrt(var_n).*randn(length(x),1); g_true = a_true*x + b_true g = g_true + n; figure(1), clf plot(x,g,'k*') K = [x ones(size(x))]; % LS estimate f_LS = inv(K'*K)*K'*g; % LS-solution (could be written also as f_LS = K\g; a_LS = f_LS(1) b_LS = f_LS(2) % Gauss-Markov estimate Gamma_n = diag(var_n); W = inv(Gamma_n); f_GM = inv(K'*W*K)*K'*W*g; % LS-solution (could be written also as f_LS = K\g; a_GM = f_GM(1) b_GM = f_GM(2) g_LS = a_LS*x + b_LS; g_GM = a_GM*x + b_GM; figure(1), hold on, plot(x,g_true,'g-') plot(x,g_LS,'b--') plot(x,g_GM,'r--') title(['a_{true}=',num2str(a_true),', b_{true}=',num2str(b_true),... ', a_{LS}=',num2str(a_LS),', b_{LS}=',num2str(b_LS)... ', a_{GM}=',num2str(a_GM),', b_{GM}=',num2str(b_GM)]); legend('Observations','True line','LS-solution','Gauss-Markov solution') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example 3.5: Gauss-Markov estimate; line-fitting; correlated noise % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a_true = 1; b_true = 2.09; Gamma_n = [.0001 0 0 ;... 0 10.0001 10;... 0 10 10.0001]; % Observations 2 & 3 have large variances and are heavily correlated x = [0 1 2]'; n = chol(inv(Gamma_n))\randn(length(x),1); g_nonoise = a_true*x + b_true g = g_nonoise + n; figure(1), clf plot(x,g,'k*'), hold on plot(x,g_nonoise,'g-') K = [x ones(size(x))]; % LS estimate f_LS = inv(K'*K)*K'*g; % LS-solution (could be written also as f_LS = K\g; a_LS = f_LS(1) b_LS = f_LS(2) % Gauss-Markov estimate W = inv(Gamma_n); f_GM = inv(K'*W*K)*K'*W*g; % LS-solution (could be written also as f_LS = K\g; a_GM = f_GM(1) b_GM = f_GM(2) g_LS = a_LS*x + b_LS; g_GM = a_GM*x + b_GM; figure(1), hold on, plot(x,g_LS,'b--') plot(x,g_GM,'r--') title(['a_{true}=',num2str(a_true),', b_{true}=',num2str(b_true),... ', a_{LS}=',num2str(a_LS),', b_{LS}=',num2str(b_LS)... ', a_{GM}=',num2str(a_GM),', b_{GM}=',num2str(b_GM)]); legend('Observations','True g','LS-solution','Gauss-Markov solution') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example 3.6: Non-linear LS % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear a_true = 1; b_true = 5; x = linspace(0,1,20)'; var_n = .002; Gamma_n = var_n*eye(length(x)); W = inv(Gamma_n); n = sqrt(var_n)*randn(length(x),1); g = a_true*exp(-b_true*x) + n; x_dense = linspace(0,1,200)'; % dense grid for visualization g_dense = a_true*exp(-b_true*x_dense); figure(1), clf, plot(x_dense,g_dense,'g',x,g,'k+') a_init = 10; % initial guess for a b_init = .1; % initial guess for b steplength = .5; niter = 50; % max number of iterations f_LS = zeros(2,niter); f_LS(:,1) = [a_init; b_init]; % Gauss-Newton iteration for k = 1:niter-1 f_k = f_LS(:,k); h_f_k = f_k(1)*exp(-f_k(2)*x); J_k = [exp(-f_k(2)*x) -f_k(1)*x.*exp(-f_k(2)*x)]; f_LS(:,k+1) = f_k + steplength*inv(J_k'*W*J_k)*(J_k'*W*(g-h_f_k)); g_LS_dense = f_LS(1,k+1)*exp(-f_LS(2,k+1)*x_dense); figure(2), clf, plot(x_dense,g_dense,'g',x,g,'k+') hold on plot(x_dense,g_LS_dense,'r--') title(['Iteration k = ',num2str(k),', a = ',num2str(f_LS(1,k+1)),', b = ',num2str(f_LS(2,k+1))]) legend('True g','Observations','LS-solution') pause(.05) end