% gdama11
% W. Menke (2024), Geophysical Data Analysis and Inverse Theory
% with MATLAB(R) or Python, 5th Edition
% Code History:
% 2022-10-01, Imported from 4th Edition by W. Menke
% 2023-04-07, Edit and integration with text
% gdama11_01
%
% Least squares fit of straight line to d(z) and z(d).
clear all;
fprintf('gdama11_01\n');
% auxially variable z and data d
z = [1, 2, 3, 4]';
d = [1, 2, 3, 5]';
N=length(z);
% least squares fit to d(z)
M=2;
G=[ones(N,1), z];
mest = (G'*G)\(G'*d);
dpre = G*mest;
fprintf('d(z): intercept %f slope %f\n', mest(1), mest(2) );
% least squares fit to z(d)
G2=[ones(N,1), d];
mest2 = (G2'*G2)\(G2'*z);
dpre2 = G2*mest2;
% convert model parameters for z(d) case to d(z)
mest3=zeros(2,1);
mest3(1)=-mest2(1)/mest2(2);
mest3(2)=1/mest2(2);
dpre3=G*mest3;
fprintf('dp(zp): intercept %f slope %f\n', mest2(1), mest2(2) );
fprintf('zp(dp): intercept %f slope %f\n', mest3(1), mest3(2) );
% plot d(z) case
figure();
clf;
subplot(1,3,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 6, 0, 6] );
plot( z, d, 'ko', 'LineWidth', 3);
plot( z, dpre, 'r-', 'LineWidth', 3);
xlabel('z');
ylabel('d');
title('d(z)');
% plot z(d) case
subplot(1,3,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 6, 0, 6] );
plot( d, z, 'ko', 'LineWidth', 3);
plot( d, dpre2, 'g-', 'LineWidth', 3);
xlabel('d');
ylabel('z');
title('z(d)');
% plot z(d) and transformed z(d) cases
subplot(1,3,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 6, 0, 6] );
plot( z, d, 'ko', 'LineWidth', 3);
plot( z, dpre3, 'g-', 'LineWidth', 3);
plot( z, dpre, 'r-', 'LineWidth', 2);
xlabel('z');
ylabel('d');
title('d(z) and inv z(d)');
fprintf('Caption: Linear fits to z(d) to data (circles): (left) fits to d(z),\n');
fprintf('(middle) fit to z(d), (right) d(z) implied by fit to z(d)\n');
% gdama11_02
%
% example of transforming distributions
clear all;
fprintf('gdama11_02\n');
% model parameter, m
M=2000;
mmin = 0;
mmax = 1;
Dm = (mmax-mmin)/(M-1);
m = mmin + Dm*[0:M-1]';
% uniform distribution
Pm = ones(M,1);
% transformation mp = m^2
mp=m;
Pmp = 0.5 ./ sqrt(mp);
Pmp(1)=0.0; % set singularity to zero
% check
Am = Dm*sum(Pm);
Amp = Dm*sum(Pmp);
fprintf('areas: %.2f %.2f\n', Am, Amp );
figure();
clf;
% plot p(m)
subplot(1,2,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, 2] );
plot( m, Pm, 'r-', 'LineWidth', 3 );
xlabel('m');
ylabel('p(m)');
% plot p(m') where m'=m^2
subplot(1,2,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, 2] );
plot( mp(2:M), Pmp(2:M), 'b-', 'LineWidth', 3 );
xlabel('m''');
ylabel('p(m'')');
fprintf('Captiom: (left) Pdf p(m) is uniform, Pdf m(mp) where mp=m^2 is not uniform\n');
% gdama11_03
% An inverse problem to determine the model parameter
% in a model with an exponential function, d(z)=m1*exp(m2 z)
% using (Case 1) a linearizing transformation
% and (Case 2) Newton's Method.
clear all;
fprintf('gdama11_03\n');
% auxiallary variable z
N=20;
zmin = 0;
zmax = 1;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
% data, a exponential d=m1*exp(m2 z)
M=2;
mtrue = [0.9, -4]';
dtrue = mtrue(1)*exp(mtrue(2)*z);
% additive noise
sd=0.02;
bottom= 0.02;
dobs = dtrue + random('Normal',0,sd,N,1);
% but don't let data become negative
i=find(dobs<0);
dobs(i)=bottom;
% linearizing transformation, solve by standard least squares
% d = m1 exp( m2 z)
% ln(d) = ln(m1) + m2 * z
G=zeros(N,M);
G=[ones(N,1), z];
lndobs = log(dobs);
mestlin = (G'*G)\(G'*lndobs);
mestlin(1) = exp(mestlin(1));
dprelin = mestlin(1)*exp(mestlin(2)*z);
% Newton's Method iterative solution
% d = m1 exp( m2 z)
% dd/dm1 = exp( m2 z)
% dd/dm2 = m1 z exp( m2 z)
% initial guess is previous solution
mest = mestlin;
Nit=10; % maximum number of iterations
sdmmin = 1e-5; % terminate iteration early when increment dm
% has sqrt(variance(dm))<sdmmin
for it=(1:Nit)
% Newton's Method
dd = dobs - (mest(1)*exp(mest(2)*z));
G = [exp(mest(2)*z), mest(1)*z.*exp(mest(2)*z)];
dm = (G'*G)\(G'*dd);
mest = mest+dm;
sdm = sqrt((dm'*dm)/M); % sqrt(variance(dm))
if( sdm < sdmmin )
break; % terminate iterations early
end
end
% least squares prediction
dpre = mest(1)*exp(mest(2)*z);
figure();
clf;
% linear plot
subplot(1,2,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [zmin, zmax, 0, 1] );
plot( z, dobs, 'ro', 'LineWidth', 3 );
plot( z, dtrue, 'r-', 'LineWidth', 3 );
plot( z, dprelin, 'b-', 'LineWidth', 3 );
plot( z, dpre, 'g-', 'LineWidth', 3 );
xlabel('z');
ylabel('d');
% log plot
subplot(1,2,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [zmin, zmax, -3, 0] );
plot( z, log10(dobs), 'ro', 'LineWidth', 3 );
plot( z, log10(dtrue), 'r-', 'LineWidth', 3 );
plot( z, log10(dprelin), 'b-', 'LineWidth', 3 );
plot( z, log10(dpre), 'g-', 'LineWidth', 3 );
xlabel('z');
ylabel('log10(d)');
fprintf('Caption (left) Fit to exponentially decating data (circles), showing\n');
fprintf('true data (red), linearized fit (blue), iterative fit (green). (Right)\n');
fprintf('same as left, except plot is semi-logarithmic.\n');
fprintf('\n');
fprintf('True m1: %.3f m2: %.3f\n', mtrue(1), mtrue(2) );
fprintf('Newtons m1: %.3f m2: %.3f\n', mest(1), mest(2) );
fprintf('Linearizing m1: %.3f m2: %.3f\n', mestlin(1), mestlin(2) );
% gdama11_04
%
% 1D grid search for the one-parameter linear problem d=m1*z
% supports Figure 9.3
clear all;
fprintf('gdama11_04\n');
% auxiliary parameter z
N = 11;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
% only one model parameter, m1
M=1;
% linear model: d = m1*z
mtrue=2.5; % true model
dtrue=mtrue*z; % true data
% observed data is true data plus random noise
sd=2;
dobs=dtrue+random('Normal',0,sd,N,1);
% set up grid
Mg = 101;
mmin = 0;
mmax = 5;
Dm = (mmax-mmin)/(Mg-1);
m = mmin + Dm*[0:Mg-1];
% tabulate error E on a grid
E=zeros(Mg,1);
for i=(1:Mg)
dpre = m(i)*z;
e = dobs - dpre;
E(i) = e'*e;
end
% find point of minimum Error
[Emin, iEmin] = min(E);
mest=m(iEmin);
% plot Error and its minimum
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, max(E)] );
axis xy;
plot( m, E, 'k-', 'Linewidth', 3 );
plot( mest, Emin, 'ko', 'LineWidth', 3);
plot( [mest, mest], [0, max(E)/50], 'k-', 'LineWidth', 2);
xlabel('m');
ylabel('E');
fprintf('Caption: One dimensional grid search of error curve E(m) (black curve)\n');
fprintf('for minimum value (black circle)\n');
% gdama11_05
%
% E(m) for several hypothetial 1D non-linear problems
clear all;
fprintf('gdama11_05\n');
% auxiliary parameter z
N = 11;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
% set up 1D grid search
Mg = 501;
mmin = 0;
mmax = 5;
Dm = (mmax-mmin)/(Mg-1);
m = mmin + Dm*[0:Mg-1];
% only one model parameter, m1
M=1;
% model 1
m1true=2.5;
d1true = sin(5*(m1true-2.5)*z)-((m1true-2.5))*z;
sd=2;
d1obs=d1true+random('Normal',0,sd,N,1);
% tabulate error on the grid
E1=zeros(Mg,1);
for i=(1:Mg)
d1pre = sin(5*(m(i)-2.5)*z)-((m(i)-2.5))*z;
e = d1obs - d1pre;
E1(i) = (e'*e);
end
% find point of minimum error
[E1min, iE1min] = min(E1);
m1est=m(iE1min);
figure();
clf;
% plot E(m)
subplot(2,2,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, max(E1)] );
axis xy;
plot( m, E1, 'k-', 'Linewidth', 3 );
xlabel('m');
ylabel('E');
% model 2
m1true=2.5;
d1true = ((abs((m1true-2))-0.5)*z);
sd=2;
d1obs=d1true+random('Normal',0,sd,N,1);
% tabulate error on the grid
E1=zeros(Mg,1);
for i=[1:Mg]
d1pre = ((abs((m(i)-2))-0.5)*z);
e = d1obs - d1pre;
E1(i) = (e'*e);
end
% find point of minimum error
[E1min, iE1min] = min(E1);
m1est=m(iE1min);
% plot E(m)
subplot(2,2,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, max(E1)] );
axis xy;
plot( m, E1, 'k-', 'Linewidth', 3 );
xlabel('m');
ylabel('E');
% model 3
m1true=2.5;
d1true = sin(5*(m1true+z));
sd=1;
d1obs=d1true+random('Normal',0,sd,N,1);
% tabulate error on the grid
E1=zeros(Mg,1);
for i=[1:Mg]
d1pre = sin(5*(m(i)+z));
e = d1obs - d1pre;
E1(i) = (e'*e);
end
% find point of minumum error
[E1min, iE1min] = min(E1);
m1est=m(iE1min);
% plor E(m)
subplot(2,2,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, max(E1)] );
axis xy;
plot( m, E1, 'k-', 'Linewidth', 3 );
xlabel('m');
ylabel('E');
% model 4
m1true=2.5;
t = abs(m1true-2.0) + abs(m1true-3.0);
d1true = exp( -(t*z/10).^2 );
sd=1;
d1obs=d1true+random('Normal',0,sd,N,1);
% tabulate error on the grid
E1=zeros(Mg,1);
for i=(1:Mg)
t = abs(m(i)-2.0) + abs(m(i)-3.0);
d1pre = exp( -(t*z/10).^2 );
e = d1obs - d1pre;
E1(i) = (e'*e);
end
% find point of minimum error
[E1min, iE1min] = min(E1);
m1est=m(iE1min);
% plot E(m)
subplot(2,2,4);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, max(E1)] );
axis xy;
plot( m, E1, 'k-', 'Linewidth', 3 );
xlabel('m');
ylabel('E');
fprintf('Caption: Hypothetical error surfaces\n');
% gdama11_06
%
% E(m) for several 2D non-linear problems
clear all;
fprintf('gdama11_06\n');
% auxiliary parameter z
N = 11;
zmin = 0;
zmax = 1.0;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
za = z.^0.10;
zb = z.^0.20;
% 2D grid of model parameters m1 and m2
% variable m1
Mg1 = 101;
m1min = 0;
m1max = 1;
Dm1 = (m1max-m1min)/(Mg1-1);
m1 = m1min + Dm1*[0:Mg1-1];
% variable m2
Mg2 = 101;
m2min = 0;
m2max = 1;
Dm2 = (m2max-m2min)/(Mg2-1);
m2 = m2min + Dm2*[0:Mg2-1];
% Case 1
% true model
m1true = 0.4;
m2true = 0.6;
% true data
dtrue = ((m1true-0.2)*za-((m2true-0.3)*zb).^2).*z;
% observed data is true data plus random noise
sd=0.01;
dobs = dtrue + random('Normal',0,sd,N,1);
% tabulate error on the grid, keeping track of
% minimum error along the way
E=zeros(Mg1,Mg2);
first=1;
Emini=1;
Eminj=1;
for i=(1:Mg1)
for j=(1:Mg2)
dpre = ((m1(i)-0.2)*za-((m2(j)-0.3)*zb).^2).*z;
e = dobs - dpre;
E(i,j) = (e'*e);
if(first==1)
Emin=E(i,j);
Emini=i;
Eminj=j;
first=0;
dpresave = dpre;
elseif( E(i,j)<Emin )
Emin=E(i,j);
Emini=i;
Eminj=j;
dpresave = dpre;
end
end
end
fprintf('Emin %f at m1(%d) m2(%d) %f %f\n', Emin, Emini, Eminj, m1(Emini), m2(Eminj) );
% plot error surface E(m1,m2)
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [m2min, m2max, m1min, m1max] );
axis xy;
imagesc( [m2min, m2max], [m1min, m1max], E );
plot( m2(Eminj), m1(Emini), 'wo', 'LineWidth', 3 );
xlabel('m2');
ylabel('m1');
colorbar();
fprintf('Caption: Hypthetical error surface, (colors)\n');
fprintf('with minimum (white circle) superimposed.\n');
fprintf(' \n');
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
dmin = min(dobs);
dmax = max(dobs);
axis( [zmin, zmax, dmin, dmax] );
plot(z,dobs,'ko','LineWidth',2);
plot(z,dtrue,'k-','LineWidth',2);
plot(z,dpresave,'r-','LineWidth',2);
xlabel('z');
ylabel('d(z)');
colorbar();
fprintf('Caption: True data (black line), observed (circles) and predicted (red line)).\n');
fprintf('with minimum (white circle) superimposed.\n');
fprintf(' \n');
% Case 2
% true model parameters
m1true = 0.4;
m2true = 0.6;
% true data
dtrue = sin(12*pi*(m1true-0.5))*z+(m2true*z).^0.5;
% observed data is true data plus random noise
sd=0.1;
dobs = dtrue + random('Normal',0,sd,N,1);
% Tabulate error E(m1,m2) on the grid, searching
% for the point of minimum error along the way
E=zeros(Mg1,Mg2);
first=1;
for i=(1:Mg1)
for j=(1:Mg2)
dpre = sin(12*pi*(m1(i)-0.5))*z+(m2(j)*z).^0.5;
e = dobs - dpre;
E(i,j) = (e'*e);
if( first )
Emin=E(i,j);
Emini=i;
Eminj=j;
first=0;
dpresave = dpre;
elseif( E(i,j)<Emin )
Emin=E(i,j);
Emini=i;
Eminj=j;
dpresave = dpre;
end
end
end
fprintf('Emin %f at m1(%d) m2(%d) %f %f\n', Emin, Emini, Eminj, m1(Emini), m2(Eminj) );
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [m2min, m2max, m1min, m1max] );
axis xy;
imagesc( [m2min, m2max], [m1min, m1max], E );
plot( m2(Eminj), m1(Emini), 'wo', 'LineWidth', 3 );
xlabel('m2');
ylabel('m1');
colorbar();
fprintf('Caption: Hypthetical error surface, (colors)\n');
fprintf('with minimum (white circle) superimposed.\n');
fprintf(' \n');
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
dmin = min(dobs);
dmax = max(dobs);
axis( [zmin, zmax, dmin, dmax] );
plot(z,dobs,'ko','LineWidth',2);
plot(z,dtrue,'k-','LineWidth',2);
plot(z,dpresave,'r-','LineWidth',2);
xlabel('z');
ylabel('d(z)');
colorbar();
fprintf('Caption: True data (black line), observed (circles) and predicted (red line)).\n');
fprintf('with minimum (white circle) superimposed.\n');
fprintf(' \n');
% gdama11_07
% Grid Search method applied to inverse problem
% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)
clear all;
fprintf('gdama11_07\n');
% data are in a simple auxillary variable, x
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = Dx*[0:N-1]';
% two model parameters
M=2;
% true model parameters
mt = [1.21, 1.54]';
% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);
w0=20;
dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2); % true data
sd=0.4;
dobs = dtrue + random('Normal',0,sd,N,1); % observed data
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, xmax, 0, 4 ] );
plot(x,dtrue,'k-','LineWidth',3);
plot(x,dobs,'ko','LineWidth',3);
xlabel('x');
ylabel('d');
% Define 2D grid
L = 101;
Dm = 0.02;
m1min=0;
m2min=0;
m1a = m1min+Dm*[0:L-1]';
m2a = m2min+Dm*[0:L-1]';
m1max = m1a(L);
m2max = m2a(L);
% tabulate error E on grid.
% Note m1 varies along rows of E and m2 along columns
E = zeros(L,L);
for j = [1:L]
for k = [1:L]
dpre = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);
E(j,k) = (dobs-dpre)'*(dobs-dpre);
end
end
% Example of use of Esurface() function
% It refines the minimum using a quadratic approximation
% and provides an error estimate using the 2nd derivative
% note that the "column" variable is sent and returned
% before the "row" variable in this function
[ m2est, m1est, E0, cov21, status ] = gda_Esurface( m2a, m1a, E, sd^2 );
sigmam1 = sqrt( cov21(2,2) );
sigmam2 = sqrt( cov21(1,1) );
% evaluate prediction and plot it
dpre = sin(w0*m1est*x) + m1est*m2est;
plot(x,dpre,'r-','LineWidth',3);
fprintf('Caption: Exemplary curve fit wia grid search, showing true data (black curve)\n');
fprintf('observed data (circles) and predicted data (red curve)\n');
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [m2min, m2max, m1min, m1max] );
axis ij;
imagesc( [m2min, m2max], [m1min, m1max], E);
colorbar;
xlabel('m2');
ylabel('m1');
% plot true and estimate model
plot( m2est, m1est, 'wo', 'LineWidth', 3 );
plot( mt(2), mt(1), 'go', 'LineWidth', 3 );
% plot error bars
plot( [m2est, m2est]', [m1est-2*sigmam1, m1est+2*sigmam1]', 'r-', 'LineWidth', 2 );
plot( [m2est-2*sigmam2, m2est+2*sigmam2]', [m1est, m1est]', 'r-', 'LineWidth', 2 );
fprintf('Caption: Error surface, with true solution (green circle), estimate solution\n');
fprintf('(white circles) and 95 percent confidence intervals (red bars), superimposed\n');
% gdama11_08
%
% Newton's Method solution to exemplary 1D non-linear
% inverse problem, d(z) = sin(k*m*z)*z;
% This version uses a starting guess for the model
% parameter m1 that leads to convergence to the global
% minimum of the error surface (and hence to the correct
% value of the model parameter).
clear all;
fprintf('gdama11_08\n');
% z-axis
N = 21;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
zp = z.^0.5;
% define 1D grid for model parameter m1
Mg = 501;
mmin = 0;
mmax = 5;
Dm = (mmax-mmin)/(Mg-1);
m = mmin + Dm*[0:Mg-1];
% only one model parameter, m1
M=1;
% true model parameter
mtrue=3;
% true data
k=1;
dtrue = sin(k*mtrue.*zp).*z;
% observed data is true data plus random noise
sd=2;
dobs=dtrue+random('Normal',0,sd,N,1);
% tabulate error E(m) on the grid
E1=zeros(Mg,1);
for i=[1:Mg]
dpre = sin(k*m(i).*zp).*z;
e = dobs - dpre;
E1(i) = (e'*e);
end
% find global minimum
[E1min, iE1min] = min(E1);
m1est=m(iE1min);
% Newton's method
% derivative
% d = sin(k*m(i).*zp).*z;
% dd/dm = k.*zp.*z.*cos(k*m(i).*zp);
% d = d(m0) + (dd/dm)|m0 (m-m0)
% tabulate a quadratic version of the error
% correspinding to a linearized version of
% the inverse problem (linearized around a point m0)
E2=zeros(Mg,1);
m0 = 2.5;
im0 = floor((m0-mmin)/Dm)+1;
d0 = sin(k*m0.*zp).*z;
dddm = k.*zp.*z.*cos(k*m0.*zp);
for i=[1:Mg]
dpre = d0 + dddm * (m(i)-m0); % linearized approximation
e = dobs - dpre;
E2(i) = (e'*e);
end
% find minimum of error surface
[E2min, iE2min] = min(E2);
m2est=m(iE2min);
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
mminp=1; mmaxp=4;
axis( [mmin, mmax, 0, max(E1)] );
axis xy;
% true error surface, with circle at minimum, and line dropped to axis
plot( m, E1, 'k-', 'Linewidth', 3 );
plot( m1est, E1min, 'ko', 'Linewidth', 3 );
plot( [m1est, m1est], [0, E1min], 'k:', 'Linewidth', 2 );
% parabiolic approximation, with circle at minimum, and line dropped to
% axis
plot( m, E2, 'r:', 'Linewidth', 3 );
plot( m2est, E2min, 'ro', 'Linewidth', 3 );
plot( [m2est, m2est], [0, E2min], 'r:', 'Linewidth', 2 );
plot( m0, E2(im0), 'ko', 'Linewidth', 3 );
plot( [m0, m0], [0, E2(im0)], 'k:', 'Linewidth', 2 );
xlabel('m');
ylabel('E');
fprintf('Caption: Exemplary error surface E(m) (back curve) with minimum value at m=3 (black circle).');
fprintf('showing parabola tangegent to curve at m=2.5 (red curve), with minimum at m=3.3 (red circle)\n');
% gdama11_09
%
% Newton's Method solution to exemplary 1D non-linear
% inverse problem, d(z) = sin(k*m*z)*z
% This version uses a starting guess for the model
% parameter m1 that leads to convergence to a local
% minimum of the error surface (and hence to an incorrect
% value of the model parameter).
clear all;
fprintf('gdama11_09\n');
% auxiliary parameter z
N = 21;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
zp = z.^0.5;
% define 1D grid for model parameter m1
Mg = 501;
mmin = 0;
mmax = 5;
Dm = (mmax-mmin)/(Mg-1);
m = mmin + Dm*[0:Mg-1];
% only one model parameter, m1
M=1;
% true model parameter
mtrue=3;
% true data
k=1;
dtrue = sin(k*mtrue.*zp).*z;
% observed data is true data plus random noise
sd=2;
dobs=dtrue+random('Normal',0,sd,N,1);
% tabulate error E(m) on the grid
E1=zeros(Mg,1);
for i=[1:Mg]
dpre = sin(k*m(i).*zp).*z;
e = dobs - dpre;
E1(i) = (e'*e);
end
% find global minimum
[E1min, iE1min] = min(E1);
m1est=m(iE1min);
% Newton's method
% derivative
% d = sin(k*m(i).*zp).*z;
% dd/dm = k.*zp.*z.*cos(k*m(i).*zp);
% d = d(m0) + (dd/dm)|m0 (m-m0)
% tabulate a quadratic version of the error
% correspinding to a linearized version of
% the inverse problem (linearized around a point m0)E2=zeros(Mg,1);
m0 = 1.0;
im0 = floor((m0-mmin)/Dm)+1;
d0 = sin(k*m0.*zp).*z;
dddm = k.*zp.*z.*cos(k*m0.*zp);
for i=[1:Mg]
dpre = d0 + dddm * (m(i)-m0);
e = dobs - dpre;
E2(i) = (e'*e);
end
% find minimum of error surface
[E2min, iE2min] = min(E2);
m2est=m(iE2min);
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
mminp=1; mmaxp=4;
axis( [mmin, mmax, 0, max(E1)] );
axis xy;
% true error surface, with circle at minimum, and line dropped to axis
plot( m, E1, 'k-', 'Linewidth', 3 );
plot( m1est, E1min, 'ko', 'Linewidth', 3 );
plot( [m1est, m1est], [0, E1min], 'k:', 'Linewidth', 2 );
% parabiolic approximation, with circle at minimum, and line dropped to
% axis
plot( m, E2, 'r:', 'Linewidth', 3 );
plot( m2est, E2min, 'ro', 'Linewidth', 3 );
plot( [m2est, m2est], [0, E2min], 'r:', 'Linewidth', 2 );
plot( m0, E2(im0), 'ko', 'Linewidth', 3 );
plot( [m0, m0], [0, E2(im0)], 'k:', 'Linewidth', 2 );
xlabel('m');
ylabel('E');
fprintf('Caption: Exemplary error surface E(m) (back curve) with minimum value at m=3 (black circle).');
fprintf('showing parabola tangegent to curve at m=1 (red curve), with minimum at m=3.3 (red circle)\n');
% gdama11_10
% Newton's method applied to inverse problem
% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)
% supports Figure 9.9
clear all;
fprintf('gdama11_10\n');
% x-axis
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = Dx*[0:N-1]';
% true model parameters
M=2;
mt = [1.21, 1.54]';
% d(x)=g(x, m1, m2) with g=sin(w0*m(1)*x) + m(1)*m(2)
w0=20;
dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2); % true data
sd=0.4;
dobs = dtrue + random('Normal',0,sd,N,1); % observed data
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, xmax, 0, 4 ] );
plot(x,dtrue,'k-','LineWidth',3);
plot(x,dobs,'ko','LineWidth',3);
xlabel('x');
ylabel('d');
% 2D grid, for plotting purposes only
L = 101;
Dm = 0.02;
m1min=0;
m2min=0;
m1a = m1min+Dm*[0:L-1]';
m2a = m2min+Dm*[0:L-1]';
m1max = m1a(L);
m2max = m2a(L);
% tabulate error, E, on grid for plotting purposed only
E = zeros(L,L);
for j = (1:L)
for k = (1:L)
dpre = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);
E(j,k) = (dobs-dpre)'*(dobs-dpre);
end
end
% Newton's method, calculate derivatives
% y = sin( w0 m1 x) + m1 m2;
% dy/dm1 = w0 x cos( w0 m1 x) + m2
% dy/dm2 = m2
% initial guess and corresponding error
mg=[1.1,0.95]';
dg = sin(w0*mg(1)*x) + mg(1)*mg(2);
Eg = (dobs-dg)'*(dobs-dg);
plot( mg(2), mg(1), 'ko', 'LineWidth', 3 );
% save solution and minimum error as a function of iteration number
Niter=20;
Ehis=zeros(Niter,1);
m1his=zeros(Niter,1);
m2his=zeros(Niter,1);
Ehis(1)=Eg;
m1his(1)=mg(1);
m2his(1)=mg(2);
% iterate to improve initial guess
G = zeros(N,M);
for k = (2:Niter)
% predict data at current guess for solution
dg = sin( w0*mg(1)*x) + mg(1)*mg(2);
dd = dobs-dg;
Eg=dd'*dd; % current error
Ehis(k)=Eg; % save error history
% build linearized data kernel
G = zeros(N,2);
G(:,1) = w0 * x .* cos( w0 * mg(1) * x ) + mg(2);
G(:,2) = mg(2)*ones(N,1);
% least squares solution for improvement to m
dm = (G'*G)\(G'*dd);
% update solution
mg = mg+dm;
plot( mg(2), mg(1), 'wo', 'LineWidth', 2 );
% save solution history
m1his(k)=mg(1);
m2his(k)=mg(2);
end
% estimated solution is from final iteration
m1est = m1his(Niter);
m2est = m2his(Niter);
% evaluate prediction and plot it
dpre = sin(w0*m1est*x) + m1est*m2est;
plot(x,dpre,'r-','LineWidth',3);
fprintf('Caption: Exemplary curve fit with Newtons methiod, showing observed data (circles),\n');
fprintf('true data (black curve) and predicted data (red curve)\n')
% plot error surface
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [m2min, m2max, m1min, m1max] );
axis ij;
imagesc( [m2min, m2max], [m1min, m1max], E);
colorbar;
xlabel('m2');
ylabel('m1');
% plot true, history, starting and final solution on error surface
plot( m2his, m1his, 'w-', 'LineWidth', 1 );
plot( mt(2), mt(1), 'go', 'LineWidth', 2 );
plot( m2est, m1est, 'ro', 'LineWidth', 2 );
plot( m2his(1), m1his(1), 'wo', 'LineWidth', 2 );
fprintf('Caption: Error surface (colors). with initial guess of solution\n');
fprintf('(white circles). trajectory (white line). estimated solution (red\n');
fprintf('circle) and true solution (green)\n');
% plots of the history of E, m1 and m2
figure();
clf;
subplot(3,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [0:Niter-1], Ehis, 'k-', 'LineWidth', 2 );
xlabel('iteration');
ylabel('E');
subplot(3,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [0, Niter-1], [mt(1), mt(1)], 'r', 'LineWidth', 2 );
plot( [0:Niter-1], m1his, 'k-', 'LineWidth', 2 );
xlabel('iteration');
ylabel('m_1');
subplot(3,1,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [0, Niter-1], [mt(2), mt(2)], 'r', 'LineWidth', 2 );
plot( [0:Niter-1], m2his, 'k-', 'LineWidth', 2 );
xlabel('iteration');
ylabel('m_2');
fprintf('Caption: (top) Error as a function of iteration number. (middle) Model parameterr\n');
fprintf('m1 as a function of iteration number (black curve) and true value (red line). (bottom) Model parameter m2\n');
fprintf('as a fuction of iteration number (black curve) and true value (red line)\n');
% gdama11_11
%
% draws a simple Normal pdf p(x1,x2)
clear all;
fprintf('gdama11_11\n');
% x1-axis
Nx1 = 101;
x1min = 0;
x1max = 5.0;
Dx1 = (x1max-x1min)/(Nx1-1);
x1 = x1min + Dx1*[0:Nx1-1]';
% x2-axis
Nx2 = 101;
x2min = 0;
x2max = 5.0;
Dx2 = (x2max-x2min)/(Nx2-1);
x2 =x2min + Dx2*[0:Nx2-1]';
% pdf
P1=zeros(Nx1,Nx2);
x1bar = 2.25;
x2bar = 2.08;
bar = [x1bar, x2bar]';
sx1 = 0.5;
sx2 = 1.0;
C1 = diag( [sx1^2, sx2^2]' );
CI1 = inv(C1);
DC1 = det(C1);
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
for i=(1:Nx1)
for j=(1:Nx2)
x =[x1(i), x2(j)]' - bar;
P1(i,j) = norm1*exp( -0.5 * x'*CI1*x );
end
end
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [x1min, x1max, x2min, x2max] );
axis ij;
imagesc( [x1min, x1max], [x2min, x2max], P1 );
plot( x2bar, x1bar, 'wo', 'LineWidth', 3 );
plot( [x2bar, x2bar], [x1min, x1min+0.15], 'w-', 'LineWidth', 3 );
plot( [x2min, x2min+0.15], [x1bar, x1bar], 'w-', 'LineWidth', 3 );
xlabel('x2');
ylabel('x1');
colorbar();
fprintf('Caption: Exemplary Normal pdf p(x1,x2) with mean (white circle) superimposed\n');
% gdama11_12
%
% Plots a simple Normal p.d.f. p(x1,x1) with a parameteric
% curve superimposed and then samples the p.d.f. under the
% curve and makes a second plot of p.d.f. as a function of
% arc-length s along the curve. In this version, the curve
% is very smooth and p(s) is bell-shaped.
% 2D grid in (x1,x2)
clear all;
fprintf('gdama11_12\n');
% x1-axis
Nx1 = 101;
x1min = 0;
x1max = 5.0;
Dx1 = (x1max-x1min)/(Nx1-1);
x1 = x1min + Dx1*[0:Nx1-1]';
% x2-axis
Nx2 = 101;
x2min = 0;
x2max = 5.0;
Dx2 = (x2max-x2min)/(Nx2-1);
x2 =x2min + Dx2*[0:Nx2-1]';
% parameters for a Normal pdf
P1=zeros(Nx1,Nx2);
x1bar = 2.25;
x2bar = 2.08;
bar = [x1bar, x2bar]';
sx1 = 0.5;
sx2 = 1.0;
C1 = diag( [sx1^2, sx2^2]' );
CI1 = inv(C1);
DC1 = det(C1);
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% tabulate the Normal pdf on the grid
for i=(1:Nx1)
for j=(1:Nx2)
x =[x1(i), x2(j)]' - bar;
P1(i,j) = norm1*exp( -0.5 * x'*CI1*x );
end
end
% axis for parametric curve
Ns = 51;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% parameric curve
x2p = 1+s+sin(2*pi*s/smax)-2*(s/smax).^2;
x1p = s;
% pdf along parameteric curve
ipx1 = 1+floor((x1p-x1min)/Dx1);
ipx2 = 1+floor((x2p-x2min)/Dx2);
% insure indices in range
i=find(ipx1<1);
ipx1(i)=1;
i=find(ipx1>Nx1);
ipx1(i)=Nx1;
i=find(ipx2<1);
ipx2(i)=1;
i=find(ipx2>Nx2);
ipx2(i)=Nx2;
Ps=zeros(Ns,1);
% evaluate P at indices
for i = (1:Ns)
Ps(i) = P1(ipx1(i), ipx2(i));
end
% maximum along curve
[Pmax, ismax]=max(Ps);
x1smax = x1p(ismax);
x2smax = x2p(ismax);
smaxP = x1smax;
% plot of 2D pdf
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [x1min, x1max, x2min, x2max] );
axis ij;
imagesc( [x1min, x1max], [x2min, x2max], P1 );
plot( x2bar, x1bar, 'wo', 'LineWidth', 3 );
plot( [x2bar, x2bar], [x1min, x1min+0.1], 'w-', 'LineWidth', 3 );
plot( [x2min, x2min+0.1], [x1bar, x1bar], 'w-', 'LineWidth', 3 );
plot( x2p, x1p, 'w-', 'LineWidth', 3 );
plot( x2smax, x1smax, 'ko', 'LineWidth', 4);
plot( [x2min, x2smax], [x1smax, x1smax], 'w:', 'LineWidth', 2);
plot( [x2smax, x2smax], [x1min, x1smax], 'w:', 'LineWidth', 2);
xlabel('x2');
ylabel('x1');
colorbar();
fprintf('Caption: Exemplary normal pdf (colors) with mean (white circle),\n');
fprintf('curve f(x)=0 (white curve), and maximum likelihood point (black circle), superimposed.')
% plot of 1D pdf (along curve)
figure(2);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [smin, smax, 0, 0.3] );
axis xy;
plot( s, Ps, 'b-', 'Linewidth', 3 );
xlabel('s');
ylabel('p(s)');
plot( [smaxP, smaxP], [0, 0.03], 'k-', 'LineWidth', 3);
plot( smaxP, Pmax, 'bo', 'LineWidth', 3);
% calculate mean of distribution
norm = Ds*sum(Ps);
Ps = Ps/norm;
smean = Ds*sum(s.*Ps);
plot( [smean, smean], [0, 0.02], 'k-', 'LineWidth', 3);
fprintf('Caption: probability alone the curve p(x)=0, with the maximum\n');
fprintf('likelihood point shown (circle and long bar) and mean (short bar).');
% gdama11_13
%
% Plots a simple Normal p.d.f. p(x1,x1) with a parameteric
% curve superimposed and then samples the p.d.f. under the
% curve and makes a second plot of p.d.f. as a function of
% arc-length s along the curve. In this version, the curve
% is complicated and p(s) is bimodal.
clear all;
fprintf('gdama11_13\n');
% x1-axis
Nx1 = 101;
x1min = 0;
x1max = 5.0;
Dx1 = (x1max-x1min)/(Nx1-1);
x1 = x1min + Dx1*[0:Nx1-1]';
% x2-axis
Nx2 = 101;
x2min = 0;
x2max = 5.0;
Dx2 = (x2max-x2min)/(Nx2-1);
x2 =x2min + Dx2*[0:Nx2-1]';
% paameters for a Normal pdf p(x1,x2)
P1=zeros(Nx1,Nx2);
x1bar = 2.25;
x2bar = 2.08;
bar = [x1bar, x2bar]';
sx1 = 0.5;
sx2 = 1.0;
C1 = diag( [sx1^2, sx2^2]' );
CI1 = inv(C1);
DC1 = det(C1);
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% tabulate the Normal pdf on the grid
for i=(1:Nx1)
for j=(1:Nx2)
x =[x1(i), x2(j)]' - bar;
P1(i,j) = norm1*exp( -0.5 * x'*CI1*x );
end
end
% axis for parametric curve
Ns = 101;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% parameric curve (carefully tuned to be aestetically complciated)
x2p = 1+s+sin(2*pi*s/smax)-(2*(s/smax).^2)+0.3*exp(-0.5*((s-2.3).^2)/(0.1^2))-0.3*exp(-0.5*((s-2.8).^2)/(0.1^2));
x1p = s;
% pdf on parameteric curve
ipx1 = 1+floor((x1p-x1min)/Dx1);
ipx2 = 1+floor((x2p-x2min)/Dx2);
% insure indices in range
i=find(ipx1<1);
ipx1(i)=1;
i=find(ipx1>Nx1);
ipx1(i)=Nx1;
i=find(ipx2<1);
ipx2(i)=1;
i=find(ipx2>Nx2);
ipx2(i)=Nx2;
Ps=zeros(Ns,1);
% evaluate P at indices
for i = [1:Ns]
Ps(i) = P1(ipx1(i), ipx2(i));
end
% maximum along curve
[Pmax, ismax]=max(Ps);
x1smax = x1p(ismax);
x2smax = x2p(ismax);
smaxP = x1smax;
% 2D plot of p(x1,x2) with curve superimposed
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [x1min, x1max, x2min, x2max] );
axis ij;
imagesc( [x1min, x1max], [x2min, x2max], P1 );
plot( x2bar, x1bar, 'wo', 'LineWidth', 3 );
plot( [x2bar, x2bar], [x1min, x1min+0.1], 'w-', 'LineWidth', 3 );
plot( [x2min, x2min+0.1], [x1bar, x1bar], 'w-', 'LineWidth', 3 );
plot( x2p, x1p, 'w-', 'LineWidth', 3 );
plot( x2smax, x1smax, 'ko', 'LineWidth', 4);
plot( [x2min, x2smax], [x1smax, x1smax], 'w:', 'LineWidth', 2);
plot( [x2smax, x2smax], [x1min, x1smax], 'w:', 'LineWidth', 2);
xlabel('x2');
ylabel('x1');
colorbar();
fprintf('Caption: Exemplary normal pdf (colors) with mean (white circle),\n');
fprintf('curve f(x)=0 (white curve), and maximum likelihood point (black circle), superimposed.')
% 1D plot of p(s)
figure(2);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [smin, smax, 0, 0.3] );
axis xy;
plot( s, Ps, 'b-', 'Linewidth', 3 );
xlabel('s');
ylabel('p(s)');
plot( [smaxP, smaxP], [0, 0.03], 'k-', 'LineWidth', 3);
plot( smaxP, Pmax, 'bo', 'LineWidth', 3);
% calculate mean of distribution
norm = Ds*sum(Ps);
Ps = Ps/norm;
smean = Ds*sum(s.*Ps);
plot( [smean, smean], [0, 0.02], 'k-', 'LineWidth', 3);
fprintf('Caption: probability alonr the curve p(x)=0, with the maximum\n');
fprintf('likelihood point shown (circle and long bar) and mean (short bar).\n');
% gdama11_14
% Gradient method applied to inverse problem
% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)
clear all;
fprintf('gdama11_14\n');
% Auxially variable x
N=40;
x = [0:N-1]'/N;
% true model parameters
mt = [1.21, 1.54]';
% true data
w0=20;
dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2);
% observed data
sd=0.4;
dobs = dtrue + random('Normal',0,sd,N,1);
% 2D (m1,m2) grid
M = 101;
dm = 0.02;
m1min=0;
m2min=0;
m1a = m1min+dm*[0:M-1]';
m2a = m2min+dm*[0:M-1]';
m1max = m1a(M);
m2max = m2a(M);
% tabulate E(m1,m2) on grid (for plotting purposes only)
E = zeros(M,M);
for j = (1:M)
for k = (1:M)
yest = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);
E(j,k) = (dobs-yest)'*(dobs-yest);
end
end
% parameters for gradient method
alpha = 0.05;
c1 = 0.0001;
c2 = 0.9;
tau = 0.5;
% trial solution
mgo=[1,1]';
% save history of Error and model parameters
Niter=500;
Ehis=zeros(Niter+1,1);
m1his=zeros(Niter+1,1);
m2his=zeros(Niter+1,1);
% error and its gradient at the trial solution
ygo = sin( w0*mgo(1)*x) + mgo(1)*mgo(2);
Ego = (ygo-dobs)'*(ygo-dobs);
dydmo = zeros(N,2);
dydmo(:,1) = w0 * x .* cos( w0 * mgo(1) * x ) + mgo(2);
dydmo(:,2) = mgo(2)*ones(N,1);
dEdmo = 2*dydmo'*(ygo-dobs);
Nhis=0;
for k = (1:Niter)
% downhill direction
v = -dEdmo / sqrt(dEdmo'*dEdmo);
% save history
if( k==1 )
Nhis=Nhis+1; % counts actual iterations
Ehis(Nhis)=Ego;
m1his(Nhis)=mgo(1);
m2his(Nhis)=mgo(2);
end
% backstep
for kk=(1:10)
mg = mgo+alpha*v;
yg = sin( w0*mg(1)*x) + mg(1)*mg(2);
Eg = (yg-dobs)'*(yg-dobs);
dydm = zeros(N,2);
dydm(:,1) = w0 * x .* cos( w0 * mg(1) * x ) + mg(2);
dydm(:,2) = mg(2)*ones(N,1);
dEdm = 2*dydm'*(yg-dobs);
if( (Eg<=(Ego + c1*alpha*v'*dEdmo)) )
break;
end
alpha = tau*alpha;
end
% change in solution
Dmg = sqrt( (mg-mgo)'*(mg-mgo) );
% update
mgo=mg;
ygo = yg;
Ego = Eg;
dydmo = dydm;
dEdmo = dEdm;
% save history
Nhis=Nhis+1;
Ehis(Nhis)=Eg;
m1his(Nhis)=mg(1);
m2his(Nhis)=mg(2);
if( Dmg < 1.0e-6 )
break; % terminate iterations when change
% in solution is suffiently small
end
end
% estimated solution is last iteration
m1est = m1his(Nhis);
m2est = m2his(Nhis);
% truncate to actual lenth
Ehis = Ehis(1:Nhis);
m1his=m1his(1:Nhis);
m2his=m2his(1:Nhis);
% evaluate prediction and plot it
dpre = sin(w0*m1est*x) + m1est*m2est;
% plot true and observed data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 1, 0, 3] );
plot(x,dobs,'ko','LineWidth',3);
plot(x,dtrue,'k-','LineWidth',3);
xlabel('x');
ylabel('d(x)');
hold on;
plot(x,dpre,'r-','LineWidth',3);
fprintf('Caption: True data (black curve), observed data (black circles) and predicted data (red curve).\n');
% plot error surface
figure();
clf;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
xlabel('m2');
ylabel('m1');
axis( [m2min, m2max, m1min, m1max] );
axis ij;
imagesc( [m2min, m2max], [m1min, m1max], E);
colorbar;
% plot true, history, starting and final solution on error surface
plot( m2his, m1his, 'w-', 'LineWidth', 1 );
plot( mt(2), mt(1), 'go', 'LineWidth', 2 );
plot( m2est, m1est, 'ro', 'LineWidth', 2 );
plot( m2his(1), m1his(1), 'wo', 'LineWidth', 2 )
fprintf('Caption: Error surface E(m1,m2) (colors), with initial guess (white circle)\n');
fprintf('minimization path (white line), estimated solution (red circle)\n');
fprintf('and true solution (green circle), superimposed.\n');
% plot history of E, m1, m2
figure();
clf;
subplot(3,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [0:Nhis-1], Ehis, 'k-', 'LineWidth', 3 );
xlabel('iteration i');
ylabel('E(i)');
subplot(3,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [0:Nhis-1], m1his, 'k-', 'LineWidth', 3 );
plot( [0,Nhis-1], [mt(1), mt(1)], 'r:', 'LineWidth', 2);
xlabel('iteration i');
ylabel('m1(i)');
subplot(3,1,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [0:Nhis-1], m2his, 'k-', 'LineWidth', 3 );
plot( [0,Nhis-1], [mt(2), mt(2)], 'r:', 'LineWidth',2 );
xlabel('iteration i');
ylabel('m2(i)');
fprintf('Caption: (top) Error E(i) as a function of iteration number i.\n');
fprintf('(middle) Model parameter m1(i) as a function of iteration number i,\n');
fprintf('with true solution (red line) for comparison. (bottom) Same as\n');
fprintf('middel, except for model parameter m2\n');
% gdama11_15
% plots the distributions of one bit-mutations
% for both normal binary and Gray code representations
% of three exemplary 16-bit integers
clear all;
fprintf('gdama11_15\n');
global gray1 gray1index pow2;
load('gda_gray1table.mat');
figure(1);
clf;
subplot(3,2,1);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
axis( [0, 65535, 0, 1] );
i1 = 11111;
plot( [i1,i1]', [0, 0.2]', 'r-', 'LineWidth', 2 );
plot( [i1,i1]', [0.8, 1]', 'r-', 'LineWidth', 2 );
g1=gda_int2bin(i1);
for i=(1:16)
g2 = g1;
k=i;
if( g1(k) == '0' )
g2(k) = '1';
else
g2(k) = '0';
end
i2=gda_bin2int(g2);
plot( [i2,i2]', [0.2, 0.8]', 'k-', 'LineWidth', 2 );
end
subplot(3,2,3);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
axis( [0, 65535, 0, 1] );
i1 = 33333;
plot( [i1,i1]', [0, 0.2]', 'r-', 'LineWidth', 2 );
plot( [i1,i1]', [0.8, 1]', 'r-', 'LineWidth', 2 );
g1=gda_int2bin(i1);
for i=(1:16)
g2 = g1;
k=i;
if( g1(k) == '0' )
g2(k) = '1';
else
g2(k) = '0';
end
i2=gda_bin2int(g2);
plot( [i2,i2]', [0.2, 0.8]', 'k-', 'LineWidth', 2 );
end
subplot(3,2,5);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
axis( [0, 65535, 0, 1] );
i1 = 55555;
plot( [i1,i1]', [0, 0.2]', 'r-', 'LineWidth', 2 );
plot( [i1,i1]', [0.8, 1]', 'r-', 'LineWidth', 2 );
g1=gda_int2bin(i1);
for i=(1:16)
g2 = g1;
k=i;
if( g1(k) == '0' )
g2(k) = '1';
else
g2(k) = '0';
end
i2=gda_bin2int(g2);
plot( [i2,i2]', [0.2, 0.8]', 'k-', 'LineWidth', 2 );
end
subplot(3,2,2);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
axis( [0, 65535, 0, 1] );
i1 = 11111;
plot( [i1,i1]', [0, 0.2]', 'r-', 'LineWidth', 2 );
plot( [i1,i1]', [0.8, 1]', 'r-', 'LineWidth', 2 );
g1=gda_int2gray(i1);
for i=(1:16)
g2 = g1;
k=i;
if( g1(k) == '0' )
g2(k) = '1';
else
g2(k) = '0';
end
i2=gda_gray2int(g2);
plot( [i2,i2]', [0.2, 0.8]', 'k-', 'LineWidth', 2 );
end
subplot(3,2,4);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
axis( [0, 65535, 0, 1] );
i1 = 33333;
plot( [i1,i1]', [0, 0.2]', 'r-', 'LineWidth', 2 );
plot( [i1,i1]', [0.8, 1]', 'r-', 'LineWidth', 2 );
g1=gda_int2gray(i1);
for i=[1:16]
g2 = g1;
k=i;
if( g1(k) == '0' )
g2(k) = '1';
else
g2(k) = '0';
end
i2=gda_gray2int(g2);
plot( [i2,i2]', [0.2, 0.8]', 'k-', 'LineWidth', 2 );
end
subplot(3,2,6);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
axis( [0, 65535, 0, 1] );
i1 = 55555;
plot( [i1,i1]', [0, 0.2]', 'r-', 'LineWidth', 2 );
plot( [i1,i1]', [0.8, 1]', 'r-', 'LineWidth', 2 );
g1=gda_int2gray(i1);
for i=(1:16)
g2 = g1;
k=i;
if( g1(k) == '0' )
g2(k) = '1';
else
g2(k) = '0';
end
i2=gda_gray2int(g2);
plot( [i2,i2]', [0.2, 0.8]', 'k-', 'LineWidth', 2 );
end
fprintf('Caption: Distribution of integers (black bars) resulting from a bit bit mutation\n');
fprintf('of a target inter (red bar). (Left column) Standard coding, (right colummn) Gray coding.\n');
% gdama11_16
% Genetic Algorithm applied to inverse problem
% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)
clear all;
fprintf('gdama11_16\n');
global gray1 gray1index pow2;
load('gray1table.mat'); % load table of Gray codes
% two model parameters
M=2;
% Number of generations
L=100;
Esave = zeros(L,1);
Msave = zeros(L, M);
% x-axis
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = Dx*[0:N-1]';
% true model parameters
mt = [1.21, 1.54]';
% y=f(x, m1, m2);
w0=20;
dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2);
sd=0.4;
dobs = dtrue + random('Normal',0,sd,N,1);
% plot data
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, xmax, 0, 4 ] );
plot(x,dtrue,'k-','LineWidth',3);
plot(x,dobs,'ko','LineWidth',3);
xlabel('x');
ylabel('d');
% define 2D grid
% this is just for plotting purposes
% the genetic algorith does not require a grid
Lm = 101;
Dm = 0.02;
m1min=0;
m2min=0;
m1a = m1min+Dm*[0:Lm-1]';
m2a = m2min+Dm*[0:Lm-1]';
m1max = m1a(Lm);
m2max = m2a(Lm);
% compute error E on the grid
% this is just for plotting purposes
% the genetic algorith does not require a grid
E = zeros(Lm,Lm);
for j = (1:Lm)
for k = (1:Lm)
dpre = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);
E(j,k) = (dobs-dpre)'*(dobs-dpre);
end
end
K = 100; % number of individuals in the population
Kold = K;
m0 = [0.25, 0.30]; % initial guess for m1
mL = [m1min, m2min]'; % solution > mL
mR = [m1max, m2max]'; % solution < mR
dpre = sin(w0*m0(1)*x) + m0(1)*m0(2); % initial prediction
e = dobs-dpre; % individual errors
E0 = e'*e; % total error at start
mind = zeros(K, M); % solution for each individual
Eind = E0*ones(K,1); % error for each individual
F = ones(K,1); % fitness of each individual
genes = char(48*ones(K, 16*M)); % genes for each individual
for i=(1:M)
r = random('Normal',0,0.05,K,1);
mind(:,i) = m0(i)*ones(K,1)+r; % all individuals initialized to initial guess
end
for i=(1:K)
for j=(1:M)
k = floor(65535.999*(mind(i,j)-mL(j))/(mR(j)-mL(j)));
k1 = (j-1)*16+1;
k2 = k1+15;
genes(i,k1:k2) = gda_int2gray(k);
end
end
% save initial solution
mindsave = mind;
for gen=(1:L)
% each individual replicates
mind = [mind;mind];
genes = [genes;genes];
Eind = [Eind;Eind];
K = 2*K;
% mutate genes
for i=(1:K)
if( unidrnd(2,1) == 1 ) % mutate only half of individuals
continue;
end
k = unidrnd(16*M,1); % one random mutation
if( genes(i,k)=='0' )
genes(i,k)='1';
else
genes(i,k)='0';
end
end
% conjugate genes; only one conjugation per generation
j1 = unidrnd(K,1); % random individual
j2 = unidrnd(K,1); % another random individual
k = unidrnd(M*16-2,1)+1; % 2<=k<=16M-1
g1(1,1:16*M) = genes(j1,:);
g1(1,k+1:16*M) = genes(j2,k+1:16*M);
g2(1,1:16*M) = genes(j2,:);
g2(1,k+1:16*M) = genes(j1,k+1:16*M);
genes(j1,:) = g1;
genes(j2,:) = g2;
% update solution and error from genes
for i=(1:K)
for j=(1:M)
k1 = (j-1)*16+1;
k2 = k1+15;
k = gda_gray2int(genes(i,k1:k2));
mind(i,j) = (mR(j)-mL(j))*(k/65535.999)+mL(j);
end
end
for i=(1:K)
dpre = sin(w0*mind(i,1)*x) + mind(i,1)*mind(i,2);
e = dobs-dpre;
Eind(i) = e'*e;
end
Emax = max(Eind);
% survival of the fitness
c=5;
F = exp(-c*Eind/Emax).*random('uniform',0,1,K,1);
[F_sorted, k] = sort(F);
k = k(Kold+1:2*Kold,1);
K = Kold;
mind = mind(k,1:M);
Eind = Eind(k,1);
genes = genes(k,1:16*M);
[Emin, k] = min(Eind);
Esave(gen,1) = Emin;
msave(gen,1:M)= mind(k,1:M);
end
% solution is lowest error individual in last generation
m1est = msave(L,1);
m2est = msave(L,2);
% evaluate prediction and plot it
dpre = sin(w0*m1est*x) + m1est*m2est;
plot(x,dpre,'r-','LineWidth',3);
fprintf('Caption: True data (black curve), observed data (black circles)\n');
fprintf('predicted data (reed curve)\n');
% plot error surface
figure();
clf;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [m2min, m2max, m1min, m1max] );
axis ij;
imagesc( [m2min, m2max], [m1min, m1max], E);
colorbar;
xlabel('m2');
ylabel('m1');
% plor solution
plot( mindsave(:,2), mindsave(:,1), 'k.', 'LineWidth', 2 );
plot( mind(:,2), mind(:,1), 'w.', 'LineWidth', 2 );
plot( mt(2), mt(1), 'go', 'LineWidth', 3 );
plot( m2est, m1est, 'ro', 'LineWidth', 3 );
fprintf('Caption: Errir surface (colors), with initial population (black dots)\n');
fprintf('final population (white dots), mean of final population (red circle) and\n');
fprintf('true solution (green circle)\n');
% figure 3 is the solution progress
figure(3);
clf;
subplot(3,1,1); % error
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, L, 0, 1.1*E0] );
plot( [0, L]', [E0, E0]', 'r-', 'LineWidth', 2);
plot( [1:L]', Esave, 'k-', 'LineWidth', 2);
xlabel('generation');
ylabel('E');
subplot(3,1,2); % m1
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, L, 0, 2] );
plot( [0, L]', [mt(1), mt(1)]', 'r-', 'LineWidth', 2);
plot( [1:L]', msave(:,1)', 'k-', 'LineWidth', 2);
xlabel('generation');
ylabel('m1');
subplot(3,1,3); % m2
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, L, 0, 2] );
plot( [0, L]', [mt(2), mt(2)]', 'r-', 'LineWidth', 2);
plot( [1:L]', msave(:,2)', 'k-', 'LineWidth', 2);
xlabel('generation');
ylabel('m2');
fprintf('Caption: (top) Error as a function of generation number (black curve), with\n');
fprintf('inital error (red line). (middle) Model parametre m1 as a function of\n');
fprintf('generation number (black curve), with true value shown form comparison (red line).\n');
fprintf('(bottom) Same as middle except for model parameter m2.\n');
% gdama11_17
% Genetic algorithm applied to an "approximately
% separable inverse problem" where the M=20 model parameters
% have a natural ordering and the data are localized
% averages of the model parameters.
% This code compares error of two cases, A and B,
% where the probability of mutation and recombination
% differ between the cases. Each case has 5 trials.
clear all;
fprintf('gdama11_17\n');
global gray1 gray1index pow2;
load('gray1table.mat');
fprintf('gdama11_16\n');
fprintf('This calculation is pretty slow, so progress shown (ends at B5)\n');
% model parameters
M=20;
% Number of generations
L=50;
Esave = zeros(L,1);
Msave = zeros(L, M);
% plot error
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1, L, 0, 1.5 ] );
xlabel('generation');
ylabel('E');
Ntrials=5;
% part 1 --------------------------------------------
Nrecom = 0;
imutate = 2;
for itrial=[1:Ntrials]
% number of recombinations per generation
% fitness coefficient
c=5;
% true model parameters are linear ramp
mt = 0.5+[0:M-1]'/(M-1);
% data average of 3 neighboring model parameters
N=M;
x = [1:M]';
dtrue = (mt+[0;mt(1:M-1)]+[0;0;mt(1:M-2)])/3;
sd=0.025;
dobs = dtrue + random('Normal',0,sd,N,1);
K = 100; % number of individuals in the population
Kold = K;
m0 = ones(M,1); % initial guess for m1
mL = zeros(M,1); % solution > mL
mR = 2*ones(M,1)'; % solution < mR
mind = zeros(K, M); % solution for each individual
Eind = zeros(K,1); % error for each individual
Ftind = zeros(K,1); % fitness of each individual
genes = char(48*ones(K, 16*M)); % genes for each individual
% all individuals initialized to initial guess + random number
for i=(1:M)
r = random('Normal',0,0.025,K,1);
mind(:,i) = m0(i)*ones(K,1)+r;
end
% compute initial error
for i=(1:K)
dpre = (mind(i,:)'+[0;mind(i,1:M-1)']+[0;0;mind(i,1:M-2)'])/3; % predicetd data
e = dobs-dpre;
Eind(i) = e'*e;
end
Emax = max(Eind);
Emaxstart = Emax;
Ftind = exp(-c*Eind/Emax).*random('uniform',0,1,K,1); % initial fitness
for i=(1:K) % convert model parameters to genes
for j=(1:M)
k = floor(65535.999*(mind(i,j)-mL(j))/(mR(j)-mL(j)));
k1 = (j-1)*16+1;
k2 = k1+15;
genes(i,k1:k2) = gda_int2gray(k);
end
end
for gen=(1:L) % loop over generations
% each individual replicates
mind = [mind;mind];
genes = [genes;genes];
Eind = [Eind;Eind];
K = 2*K;
% mutate genes
for i=(1:K)
if( unidrnd(imutate,1) == 1 ) % mutate 1/imutate of individuals
continue;
end
k = unidrnd(16*M,1); % one random mutation
if( genes(i,k)=='0' )
genes(i,k)='1';
else
genes(i,k)='0';
end
end
% recominate genes; Nrecom per generation
for rc = (1:Nrecom)
j1 = unidrnd(K,1); % random individual
j2 = unidrnd(K,1); % another random individual
k = 16*(unidrnd(M-2,1)+1)+1; % always at start of word
g1(1,1:16*M) = genes(j1,:);
g1(1,k:16*M) = genes(j2,k:16*M);
g2(1,1:16*M) = genes(j2,:);
g2(1,k:16*M) = genes(j1,k:16*M);
genes(j1,:) = g1;
genes(j2,:) = g2;
end
% update solution and error from genes
for i=(1:K)
for j=(1:M)
k1 = (j-1)*16+1;
k2 = k1+15;
k = gda_gray2int(genes(i,k1:k2));
mind(i,j) = (mR(j)-mL(j))*(k/65535.999)+mL(j);
end
end
for i=(1:K)
dpre = (mind(i,:)'+[0;mind(i,1:M-1)']+[0;0;mind(i,1:M-2)'])/3; % predicetd data
e = dobs-dpre;
Eind(i) = e'*e;
end
Emax = max(Eind);
% survival of the fitness
F = exp(-c*Eind/Emax).*random('uniform',0,1,K,1);
[F_sorted, k] = sort(F);
k = k(Kold+1:2*Kold,1);
K = Kold;
mind = mind(k,1:M);
Eind = Eind(k,1);
genes = genes(k,1:16*M);
[Emin, k] = min(Eind);
Esave(gen,1) = Emin;
msave(gen,1:M)= mind(k,1:M);
end
% evaluate prediction and plot it
[Emin, i] = min(Eind);
dpre = (mind(i,:)'+[0;mind(i,1:M-1)']+[0;0;mind(i,1:M-2)'])/3; % predicetd data
plot([1:L]',Esave,'k-','LineWidth',2);
fprintf('Done with A %d\n', itrial);
end
% part 2 --------------------------------------------
Nrecom = 5;
imutate = 8;
for itrial=(1:Ntrials)
% fitness coefficient
c=5;
% true model parameters are linear ramp
mt = 0.5+[0:M-1]'/(M-1);
% data average of 3 neighboring model parameters
N=M;
x = [1:M]';
dtrue = (mt+[0;mt(1:M-1)]+[0;0;mt(1:M-2)])/3;
sd=0.025;
dobs = dtrue + random('Normal',0,sd,N,1);
K = 100; % number of individuals in the population
Kold = K;
m0 = ones(M,1); % initial guess for m1
mL = zeros(M,1); % solution > mL
mR = 2*ones(M,1)'; % solution < mR
mind = zeros(K, M); % solution for each individual
Eind = zeros(K,1); % error for each individual
F = zeros(K,1); % fitness of each individual
genes = char(48*ones(K, 16*M)); % genes for each individual
% all individuals initialized to initial guess + random number
for i=(1:M)
r = random('Normal',0,0.025,K,1);
mind(:,i) = m0(i)*ones(K,1)+r;
end
% compute initial error
for i=(1:K)
dpre = (mind(i,:)'+[0;mind(i,1:M-1)']+[0;0;mind(i,1:M-2)'])/3; % predicetd data
e = dobs-dpre;
Eind(i) = e'*e;
end
Emax = max(Eind);
Emaxstart = Emax;
F = exp(-c*Eind/Emax).*random('uniform',0,1,K,1); % initial fitness
for i=(1:K) % convert model parameters to genes
for j=(1:M)
k = floor(65535.999*(mind(i,j)-mL(j))/(mR(j)-mL(j)));
k1 = (j-1)*16+1;
k2 = k1+15;
genes(i,k1:k2) = gda_int2gray(k);
end
end
for gen=(1:L)
% each individual replicates
mind = [mind;mind];
genes = [genes;genes];
Eind = [Eind;Eind];
K = 2*K;
% mutate genes
for i=(1:K)
if( unidrnd(imutate,1) == 1 ) % mutate 1/imutate of individuals
continue;
end
k = unidrnd(16*M,1); % one random mutation
if( genes(i,k)=='0' )
genes(i,k)='1';
else
genes(i,k)='0';
end
end
% recominate genes; Nrecom per generation
for rc = (1:Nrecom)
j1 = unidrnd(K,1); % random individual
j2 = unidrnd(K,1); % another random individual
k = 16*(unidrnd(M-2,1)+1)+1; % always at start of word
g1(1,1:16*M) = genes(j1,:);
g1(1,k:16*M) = genes(j2,k:16*M);
g2(1,1:16*M) = genes(j2,:);
g2(1,k:16*M) = genes(j1,k:16*M);
genes(j1,:) = g1;
genes(j2,:) = g2;
end
% update solution and error from genes
for i=(1:K)
for j=(1:M)
k1 = (j-1)*16+1;
k2 = k1+15;
k = gda_gray2int(genes(i,k1:k2));
mind(i,j) = (mR(j)-mL(j))*(k/65535.999)+mL(j);
end
end
for i=(1:K)
dpre = (mind(i,:)'+[0;mind(i,1:M-1)']+[0;0;mind(i,1:M-2)'])/3; % predicetd data
e = dobs-dpre;
Eind(i) = e'*e;
end
Emax = max(Eind);
% survival of the fitness
F = exp(-c*Eind/Emax).*random('uniform',0,1,K,1);
[F_sorted, k] = sort(F);
k = k(Kold+1:2*Kold,1);
K = Kold;
mind = mind(k,1:M);
Eind = Eind(k,1);
genes = genes(k,1:16*M);
[Emin, k] = min(Eind);
Esave(gen,1) = Emin;
msave(gen,1:M)= mind(k,1:M);
end
fprintf('Done with B %d\n', itrial);
% evaluate prediction and plot it
[Emin, i] = min(Eind);
dpre = (mind(i,:)'+[0;mind(i,1:M-1)']+[0;0;mind(i,1:M-2)'])/3; % predicetd data
plot([1:L]',Esave,'r-','LineWidth',2);
end
fprintf('Caption: Error as a function of generation number for five trials\n');
fprintf('that omit conjugation (black curves) and five that include it (red curves)\n');
% gdama11_18
% Bootstrap estimate of confidence interval
% for solution computed via Newton's method
% This function with these derivatives is being solved
% d(x) = sin( w0 m1 x) + m1 m2;
% dy/dm1 = w0 x cos( w0 m1 x) + m2
% dy/dm2 = m2
% supports Figure 9.18
clear all;
fprintf('gdama11_18\n');
% data are in an auxillary variable, x
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = Dx*[0:N-1]';
% true model parameters
M=2;
mt = [1.21, 1.54]';
% y=f(x, m1, m2);
w0=20;
dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2);
sd=0.4;
dobs = dtrue + random('Normal',0,sd,N,1);
% 2D grid, for plotting purposes only
L = 101;
Dm = 0.02;
m1min=0;
m2min=0;
m1a = m1min+Dm*[0:L-1]';
m2a = m2min+Dm*[0:L-1]';
m1max = m1a(L);
m2max = m2a(L);
% compute error, E, on grid for plotting purposed only
E = zeros(L,L);
for j = [1:L]
for k = [1:L]
dpre = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);
E(j,k) = (dobs-dpre)'*(dobs-dpre);
end
end
% plot error surface
figure();
clf;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [m2min, m2max, m1min, m1max] );
axis ij;
imagesc( [m2min, m2max], [m1min, m1max], E);
colorbar;
xlabel('m2');
ylabel('m1');
plot( mt(2), mt(1), 'go', 'LineWidth', 3 );
Nresamples = 1000;
m1save = zeros(Nresamples,1);
m2save = zeros(Nresamples,1);
for ir = (1:Nresamples)
% resampling with duplications of data
% (first estimate is without resampling)
if( ir==1 )
dresampled = dobs;
xresampled = x;
else
rowindex = unidrnd(N,N,1);
xresampled = x( rowindex );
dresampled = dobs( rowindex );
end
% initial guess and corresponding error
mg = [1.3, 1.5]';
dg = sin(w0*mg(1)*xresampled) + mg(1)*mg(2);
Eg = (dobs-dg)'*(dobs-dg);
% iterate to improve initial guess
Niter=100;
for k = (1:Niter)
dg = sin( w0*mg(1)*xresampled) + mg(1)*mg(2);
dd = dresampled-dg;
Eg=dd'*dd;
Ehis(k+1)=Eg;
G = zeros(N,2);
G(:,1) = w0 * xresampled .* cos( w0 * mg(1) * xresampled ) + mg(2);
G(:,2) = mg(2)*ones(N,1);
% least squares solution
dm = (G'*G)\(G'*dd);
% update
mg = mg+dm;
end
m1save(ir) = mg(1);
m2save(ir) = mg(2);
if( ir==1 )
% plot estimated solution
m1est = mg(1);
m2est = mg(2);
plot( m2est, m1est, 'ro', 'LineWidth', 2 );
else
% plot resampled solutions
plot( mg(2), mg(1), 'w.', 'LineWidth', 2 );
end
end
% mean solution
m1mean = mean(m1save);
m2mean = mean(m2save);
% plot again so is on top
plot( mt(2), mt(1), 'go', 'LineWidth', 3 );
plot( m2est, m1est, 'ro', 'LineWidth', 2 );
fprintf('Caption: Error surface (colors), with bootstrap solutions (white dots),\n');
fprintf('estimated solution (red circle) and true solution (green circle), superimposes\n');
% histogram plot
figure();
clf;
% histogram2 for m1 and m2
Nbins=50;
m1hmin=min(m1save);
m1hmax=max(m1save);
Dm1bins = (m1hmax-m1hmin)/(Nbins-1);
m1bins=m1hmin+Dm1bins*[0:Nbins-1]';
m1hist = hist(m1save,m1bins);
pm1 = m1hist/(Dm1bins*sum(m1hist)); % norm al to pdf
Pm1 = Dm1bins*cumsum(pm1);
m1low=m1bins(find(Pm1>0.025,1));
m1high=m1bins(find(Pm1>0.975,1));
subplot(2,1,1);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
xlabel('m1');
ylabel('p(m1)');
maxpm1=max(pm1);
axis( [m1hmin, m1hmax, 0, maxpm1] );
plot( m1bins, pm1, 'k-', 'LineWidth', 3 );
plot( [m1mean, m1mean]', [0 maxpm1/10]', 'r-', 'LineWidth', 3 );
plot( [m1low, m1low]', [0 maxpm1/20]', 'r-', 'LineWidth', 3 );
plot( [m1high, m1high]', [0 maxpm1/20]', 'r-', 'LineWidth', 3 );
fprintf('estimated m1: %f < %f < %f (95 percent confidence)\n', m1low, m1est, m1high );
m2hmin=min(m2save);
m2hmax=max(m2save);
Dm2bins = (m2hmax-m2hmin)/(Nbins-1);
m2bins=m2hmin+Dm2bins*[0:Nbins-1]';
m2hist = hist(m2save,m2bins);
pm2 = m2hist/(Dm2bins*sum(m2hist)); % norm al to pdf
Pm2 = Dm2bins*cumsum(pm2);
m2low=m2bins(find(Pm2>0.025,1));
m2high=m2bins(find(Pm2>0.975,1));
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
xlabel('m2');
ylabel('p(m2)');
maxpm2=max(pm2);
axis( [m2hmin, m2hmax, 0, maxpm2] );
plot( m2bins, pm2, 'k-', 'LineWidth', 3 );
plot( [m2est, m2est]', [0 maxpm2/10]', 'r-', 'LineWidth', 3 );
plot( [m2low, m2low]', [0 maxpm2/20]', 'r-', 'LineWidth', 3 );
plot( [m2high, m2high]', [0 maxpm2/20]', 'r-', 'LineWidth', 3 );
fprintf('Caption: (top) Empirical pdfp(m1) for model parameter m1. (black curv) with mean (large red bar) and\n');
fprintf('95 percent confidence intervals (small red bars).\n')
fprintf('\n');
fprintf('estimated m2: %f < %f < %f (95 percent confidence)\n', m2low, m2est, m2high );