Untitled

% gdama04
% 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-04, Edit and integration with text

% gdama04_01

% least squares fit to synthetic data

clear all;

fprintf('gdama04_01\n');

gdama04_01

% z's

N=30;

zmin=0;

zmax=10;

z = sort(random('Uniform',zmin,zmax,N,1));

% d = a + b*z + noise

a=2.0;

b=1.0;

sigmad=1.0;

sigmad2 = sigmad^2;

dobs = a+b*z+random('Normal',0,sigmad,N,1);

% Data kernel

M=2;

G=[ones(N,1), z];

% least squares fit

GTG = G'*G;

mest = GTG\(G'*dobs);

covm = sigmad2 * inv(GTG);

% print solution

fprintf('Estimated solution:\n');

Estimated solution:

fprintf('m1: %.2f +/- %.2f (95%%)', mest(1), 2*sqrt(covm(1,1)));

m1: 1.96 +/- 0.79 (95%)

fprintf('m2: %.2f +/- %.2f (95%%)', mest(2), 2*sqrt(covm(2,2)));

m2: 0.98 +/- 0.14 (95%)

fprintf(' \n');

% predicted data

dpre = G*mest;

e = dobs - dpre;

[emax, iemax] = max(abs(e));

% plot

figure(1);

clf;

% plot scale

pdmin=0;

pdmax=15;

% plot of data & straight line fit

subplot(1,2,1);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, pdmin, pdmax ]' );

plot( z, dobs, 'ro', 'LineWidth', 2);

plot( z, dpre, 'b-', 'LineWidth', 2);

xlabel('z');

ylabel('d');

% plot illustrating definition of prediction error

subplot(1,2,2);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, pdmin, pdmax ]' );

i=iemax;

plot( z(i), dobs(i), 'ro', 'LineWidth', 2);

plot( z, dpre, 'b-', 'LineWidth', 2);

plot( [z(i), z(i)]', [pdmin, dpre(i)]', 'k:', 'LineWidth', 2);

plot( [zmin, z(i)]', [dpre(i), dpre(i)]', 'k:', 'LineWidth', 2);

plot( [zmin, z(i)]', [dobs(i), dobs(i)]', 'k:', 'LineWidth', 2);

xlabel('z');

ylabel('d');

fprintf('Caption: (left) Observed data (circles) and predicted data (blue).\n');

Caption: (left) Observed data (circles) and predicted data (blue).

fprintf('(right) depiction of one predicted datum correponding to a single\n');

(right) depiction of one predicted datum correponding to a single

fprintf('observed datum (red).\n');

observed datum (red).

% gdama04_02

% comparison of norms

clear all;

fprintf('gdama04_02\n');

gdama04_02

% z's

N=31;

zmin=0;

zmax=10;

Dz=(zmax-zmin)/(N-1);

z=zmin+Dz*[0:N-1]';

% random errors

e = random('Uniform',-1,1,N,1);

% variery of powers of errors

e1 = abs(e);

E1 = sum(e1);

e2 = abs(e).^2;

E2 = sqrt(sum(e1));

e10 = abs(e).^10;

E10 = sum(e10)^0.1;

fprintf('Errors E1 %f E2 %f E10 %f\n', E1, E2, E10);

Errors E1 16.971244 E2 4.119617 E10 1.097731

% plot of errors

figure(1);

clf;

% plot scale

pemin=-1;

pemax=1;

% plot error

subplot(4,1,1);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, pemin, pemax ]' );

% improvise bar chart

for i=(1:N-1)

plot( [z(i), z(i), z(i+1), z(i+1) ]', [0, e(i), e(i), 0]', 'k-', 'LineWidth', 2 );

end

xlabel('z');

ylabel('e');

plot( [zmin, zmax ]', [0, 0]', 'k-', 'LineWidth', 2 );

% plot |error|

subplot(4,1,2);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, pemin, pemax ]' );

% improvise bar chart

for i=(1:N-1)

plot( [z(i), z(i), z(i+1), z(i+1) ]', [0, e1(i), e1(i), 0]', 'k-', 'LineWidth', 2 );

end

xlabel('z');

ylabel('|e|');

plot( [zmin, zmax ]', [0, 0]', 'k-', 'LineWidth', 2 );

% plot |error|^2

subplot(4,1,3);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, pemin, pemax ]' );

% improvise bar chart

for i=(1:N-1)

plot( [z(i), z(i), z(i+1), z(i+1) ]', [0, e2(i), e2(i), 0]', 'k-', 'LineWidth', 2 );

end

xlabel('z');

ylabel('|e|^2');

plot( [zmin, zmax ]', [0, 0]', 'k-', 'LineWidth', 2 );

% plot |error|^10

subplot(4,1,4);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, pemin, pemax ]' );

% improvise bar chart

for i=(1:N-1)

plot( [z(i), z(i), z(i+1), z(i+1) ]', [0, e10(i), e10(i), 0]', 'k-', 'LineWidth', 2 );

end

xlabel('z');

ylabel('|e|^1^0');

plot( [zmin, zmax ]', [0, 0]', 'k-', 'LineWidth', 2 );

fprintf('Caption: Comparisons of norms (top plot) Observed error, 1\n');

Caption: Comparisons of norms (top plot) Observed error, 1

fprintf('(second plot) |e|, (third plot) |e|^2, (fourth plot) |e|^10.');

(second plot) |e|, (third plot) |e|^2, (fourth plot) |e|^10.

% gdama04_03

% Ln fits of straight line to to synthetic data, via grid search.

clear all;

% z-axis

N=30;

zmin=0;

zmax=10;

z = sort(random('Uniform',zmin,zmax,N,1));

% linear model

% d = a + b*z + noise

a=2.0;

b=1.0;

sd=0.5;

dobs = a+b*z+random('Normal',0,sd,N,1);

% one terrible outlier

dobs(N)=random('uniform',0,dobs(N),1,1);

% grid (a=intercept, b=slope)

Na=1001;

amin=-20;

amax=20;

Da = (amax-amin)/(Na-1);

Nb=1001;

bmin=-6;

bmax=6;

Db = (bmax-bmin)/(Nb-1);

% populate grid with errors

E1 = zeros(Na,Nb);

E2 = zeros(Na,Nb);

Einf = zeros(Na,Nb);

for i=(1:Na)

for j=(1:Nb)

ao = amin+Da*(i-1);

bo = bmin+Db*(j-1);

dpre = ao + bo*z;

e = dobs-dpre;

abse = abs(e);

E1(i,j)=sum(abse);

E2(i,j)=sum(abse.^2);

Einf(i,j)=sum(abse.^20); % cheating; using large but finite power

end

% Important! Examine the error surface to make sure that

% the minimum is within the grid!

gda_draw(log10(E1),'caption L1',' ',log10(E2), 'caption L2',' ',log10(Einf),'caption Linf');

fprintf('Caption: Error surfaces (left) L1 norm, (middle) L2 norm, (right0 Linf norm.\n');

Caption: Error surfaces (left) L1 norm, (middle) L2 norm, (right0 Linf norm.

% find minimum error

[Ep, p] = min(E1);

[E1min, c1] = min(Ep);

r1 = p(c1);

a1 = amin+Da*(r1-1);

b1 = bmin+Db*(c1-1);

dpre1 = a1 + b1*z;

[Ep, p] = min(E2);

[E2min, c2] = min(Ep);

r2 = p(c2);

a2 = amin+Da*(r2-1);

b2 = bmin+Db*(c2-1);

dpre2 = a2 + b2*z;

[Ep, p] = min(Einf);

[Einfmin, cinf] = min(Ep);

rinf = p(cinf);

ainf = amin+Da*(rinf-1);

binf = bmin+Db*(cinf-1);

dpreinf = ainf + binf*z;

% write our intecept & slope

fprintf('Solution (intercept a and slope b)\n');

Solution (intercept a and slope b)

fprintf('True: a %.3f b %.3f\n', a, b );

True: a 2.000 b 1.000

fprintf('L1: a %.3f b %.3f\n', a1, b1 );

L1:   a 2.000 b 1.008

fprintf('L2: a %.3f b %.3f\n', a2, b2 );

L2:   a 2.400 b 0.888

fprintf('Linf: a %.3f b %.3f\n', ainf, binf );

Linf: a 4.200 b 0.456

fprintf('\n');

% plot results

figure(2);

clf;

% plot bounds

pdmin=0;

pdmax=15;

% plot data and fits

set(gca,'LineWidth',2);

hold on;

axis( [zmin, zmax, pdmin, pdmax ]' );

plot( z, dobs, 'ko', 'LineWidth', 2);

plot( z, dpre1, 'r-', 'LineWidth', 2);

plot( z, dpre2, 'g-', 'LineWidth', 2);

plot( z, dpreinf, 'b-', 'LineWidth', 2);

xlabel('z');

ylabel('d');

fprintf('Caption: Straight line fits to observed data (circles), showing\n');

Caption: Straight line fits to observed data (circles), showing

fprintf('L1 fit (red), L2 fit (green), Linf fit (blue).\n');

L1 fit (red), L2 fit (green), Linf fit (blue).

% gdama04_04

% comparison of long-tailed and short tailed pdf's

clear all;

fprintf('gdama04_04\n');

gdama04_04

% d-axis

Dd = 0.1;

N = 101;

d = Dd*[0:N-1]';

dmin=0;

dmax=10;

% short tailed pdf: Normal pdf

dbar=5;

d2=(d-dbar).^2;

sd = 1.0;

dbar = 5.0;

p1 = exp(-0.5*d2/(sd^2))/(sqrt(2*pi)*sd);

A1 = Dd*sum(p1);

% long-tailed distribution: Cauchy–Lorentz distribution

g = 1;

p2 = 1./(pi.*g.*(1+(d2./(g^2))));

A2 = Dd*sum(p2);

fprintf('check on areas: true %f est1 %f est2 %f\n', 1, A1, A2);

check on areas: true 1.000000 est1 1.000000 est2 0.875551

figure(1);

clf;

% plot long-tailed pdf

subplot(1,2,2);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [dmin, dmax, 0, 0.5 ] );

plot(d,p1,'r-','LineWidth',3);

xlabel('d');

ylabel('p(d)');

% plot short-tailed pdf

subplot(1,2,1);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [dmin, dmax, 0, 0.5 ] );

plot(d,p2,'b-','LineWidth',3);

xlabel('d');

ylabel('p(d)');

fprintf('Caption: (Left) long-tailed pdf, (right) short-tailed pdf\n');

Caption: (Left) long-tailed pdf, (right) short-tailed pdf

% gdama04_05

% Kepler's 3rd law example, period^2 = radius^3

clear all;

fprintf('gdama04_05\n');

gdama04_05

% read data

D=load('../data/planetary.txt');

radius=D(:,1)/1e9; % work in units of 10^9 km

period=D(:,2)/1000; % and units of 10^3 days

N=length(radius);

% take radius^2 to be the observation

% and period to be the auxillary variable

% (which is ok, since time can be measured so accurately)

dobs = radius.^3;

z=period;

% data kernel

M=3;

G=[ones(N,1), z, z.^2];

% least squares solution

GTG = G'*G;

mest = GTG\(G'*dobs);

% predicted solutiona nd its error

dpre = G*mest;

e = dobs-dpre;

sigmad2 = e'*e / (N-M); % posterior variance of data

covm = sigmad2*inv(GTG);

fprintf('Estimated solution:\n');

Estimated solution:

fprintf('m1 %f +1 %f (95%)\n', mest(1), 2*sqrt(covm(1,1)));

m1 -0.025331 +1 0.091256 (95

fprintf('m2 %f +1 %f (95%)\n', mest(2), 2*sqrt(covm(2,2)));

m2 0.009980 +1 0.008962 (95

fprintf('m3 %f +1 %f (95%)\n', mest(3), 2*sqrt(covm(3,3)));

m3 0.024996 +1 0.000105 (95

fprintf('\n');

% plot smooth parabola, so use lots of z's

zeval=[0:250]';

deval=mest(1)+mest(2)*zeval+mest(3)*zeval.^2;

figure(1);

clf;

subplot(2,1,1);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [0, 100, 0, 250 ]');

plot(z,dobs,'ro','LineWidth',3);

plot(zeval,deval,'k-','LineWidth',3);

xlabel('period, Kdays');

ylabel('radius^3, Tm^3');

subplot(2,1,2);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [0, 100, -0.5, 0.5 ]');

plot(z,e,'ro','LineWidth',3);

plot( [0, 100], [0, 0], 'k-','LineWidth',2);

xlabel('period, Kdays');

ylabel('error, Tm^3');

fprintf('Caption: Kepler''s third law example. (Top) Observed data (red)\n');

Caption: Kepler's third law example. (Top) Observed data (red)

fprintf('least squares fit (black). (Bottom) Prediction error (red).')

least squares fit (black). (Bottom) Prediction error (red).

% gdama04_06

% planar fit to depths of earthquakes in the Kurile subduction zone

clear all;

% read data

D=load('../data/kurile_eqs.txt');

lat=D(:,1);

lon=D(:,2);

depth=D(:,3);

% convert to km

x = 111.12*cos((pi/180)*mean(lat))*(lon-min(lon));

y = 111.12*(lat-min(lat));

dobs=-depth;

N=length(-depth);

% three-d plot of data

% use rotate button to examine it prom different perspectives!

figure(1);

clf;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

pxmin=0; pxmax=1200;

pymin=0; pymax=1200;

pzmin=-700; pzmax=0;

axis( [pxmin, pxmax, pymin, pymax, pzmin, pzmax]' );

plot3(x,y,dobs,'ko','LineWidth',2);

% improvise 3D box

plot3( [pxmin,pxmin], [pymin,pymin], [pzmin,pzmax], 'k-', 'LineWidth', 2 );

plot3( [pxmin,pxmin], [pymin,pymax], [pzmin,pzmin], 'k-', 'LineWidth', 2 );

plot3( [pxmin,pxmax], [pymin,pymin], [pzmin,pzmin], 'k-', 'LineWidth', 2 );

plot3( [pxmax,pxmax], [pymax,pymax], [pzmax,pzmin], 'k-', 'LineWidth', 2 );

plot3( [pxmax,pxmax], [pymax,pymin], [pzmax,pzmax], 'k-', 'LineWidth', 2 );

plot3( [pxmax,pxmin], [pymax,pymax], [pzmax,pzmax], 'k-', 'LineWidth', 2 );

plot3( [pxmax,pxmin], [pymin,pymin], [pzmax,pzmax], 'k-', 'LineWidth', 2 );

plot3( [pxmax,pxmax], [pymin,pymin], [pzmax,pzmin], 'k-', 'LineWidth', 2 );

plot3( [pxmin,pxmin], [pymax,pymin], [pzmax,pzmax], 'k-', 'LineWidth', 2 );

plot3( [pxmin,pxmin], [pymax,pymax], [pzmax,pzmin], 'k-', 'LineWidth', 2 );

plot3( [pxmax,pxmax], [pymax,pymin], [pzmin,pzmin], 'k-', 'LineWidth', 2 );

plot3( [pxmax,pxmin], [pymax,pymax], [pzmin,pzmin], 'k-', 'LineWidth', 2 );

% solve inverse problem of fitting a plane to data

M=3;

G = [ones(N,1), x, y];

GTG = G'*G; % Gram matrix

mest = GTG\(G'*dobs); % least squres solution

dpre=G*mest; % predicted data

e = dobs-dpre; % prediction error

sigmad2 = e'*e/(N-M); % posterior variance of data

covm = sigmad2*inv(GTG); % covariance of solution

fprintf('Estimated solution:\n');

Estimated solution:

fprintf('m1 %f +1 %f (95%)\n', mest(1), 2*sqrt(covm(1,1)));

m1 -94.887492 +1 5.932558 (95

fprintf('m2 %f +1 %f (95%)\n', mest(2), 2*sqrt(covm(2,2)));

m2 0.659768 +1 0.022970 (95

fprintf('m3 %f +1 %f (95%)\n', mest(3), 2*sqrt(covm(3,3)));

m3 -0.657462 +1 0.020261 (95

fprintf('\n');

% display fit as mesh in 3D space

M=31;

xx=pxmin+(pxmax-pxmin)*[0:M-1]'/(M-1);

yy=pymin+(pymax-pymin)*[0:M-1]'/(M-1);

[X,Y]=meshgrid( xx, yy);

Z = mest(1)+mest(2)*X+mest(3)*Y;

[Ni,Nj]=size(Z);

for i=(1:Ni)

for j=(1:Nj)

if ( (Z(i,j)>0) || (Z(i,j)<pzmin) )

Z(i,j)=nan(1);

end

mesh(X,Y,Z);

view(3);

fprintf('Caption: Subdiction zone earthquake locations (black), best-fitting\n');

Caption: Subdiction zone earthquake locations (black), best-fitting

fprintf('planar surface (colored mesh).\n');

planar surface (colored mesh).

% gdama04_07

% Least squares estimation of layer interface

% properties from reflection coefficient data

% note that this version used functions to simplicy

% the structure of the code. In MATLAB, a function

% "myfcn" that takes arguments A and B are returns

% values C=2A and D=2B is codes as follows

% function [C,D] = myfcn(A,B)

% C=2*A;

% D=2*B;

% end

clear all;

fprintf('gdama04_07\n');

gdama04_07

% q=theta is the angle of incidence

Nq = 21;

qmin = 0.0;

qmax = 25.0;

Dq = (qmax-qmin)/(Nq-1);

q = qmin + Dq*[0:Nq-1]';

% medium properties

a = 5000; % compressional velocity

b = 3000; % shear velocity

r = 2500; % density

% interface properties

Dror = 0.07; % normalized jump in compressional velocity

Daoa = 0.10; % normalized jump in shear velocity

Dbob = 0.11; % normalized jump in density velocity

mtrue = [Dror, Daoa, Dbob]'; % group into model parameters

% loop over angle of incidence

N=4*Nq; % number of data, 4 reflection coefficients times Nq angles

M=3; % number of model parameters

G = zeros( N, M ); % initialize data kernel

for iq = [1:Nq]

% Reflection coefficients for an incident P wave

% I have coded them in a function, because they are messy

[ PDPUr, PDPUa, PDPUb, PDSUr, PDSUa, PDSUb ] = gda_PD( q(iq), r, a, b );

PDPU = PDPUr*Dror + PDPUa*Daoa + PDPUb*Dbob;

PDSU = PDSUr*Dror + PDSUa*Daoa + PDSUb*Dbob;

% Reflection coefficients for an incident S wave

% I have coded them in a function, because they are messy

[ SDPUr, SDPUa, SDPUb, SDSUr, SDSUa, SDSUb ] = gda_SD( q(iq), r, a, b );

SDPU = SDPUr*Dror + SDPUa*Daoa + SDPUb*Dbob;

SDSU = SDSUr*Dror + SDSUa*Daoa + SDSUb*Dbob;

% data kernel

G(iq,:) = [PDPUr, PDPUa, PDPUb];

G(iq+Nq,:) = [PDSUr, PDSUa, PDSUb];

G(iq+2*Nq,:) = [SDPUr, SDPUa, SDPUb];

G(iq+3*Nq,:) = [SDSUr, SDSUa, SDSUb];

end

% true reflection coefficients

dtrue = G*mtrue;

% plot the reflection coefficents as a function of angle

figure(1);

clf;

set(gca,'LineWidth',3);

hold on;

axis( [qmin, qmax, -0.1, 0.1] );

plot( [qmin, qmax]', [0,0]', 'k:', 'LineWidth', 2 );

plot( q, dtrue(1:Nq), 'k-', 'LineWidth', 3 );

plot( q, dtrue(1+Nq:2*Nq), 'r-', 'LineWidth', 3 );

plot( q, dtrue(1+2*Nq:3*Nq), 'g-', 'LineWidth', 3 );

plot( q, dtrue(1+3*Nq:4*Nq), 'b-', 'LineWidth', 3 );

% make synthetic noisy data and plot it

sigmad = 0.005;

dobs = dtrue + random('Normal',0,sigmad,N,1);

plot( q, dobs(1:Nq), 'k-', 'LineWidth', 2 );

plot( q, dobs(1+Nq:2*Nq), 'r-', 'LineWidth', 2 );

plot( q, dobs(1+2*Nq:3*Nq), 'g-', 'LineWidth', 2 );

plot( q, dobs(1+3*Nq:4*Nq), 'b-', 'LineWidth', 2 );

% least squares estimate of the layer interface properties

GTGinv = inv(G'*G);

mest = GTGinv*(G'*dobs);

dpre = G*mest;

e = dobs-dpre;

E = e'*e;

sigmadest2 = E/(N-M);

sigmadest = sqrt(sigmadest2);

fprintf('sigmad true %.4f est %.4f\n', sigmad, sigmadest);

sigmad true 0.0050 est 0.0051

fprintf('\n');

covm = (sigmadest2)*GTGinv; % covariance of the estimates

fprintf('cov m\n');

cov m

fprintf('%.5f %.5f %.5f\n', covm(1,1), covm(1,2), covm(1,3) );

0.00007 -0.00007 -0.00008

fprintf('%.5f %.5f %.5f\n', covm(2,1), covm(2,2), covm(2,3) );

-0.00007 0.00008 0.00008

fprintf('%.5f %.5f %.5f\n', covm(3,1), covm(3,2), covm(3,3) );

-0.00008 0.00008 0.00009

fprintf('\n');

sigmam1 = sqrt(covm(1,1));

sigmam2 = sqrt(covm(2,2));

sigmam3 = sqrt(covm(3,3));

fprintf('m1: true %.2f est %.2f +/- %.2f (95%%)\n', mtrue(1), mest(1), 2*sigmam1 );

m1: true 0.07 est 0.06 +/- 0.02 (95%)

fprintf('m2: true %.2f est %.2f +/- %.2f (95%%)\n', mtrue(2), mest(2), 2*sigmam2 );

m2: true 0.10 est 0.11 +/- 0.02 (95%)

fprintf('m3: true %.2f est %.2f +/- %.2f (95%%)\n', mtrue(3), mest(3), 2*sigmam3 );

m3: true 0.11 est 0.12 +/- 0.02 (95%)

fprintf('\n');

plot( q, dpre(1:Nq), 'ko', 'LineWidth', 2 );

plot( q, dpre(1+Nq:2*Nq), 'ro', 'LineWidth', 2 );

plot( q, dpre(1+2*Nq:3*Nq), 'go', 'LineWidth', 2 );

plot( q, dpre(1+3*Nq:4*Nq), 'bo', 'LineWidth', 2 );

xlabel('theta (deg)');

ylabel('amplitude')

fprintf('Caption: Reflection coefficients PDPU (black), PDSU (red), SDPU (green) and SDSU (blue),\n');

Caption: Reflection coefficients PDPU (black), PDSU (red), SDPU (green) and SDSU (blue),

fprintf(' with true data (bold curves), observed data (dots) and predicted data (thin curves).\n');

 with true data (bold curves), observed data (dots) and predicted data (thin curves).

figure();

clf;

subplot(1,3,1);

set(gca,'LineWidth',3);

hold on;

axis( [0, 1, 0, 0.25] );

plot( [0, 1]', [mtrue(1),mtrue(1)]', 'k-', 'LineWidth', 5 );

plot( [0, 1]', [mest(1),mest(1)]', 'r-', 'LineWidth', 3 );

plot( [0.5, 0.5]', [mest(1)-2*sigmam1,mest(1)+2*sigmam1]', 'r-', 'LineWidth', 3 );

ylabel('density perturbation');

subplot(1,3,2);

set(gca,'LineWidth',3);

hold on;

axis( [0, 1, 0, 0.25] );

plot( [0, 1]', [mtrue(2),mtrue(2)]', 'k-', 'LineWidth', 5 );

plot( [0, 1]', [mest(2),mest(2)]', 'r-', 'LineWidth', 3 );

plot( [0.5, 0.5]', [mest(2)-2*sigmam2,mest(2)+2*sigmam2]', 'r-', 'LineWidth', 3 );

ylabel('compressional velocity perturbation');

subplot(1,3,3);

set(gca,'LineWidth',3);

hold on;

axis( [0, 1, 0, 0.25] );

plot( [0, 1]', [mtrue(3),mtrue(3)]', 'k-', 'LineWidth', 5 );

plot( [0, 1]', [mest(3),mest(3)]', 'r-', 'LineWidth', 3 );

plot( [0.5, 0.5]', [mest(3)-2*sigmam3,mest(3)+2*sigmam3]', 'r-', 'LineWidth', 3 );

ylabel('shear velocity perturbation');

fprintf('Caption: (Left) Density perturbation, (middle) compressional velocity perturbation\n');

Caption: (Left) Density perturbation, (middle) compressional velocity perturbation

fprintf('and (right) shear velocity perturbation, showing true values (black), predicted\n');

and (right) shear velocity perturbation, showing true values (black), predicted

fprintf('values with 95 percent confidence interval (red).')

values with 95 percent confidence interval (red).

% gdama04_08
% sparse matrix and bicg() solver for 5x5 anti-identity matrix
 
clear all;
fprintf('gdama04_08\n');
gdama04_08
 
% this uses the Generalized Least Squares (GLS) formalism
% in which the NxM data kernel, G, is the top part of the
% LxM GLS matrix, F.  However, in this case, there is no
% bottom part, so L=N
 
% global variable for sparse version of F
clear gdaFsparse;
global gdaFsparse;
 
N=5;
L=N;
M=5;
 
% build element arrays
P=0;
Pmax = 100;
Fr = zeros(Pmax,1);
Fc = zeros(Pmax,1);
Fv = zeros(Pmax,1);
 
% build array row-wide
for i = (1:L)    % row i
    j=N-i+1;     % column j
    P=P+1;       % increment element number
    Fr(P,1)=i;   % row index
    Fc(P,1)=j;   % column index
    Fv(P,1)=1.0; % value
end
 
% true solution and data
mtrue = [1, 2, 3, 4, 5]';
dtrue = [5, 4, 3, 2 ,1]';
 
% observed data is true data
dobs=dtrue;
 
% print element arrays
fprintf('element arrays of size P=%d\n',P);
element arrays of size P=5
fprintf('i rF    CF     vF\n');
i rF    CF     vF
for i = (1:N)    % row i
    fprintf('%d %.0f %.0f %.3f\n', i, Fr(i,1),  Fc(i,1), Fv(i,1) );
end
1 1 5 1.000
2 2 4 1.000
3 3 3 1.000
4 4 2 1.000
5 5 1 1.000
fprintf('\n');
 
% build sparse array
L=N;
 
Fr=Fr(1:P); Fc=Fc(1:P); Fv=Fv(1:P); % truncate
gdaFsparse = sparse(Fr,Fc,Fv,L,M); % build sparse matrix
clear Fr Fc Fv;  % delete no longer needed arrays
 
% print F (huge gyrations to print a nice rectangular block!)
fprintf('F:\n');
F:
F=full(gdaFsparse);
for i=(1:L)
    line = [];
    for j=(1:M)
        str = sprintf('%2.0f ', F(i,j) );
        line = [line, str];
    end
    fprintf('%s\n',line);
end
 0  0  0  0  1 
 0  0  0  1  0 
 0  0  1  0  0 
 0  1  0  0  0 
 1  0  0  0  0 
fprintf('\n');
 
% least squared solution
tol = 1e-12;     % tolerance
maxit = 3*(M+N); % maximum number of iterations allowed
FTf = gda_FTFrhs(dobs);
mest=bicg(@gda_FTFmul,FTf,tol ,maxit);
bicg converged at iteration 1 to a solution with relative residual 0.
 
fprintf('i    mtrue  mest\n');
i    mtrue  mest
for i=(1:M)
    fprintf('%d %.2f %.2f\n', i, mtrue(i), mest(i) );
end
1 1.00 1.00
2 2.00 2.00
3 3.00 3.00
4 4.00 4.00
5 5.00 5.00
fprintf('\n');

% gdama04_09
% sparse matrix and bicg() solver for case where number of
% non-zero elements on a row of F is unknown.
 
clear all;
fprintf('gdama04_09\n');
gdama04_09
 
% this uses the Generalized Least Squares (GLS) formalism
% in which the NxM data kernel, G, is the top part of the
% LxM GLS matrix, F.  However, in this case, there is no
% bottom part, so L=N
 
% global variable for sparse version of F
clear gdaFsparse;
global gdaFsparse;
 
N=5;
L=N;
M=5;
 
% build element arrays
P=0;
Pmax = 100;
Fr = zeros(Pmax,1);
Fc = zeros(Pmax,1);
Fv = zeros(Pmax,1);
 
for i = (1:L)    % row i
    Frow = zeros(M,1);  % i-th row of F
    
    % now set elements of Frow here
    % I cooked up a case where F is an
    % identity matrix except that 1 or 2
    % randomly-chosen elements of each row
    % have been overwritten by a random value
    Frow(i,1)=1.0;
    k = randi(2);
    for j=(1:k)
        kk = randi(M);  % randonm integer 
        v = randn();    % random value
        Frow(kk)=v;
    end
 
    % now extract non-zero values from Frow
    for j = (1:M)
        if( Frow(j) ~= 0.0 )
            P=P+1; % increment element number
            Fr(P,1)=i;         % row index
            Fc(P,1)=j;         % column index
            Fv(P,1)=Frow(j,1); % value
        end
    end
end
 
% print element arrays
fprintf('element arrays of size P=%d\n',P);
element arrays of size P=11
fprintf('i rF    CF     vF\n');
i rF    CF     vF
for i = (1:N)    % row i
    fprintf('%d %.0f %.0f %.3f\n', i, Fr(i,1),  Fc(i,1), Fv(i,1) );
end
1 1 1 1.000
2 1 3 -2.486
3 1 4 -2.192
4 2 2 1.000
5 2 3 -0.948
fprintf('\n');
 
 
L=N;
% delete unused parts of arrays
Fr=Fr(1:P);
Fc=Fc(1:P);
Fv=Fv(1:P);
gdaFsparse = sparse(Fr,Fc,Fv,L,M);  % build sparse array
clear Fr Fc Fv;  % delete no longer needed arrays
 
% print F (huge gyrations to print a nice rectangular block!)
fprintf('F:\n');
F:
F=full(gdaFsparse);
for i=(1:L)
    line = [];
    for j=(1:M)
        str = sprintf('%2.0f ', F(i,j) );
        line = [line, str];
    end
    fprintf('%s\n',line);
end
 1  0 -2 -2  0 
 0  1 -1  0  0 
 0  0  1  0  1 
 0  0  0  1 -0 
 0  0  0  0 -2 
fprintf('\n');
 
% true solution and data
mtrue = [1, 2, 3, 4, 5]';
dtrue = gdaFsparse*mtrue;
 
% observed data is true data
dobs=dtrue;
 
% least squared solution
tol = 1e-12;     % tolerance
maxit = 3*(M+N); % maximum number of iterations allowed
FTf = gda_FTFrhs(dobs);
mest=bicg(@gda_FTFmul,FTf,tol ,maxit);
bicg converged at iteration 5 to a solution with relative residual 9.7e-13.
 
fprintf('i    mtrue  mest\n');
i    mtrue  mest
for i=(1:M)
    fprintf('%d %.2f %.2f\n', i, mtrue(i), mest(i) );
end
1 1.00 1.00
2 2.00 2.00
3 3.00 3.00
4 4.00 4.00
5 5.00 5.00
fprintf('\n');

% gdama04_10

% straight line problem

% least squares fit to synthetic data

% using two solution methods

% standard solver

% bicg() solver

clear all;

fprintf('gdama04_10\n');

gdama04_10

% this uses the Generalized Least Squares (GLS) formalism

% in which the NxM data kernel, G, is the top part of the

% LxM GLS matrix, F. However, in this case, there is no

% bottom part, so L=N

% global variable for sparse version of F

clear gdaFsparse;

global gdaFsparse;

% z's

N=30;

zmin=0;

zmax=10;

z = sort(random('Uniform',zmin,zmax,N,1));

% d = a + b*z + noise

a=2.0;

b=1.0;

mtrue=[a,b]';

sd=1.0;

dobs = a+b*z+random('Normal',0,sd,N,1);

% data kernel

M=2;

G=[ones(N,1), z];

% standard matrix solution

mest1 = (G'*G)\(G'*dobs);

dpre1 = G*mest1;

% build element arrays

P=0;

Pmax = 100;

Fr = zeros(Pmax,1);

Fc = zeros(Pmax,1);

Fv = zeros(Pmax,1);

for i = (1:N) % row i

j=1;

P=P+1; % increment element number

Fr(P,1)=i; % row index

Fc(P,1)=j; % column index

Fv(P,1)=1.0; % value

j=2;

P=P+1; % increment element number

Fr(P,1)=i; % row index

Fc(P,1)=j; % column index

Fv(P,1)=z(i); % value

end

% print element arrays

fprintf('element arrays of size P=%d\n',P);

element arrays of size P=60

fprintf('i rF CF vF\n');

i rF    CF     vF

for i = (1:N) % row i

fprintf('%d %.0f %.0f %.3f\n', i, Fr(i,1), Fc(i,1), Fv(i,1) );

end

fprintf('\n');

% build sparse array

L=N;

% delete unused parts of arrays

Fr=Fr(1:P);

Fc=Fc(1:P);

Fv=Fv(1:P);

gdaFsparse = sparse(Fr,Fc,Fv,L,M);

clear Fr Fc Fv; % delete no longer needed arrays

% print F (huge gyrations to print a nice rectangular block!)

fprintf('F:\n');

F:

F=full(gdaFsparse);

for i=(1:L)

line = [];

for j=(1:M)

str = sprintf('%2.0f ', F(i,j) );

line = [line, str];

end

fprintf('%s\n',line);

end

fprintf('\n');

% least squared solution

tol = 1e-12; % tolerance

maxit = 3*(M+N); % maximum number of iterations allowed

FTf=gda_FTFrhs(dobs);

mest2=bicg(@gda_FTFmul,FTf,tol,maxit);

bicg converged at iteration 2 to a solution with relative residual 1.2e-14.

fprintf('i mtrue mest1 mest2\n');

i    mtrue  mest1  mest2

for i=(1:M)

fprintf('%d %.2f %.2f %.2f\n', i, mtrue(i), mest1(i), mest2(i) );

end

1 2.00 2.35 2.35
2 1.00 0.92 0.92

fprintf('\n');

% predicted data

dpre2 = gdaFsparse*mest2;

% plot

figure();

clf;

hold on;

pdmin=0;

pdmax=15;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

axis( [zmin, zmax, pdmin, pdmax ]' );

plot( z, dobs, 'ko', 'LineWidth', 2);

plot( z, dpre1, 'r-', 'LineWidth', 3);

plot( z, dpre2, 'g--', 'LineWidth', 2);

xlabel('z');

ylabel('d');

fprintf('Caption: Data, with observed (circles), predicted (non-sparse method)\n');

Caption: Data, with observed (circles), predicted (non-sparse method)

fprintf('(red), prediced (sparse methid) (green)\n');

(red), prediced (sparse methid) (green)

% gdama04_11

% example of data smoothing and gap filling with

% weighted damped least squares

clear all;

fprintf('gdama04_11\n');

gdama04_11

% global variable for sparse version of F

clear gdaFsparse;

global gdaFsparse;

% the data is a sinusoid in an auxillary varible z;

M=101;

Dz = 1.0;

z = Dz*[0:M-1]'; % z is an auxiallty variable

zmax = max(z);

mtrue = sin( 3*pi*z/zmax );

% Evaluate sine wave at the following randomly-chosen z indices

index = [1, unidrnd(M,1,floor(M/2)),M]';

N = length(index);

zobs=z(index);

% observed data is sinusoid plus random noise

sigmad=0.05;

dobs = sin( 3*pi*zobs/zmax ) + random('Normal',0,sigmad,N,1);

% PART 1: First derivative smoothing

% build element arrays

L=N+M;

P=0;

Pmax = L*2*(M-2)+100;

Fr = zeros(Pmax,1);

Fc = zeros(Pmax,1);

Fv = zeros(Pmax,1);

f = zeros(L,1);

% the N data equations are just m = dobs. The only trick is lining

% up the corresponding elements of m and dobs, since they are not of

% the same length

nr=0; % row counter

for i = (1:N)

j=index(i);

nr=nr+1; % increment row counter

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = j; % column index

Fv(P,1) = 1.0; % value

f(nr,1) = dobs(i,1); % r.h.s. is data

end

% now implement 1st derivative flatness constraints but the last m

% (that makes M-1 rows)

epsilon=1.0;

rDz = 1/Dz;

rDz2 = 1/Dz^2;

for i = (1:M-1)

j=i;

nr=nr+1; % increment row counter

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = j; % column index

Fv(P,1) = -rDz; % value

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = j+1; % column index

Fv(P,1) = rDz; % value

f(nr,1) = 0.0; % r.h.s. is data

end

fprintf('First-derivative smoothng\n');

First-derivative smoothng

fprintf('Expected P=%d elements, got %d\n', N+2*(M-1), P);

Expected P=252 elements, got 252

fprintf('Expected L=%d rows, got %d\n', N+M-1, nr);

Expected L=152 rows, got 152

fprintf('\n');

% build sparse array

Fr=Fr(1:P); % truncate to actual length

Fc=Fc(1:P);

Fv=Fv(1:P);

gdaFsparse = sparse(Fr,Fc,Fv,L,M);

clear Fr Fc Fv; % delete no longer needed arrays

% least squares solution, using bicg()

tol = 1e-12; % tolerance

maxit = 3*(M+N); % maximum number of iterations allowed

FTf=gda_FTFrhs(f);

mest=bicg(@gda_FTFmul,FTf,tol ,maxit);

bicg converged at iteration 72 to a solution with relative residual 8.3e-13.

% plot

figure();

clf;

hold on;

set(gca,'LineWidth',2);set(gca,'LineWidth',3);

set(gca,'FontSize',14);hold on;

axis( [0, zmax, -1.1, 1.1 ] );

plot( z, mtrue, 'r-', 'LineWidth', 4 );

plot( z, mest, 'g-', 'LineWidth', 3 );

plot( zobs, dobs, 'ko', 'LineWidth', 2 );

xlabel('z');

ylabel('m(z)');

fprintf('Caption: First derivative smoothing. True solution (red), estimated\n');

Caption: First derivative smoothing. True solution (red), estimated

fprintf('solution (green), observed data (black)\n');

solution (green), observed data (black)

% PART 2: Second derivative smoothing

% build element arrays

L=N+M;

P=0;

Pmax = L*3*(M-2)+2*2+100;

Fr = zeros(Pmax,1);

Fc = zeros(Pmax,1);

Fv = zeros(Pmax,1);

f = zeros(L,1);

% the N data equations are just m = dobs. The only trick is lining

% up the corresponding elements of m and dobs, since they are not of

% the same length

nr=0; % row counter

for i = (1:N)

j=index(i);

nr=nr+1; % increment row counter

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = j; % column index

Fv(P,1) = 1.0; % value

f(nr,1) = dobs(i,1); % r.h.s. is data

end

% now implement 2nd derivative smoothness constraints in the

% interior of m (that makes M-2 rows)

epsilon=1.0;

rDz = 1/Dz;

rDz2 = 1/Dz^2;

for i = (1:M-2)

j=i;

nr=nr+1; % increment row counter

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = j; % column index

Fv(P,1) = rDz2; % value

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = j+1; % column index

Fv(P,1) = -2.0*rDz2; % value

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = j+2; % column index

Fv(P,1) = rDz2; % value

f(nr,1) = 0.0; % r.h.s. is data

end

% now implement 1st derivative smoothness constraints at the

% ends of m (that makes 2 rows)

% top row

nr=nr+1; % increment row counter

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = 1; % column index

Fv(P,1) = -rDz; % value

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = 2; % column index

Fv(P,1) = rDz; % value

% bottom row

nr=nr+1; % increment row counter

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = M-1; % column index

Fv(P,1) = -rDz; % value

P=P+1; % increment element number

Fr(P,1) = nr; % row index

Fc(P,1) = M; % column index

Fv(P,1) = rDz; % value

fprintf('Second-derivative smoothng\n');

Second-derivative smoothng

fprintf('Expected P=%d elements, got %d\n', N+3*(M-2)+2*2, P);

Expected P=353 elements, got 353

fprintf('Expected L=%d rows, got %d\n', N+M, nr);

Expected L=153 rows, got 153

fprintf('\n');

% build sparse array

Fr=Fr(1:P);

Fc=Fc(1:P);

Fv=Fv(1:P);

gdaFsparse = sparse(Fr,Fc,Fv,L,M);

clear Fr Fc Fv; % delete no longer needed arrays

% least squares solution, using bicg()

tol = 1e-12; % tolerance

maxit = 3*(M+N); % maximum number of iterations allowed

FTf = gda_FTFrhs(f);

mest=bicg(@gda_FTFmul,FTf,tol,maxit);

bicg converged at iteration 123 to a solution with relative residual 6.4e-13.

% plot

figure();

clf;

hold on;

set(gca,'LineWidth',2);set(gca,'LineWidth',3);

set(gca,'FontSize',14);hold on;

axis( [0, zmax, -1.1, 1.1 ] );

plot( z, mtrue, 'r-', 'LineWidth', 4 );

plot( z, mest, 'g-', 'LineWidth', 3 );

plot( zobs, dobs, 'ko', 'LineWidth', 2 );

xlabel('z');

ylabel('m(z)');

fprintf('Caption: Second derivative smoothing. True solution (red), estimated\n');

Caption: Second derivative smoothing. True solution (red), estimated

fprintf('solution (green), observed data (black)\n');

solution (green), observed data (black)

% gdama04_12 Example of smoothing operators

% using the formula of Menke and Eilon (2015)

clear all;

fprintf('gdama04_12\n');

gdama04_12

% z-axis

N=101;

Dx = 0.1;

xmin = -N*Dx/2;

x = xmin + Dx*[0:N-1]';

figure(1);

clf;

subplot(2,1,1);

set( gca, 'LineWidth', 3);

hold on;

ep = 0.3;

epi = 1/ep;

m1 = (epi/2)*exp(-epi*abs(x));

amp = max(abs(m1));

axis( [x(1), x(end), -1.1*amp, 1.1*amp ] );

plot( [x(1), x(end)]', [0,0]', 'b:', 'LineWidth', 3 );

plot( x, m1, 'k-', 'LineWidth', 3 );

ep = 0.6;

epi = 1/ep;

m1 = (epi/2)*exp(-epi*abs(x));

plot( x, m1, 'r-', 'LineWidth', 3 );

xlabel('x');

ylabel('m');

subplot(2,1,2);

set( gca, 'LineWidth', 3);

hold on;

ep = 0.3;

epi = 1/ep;

a = sqrt(2*ep);

V = (a^3)/(8*ep*ep);

m2 = V*exp(-abs(x)/a).*(cos(abs(x)/a)+sin(abs(x)/a));

amp = max(abs(m2));

axis( [x(1), x(end), -1.1*amp, 1.1*amp ] );

plot( [x(1), x(end)]', [0,0]', 'b:', 'LineWidth', 3 );

plot( x, m2, 'k-', 'LineWidth', 3 );

ep = 0.6;

epi = 1/ep;

a = sqrt(2*ep);

V = (a^3)/(8*ep*ep);

m2 = V*exp(-abs(x)/a).*(cos(abs(x)/a)+sin(abs(x)/a));

plot( x, m2, 'r-', 'LineWidth', 3 );

xlabel('x');

ylabel('m');

fprintf('Caption: Smoothing operators. (Top) First derivative, showing narrow (black)\n');

Caption: Smoothing operators. (Top) First derivative, showing narrow (black)

fprintf('wide (red). (Bottom) Second derivative, showing narrow (black) and wide (red).\n');

wide (red). (Bottom) Second derivative, showing narrow (black) and wide (red).

% gdam04_13

% constrained least squares fit to synthetic data

clear all;

% z's

N=30;

zmin=0;

zmax=10;

z = sort(random('Uniform',zmin,zmax,N,1));

% d = a + b*z + noise

a=2.0;

b=1.0;

sd=0.5;

dobs = a+b*z+random('Normal',0,sd,N,1);

% constraint

zp = 8;

dp = 6;

% constrained least squares fit

M=2;

G=[ones(N,1), z];

F = [1, zp];

A = [ [G'*G, F']', [F, 0]' ];

b = [ [G'*dobs]', dp ]';

mest = A\b;

mest=mest(1:M);

% predicted data

dpre = G*mest;

e = dobs - dpre;

[emax, iemax] = max(abs(e));

% plot

% plot scale

pdmin=0;

pdmax=15;

% plot observed and presicted data, constained point

figure();

clf;

hold on;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

axis( [zmin, zmax, pdmin, pdmax ]' );

plot( z, dobs, 'ro', 'LineWidth', 2);

plot( zp, dp, 'bo', 'LineWidth', 4);

plot( z, dpre, 'g-', 'LineWidth', 2);

xlabel('z');

ylabel('d');

fprintf('Caption: Constrained least squares. Observed data (red), predicted data (green),\n');

Caption: Constrained least squares. Observed data (red), predicted data (green),

fprintf('constraint (blue)');

constraint (blue)

% Weighted least squares fit as a check.

% The idea is to set the weight on the

% prior information to a very large number

% so that it functions as a hard constraint.

wd = 1;

wh = 1e6;

FF = [G*wd; F*wh ];

ff = [dobs*wd; dp*wh];

mest2 = (FF'*FF)\(FF'*ff);

dpre2 = G*mest2;

figure();

clf;

hold on;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

axis( [zmin, zmax, pdmin, pdmax ]' );

plot( z, dobs, 'ro', 'LineWidth', 2);

plot( zp, dp, 'bo', 'LineWidth', 4);

plot( z, dpre, 'g-', 'LineWidth', 3);

plot( z, dpre2, 'y--', 'LineWidth', 1);

xlabel('z');

ylabel('d');

fprintf('Caption: Weighted least squares (with strong weighting). Observed data (red),\n');

Caption: Weighted least squares (with strong weighting). Observed data (red),

fprintf('predicted data (green), predicted data from weighted soln (yelloe), constraint (blue)');

predicted data (green), predicted data from weighted soln (yelloe), constraint (blue)

% Explicit version of solution

GTG=G'*G;

GTd=G'*dobs;

mls = GTG\GTd;

Z = F*(GTG\F');

Dh = dp-F*mls;

q = Z\Dh;

mex = mls + GTG\(F'*q);

D = [mest,mex,mest-mex];

disp('comparison with explicit solution');

comparison with explicit solution

disp('mest mex mest-mes');

mest  mex  mest-mes

disp(D);

    3.8226    3.8226   -0.0000
    0.2722    0.2722    0.0000

% gdama04_14

% two examples of variance

clear all;

fprintf('gdama04_14\n');

gdama04_14

% z-axis

N=101;

zmin=0;

zmax=100;

Dz=(zmax-zmin)/(N-1);

z=zmin+Dz*[0:N-1]';

% number of model parameters = number of model data

M=N;

% true model all ones

mtrue = ones(M,1);

% two different problems

% 1: sum of first K blocks

% 2: sum of pairs of blocks

G1=toeplitz( ones(N,1), [1, zeros(1,N-1)] );

G2=toeplitz( [1, 1, zeros(1, N-2)]', [1, zeros(1,N-1)] );

d1true = G1*mtrue;

d2true = G2*mtrue;

sd=0.1;

% synthetic observed data

d1obs = d1true + random('Normal',0,sd,N,1);

d2obs = d2true + random('Normal',0,sd,N,1);

% least squares solutions

m1est = (G1'*G1)\(G1'*d1obs);

m2est = (G2'*G2)\(G2'*d2obs);

% variances

C1 = (sd^2) * inv(G1'*G1);

sm1 = sqrt(diag(C1));

C2 = (sd^2) * inv(G2'*G2);

sm2 = sqrt(diag(C2));

% plot

figure();

clf;

% plot scale

pmmin=-3;

pmmax=3;

psmin=0;

psmax=1;

% plot solutions

subplot(2,1,1);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, pmmin, pmmax ]' );

plot( z, mtrue, 'k-', 'LineWidth', 2 );

ylabel('m');

xlabel('z');

plot( [zmin, zmax ]', [0, 0]', 'k-', 'LineWidth', 2 );

plot( z, m2est, 'b-', 'LineWidth', 3 );

plot( z, m1est, 'k-', 'LineWidth', 3 );

plot( z, m1est, 'r-', 'LineWidth', 3 );

% plot sqrt(variance)

subplot(2,1,2);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

hold on;

axis( [zmin, zmax, psmin, psmax ]' );

plot( z, sm1, 'r-', 'LineWidth', 3 );

plot( z, sm2, 'b-', 'LineWidth', 3 );

ylabel('sm');

xlabel('z');

fprintf('Caption (Top) Predicted models 1 (red) and blue (2). (Bottom) Posterior\n');

Caption (Top) Predicted models 1 (red) and blue (2). (Bottom) Posterior

fprintf('variance of model 1 (red) and 2 (blue)\n');

variance of model 1 (red) and 2 (blue)

% gdama04_15

% Examine error surface of intercept and slope in a straight

% line fit to synthetic data

clear all;

fprintf('gdama04_15\n')'

gdama04_15

ans = 11

% z's

N=30;

zmin=-5;

zmax=5;

z = sort(random('Uniform',zmin,zmax,N,1));

% PART 1: data evenly spread along interval

% d = a + b*z + noise

a=2.0;

b=1.0;

sd=0.51;

dobs = a+b*z+random('Normal',0,sd,N,1);

% grid

Na=101;

amin=0;

amax=4;

Da = (amax-amin)/(Na-1);

a = amin + Da*[0:Na-1]';

Nb=101;

bmin=0;

bmax=4;

Db = (bmax-bmin)/(Nb-1);

b = bmin + Db*[0:Nb-1]';

% populate grid with errors

EA = zeros(Na,Nb);

for i=(1:Na)

for j=(1:Nb)

ao = amin+Da*(i-1);

bo = bmin+Db*(j-1);

dpre = ao + bo*z;

e = dobs-dpre;

EA(i,j)=e'*e;

end

% find minimum error

[Ep, p] = min(EA);

[EAmin, c1] = min(Ep);

r1 = p(c1);

a1 = amin+Da*(r1-1);

b1 = bmin+Db*(c1-1);

dpre = a1 + b1*z;

% covariance calculation

% least squares formula

G = [ones(N,1), z];

C1 = (sd^2)*inv(G'*G);

disp('case 1: data straddles origin');

case 1: data straddles origin

fprintf(' method 1: sm1 %f sm2 %f cov %f\n', sqrt(C1(1,1)), sqrt(C1(2,2)), C1(1,2) );

    method 1: sm1 0.093704  sm2 0.035093 cov -0.000369

% curvature of error surface formula

j=1;

d2Eda2 = (EA(r1+j,c1)-2*EA(r1,c1)+EA(r1-j,c1))/((j*Da)^2);

d2Edb2 = (EA(r1,c1+j)-2*EA(r1,c1)+EA(r1,c1-j))/((j*Db)^2);

d2Edadb = (EA(r1+j,c1+j)-EA(r1+j,c1-j)-EA(r1-j,c1+j)+EA(r1-j,c1-j))/(4*j*Da*j*Db);

DA=zeros(2,2);

DA(1,1)=d2Eda2;

DA(1,2)=d2Edadb;

DA(2,1)=d2Edadb;

DA(2,2)=d2Edb2;

C2 = (sd^2)*inv(DA/2);

fprintf(' method 2: sm1 %f sm2 %f cov %f\n', sqrt(C2(1,1)), sqrt(C2(2,2)), C2(1,2) );

    method 2: sm1 0.093704  sm2 0.035093 cov -0.000369

disp(' ');

% plot data and fit

figure(1);

clf;

subplot(1,2,1);

hold on;

subplot(1,2,1);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

axis( [zmin, zmax, -10, 10] );

plot( z, dobs, 'ro', 'LineWidth', 3);

plot( z, dpre, 'b-', 'LineWidth', 2);

xlabel('z');

ylabel('d(z)');

% plot error surface

figure(2);

clf;

subplot(1,2,1);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

hold on;

axis ij;

axis equal;

axis( [bmin, bmax, amin, amax] );

imagesc( [bmin, bmax], [amin, amax], EA );

plot( b1, a1, 'wo', 'LineWidth', 2);

EAmax=max(max(EA));

[X, Y] = meshgrid(b,a);

DeltaE=(EAmax-EAmin)/100;

contour(X,Y,EA,[EAmin+DeltaE, EAmin+DeltaE],'w-','LineWidth',3);

colorbar;

% PART 2: data bunched up at end of interval

% z's

N=30;

zmin=-5;

zmax=5;

z = sort(random('Uniform',zmin+1*(zmax-zmin)/2,zmax,N,1));

% d = a + b*z + noise

a=2.0;

b=1.0;

sd=0.5;

dobs = a+b*z+random('Normal',0,sd,N,1);

% grid

Na=101;

amin=0;

amax=4;

Da = (amax-amin)/(Na-1);

a = amin + Da*[0:Na-1]';

Nb=101;

bmin=0;

bmax=4;

Db = (bmax-bmin)/(Nb-1);

b = bmin + Db*[0:Nb-1]';

% populate grid with errors

EB = zeros(Na,Nb);

for i=(1:Na)

for j=(1:Nb)

ao = amin+Da*(i-1);

bo = bmin+Db*(j-1);

dpre = ao + bo*z;

e = dobs-dpre;

EB(i,j)=e'*e;

end

% find minimum error

[Ep, p] = min(EB);

[EBmin, c1] = min(Ep);

r1 = p(c1);

a1 = amin+Da*(r1-1);

b1 = bmin+Db*(c1-1);

dpre = a1 + b1*z;

% covariance calculation

% least squares formula

G = [ones(N,1), z];

C1 = (sd^2)*inv(G'*G);

disp('case 2: data to one side of origin');

case 2: data to one side of origin

fprintf(' method 1: sm1 %f sm2 %f cov %f\n', sqrt(C1(1,1)), sqrt(C1(2,2)), C1(1,2) );

    method 1: sm1 0.228829  sm2 0.068425 cov -0.014358

% curvature of error surface formula

j=1;

d2Eda2 = (EB(r1+j,c1)-2*EB(r1,c1)+EB(r1-j,c1))/((j*Da)^2);

d2Edb2 = (EB(r1,c1+j)-2*EB(r1,c1)+EB(r1,c1-j))/((j*Db)^2);

d2Edadb = (EB(r1+j,c1+j)-EB(r1+j,c1-j)-EB(r1-j,c1+j)+EB(r1-j,c1-j))/(4*j*Da*j*Db);

DB=zeros(2,2);

DB(1,1)=d2Eda2;

DB(1,2)=d2Edadb;

DB(2,1)=d2Edadb;

DB(2,2)=d2Edb2;

C2 = (sd^2)*inv(DB/2);

fprintf(' method 2: sm1 %f sm2 %f cov %f\n', sqrt(C2(1,1)), sqrt(C2(2,2)), C2(1,2) );

    method 2: sm1 0.228829  sm2 0.068425 cov -0.014358

disp(' ');

% plot data and fit

figure(1);

subplot(1,2,2);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

hold on;

axis( [zmin, zmax, -10, 10] );

plot( z, dobs, 'ro', 'LineWidth', 3);

plot( z, dpre, 'b-', 'LineWidth', 2);

fprintf('Caption: Linear fit (blue) to data (red). (left) Full dataset. (Right) data for z>0, ony')

Caption: Linear fit (blue) to data (red). (left) Full dataset. (Right) data for z>0, ony

% plot error surface

figure(2)

subplot(1,2,2);

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

hold on;

axis ij;

axis equal;

axis( [bmin, bmax, amin, amax] );

imagesc( [bmin, bmax], [amin, amax], EB );

plot( b1, a1, 'wo', 'LineWidth', 2);

EBmax=max(max(EB));

[X, Y] = meshgrid(b,a);

% note: use same DeltaE as previous plot

contour(X,Y,EB,[EBmin+DeltaE, EBmin+DeltaE],'w-','LineWidth',3);

colorbar;

fprintf('Caption: Linear fit (blue) to data (red). (left) Error surface of full dataset.\n');

Caption: Linear fit (blue) to data (red). (left) Error surface of full dataset.

fprintf('(Right) Errpr for data with z>0, only');

(Right) Errpr for data with z>0, only