% gdama06
% 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-05, Edit and integration with text
% gdama06_01
% data point inside of 3D box with 1:1:1 refernce line
clear all;
fprintf('gdama06_01\n');
% x-axis
Nx = 31;
xmin = 0;
xmax = 5;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
% y-axis
Ny = 31;
ymin = 0;
ymax = 5;
Dy = (ymax-ymin)/(Ny-1);
y = ymin + Dy*[0:Ny-1]';
% z-axis
Nz = 31;
zmin = 0;
zmax = 5;
Dz = (zmax-zmin)/(Nz-1);
z = zmin + Dz*[0:Nz-1]';
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% (x,y,z) grid
[XX, YY, ZZ] = meshgrid( x, y, z );
% the point is at a randomly chosen point in d-space
% using a Normal distribution with specified mean and variance
rbar = random('Normal',2.5,0.5,3,1);
% visualize the point in 3D using an isosurface of a normal distribution
sd=0.5;
C = (sd^2)*eye(3,3);
CI = inv(C);
DC = det(C);
norm = ( ((2*pi)^(3/2)) * sqrt(DC) );
% normal distribution
PP = zeros(Nx,Ny,Nz);
for i=(1:Nx)
for j=(1:Ny)
for k=(1:Nz)
r = [XX(i,j,k), YY(i,j,k), ZZ(i,j,k)]';
PP(i,j,k) = exp(-0.5*(r-rbar)'*CI*(r-rbar))/norm;
end
end
end
% get max for scale
maxP=max(max(max(PP)));
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [0, 5], [0, 5], [0, 5], 'b-', 'LineWidth', 3 );
% plot isosurface to show point in 3-space
p=patch(isosurface( XX, YY, ZZ, PP, 0.9*maxP ));
isonormals(XX,YY,ZZ,PP,p)
set(p, 'FaceColor', 'black', 'EdgeColor', 'none');
% set view direction & lighting
daspect([1 1 1])
view(3)
camlight; lighting phong
fprintf('Caption: Data point [d1,d2,d3] (circle) near the d2=d2=d3 line(blue)\n');
% gdama06_02
% isosurfaace of a 3D Normal distribution
clear all;
fprintf('gdama06_02\n');
% x0axis
Nx = 31;
xmin = 0;
xmax = 5;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
% y-axis
Ny = 31;
ymin = 0;
ymax = 5;
Dy = (ymax-ymin)/(Ny-1);
y = ymin + Dy*[0:Ny-1]';
% z-axis
Nz = 31;
zmin = 0;
zmax = 5;
Dz = (zmax-zmin)/(Nz-1);
z = zmin + Dz*[0:Nz-1]';
% (x,y,z) grid
[XX, YY, ZZ] = meshgrid( x, y, z );
% parameters for Normal distribution
rbar = [2.5, 2.5, 2.5]';
sd=0.5;
C = (sd^2)*eye(3,3);
CI = inv(C);
DC = det(C);
norm = ( ((2*pi)^(3/2)) * sqrt(DC) );
% normal distribution in 3D space
PP = zeros(Nx,Ny,Nz);
for i=(1:Nx)
for j=(1:Ny)
for k=(1:Nz)
r = [XX(i,j,k), YY(i,j,k), ZZ(i,j,k)]';
PP(i,j,k) = exp(-0.5*(r-rbar)'*CI*(r-rbar))/norm;
end
end
end
% for test purposes
% A = Dx*Dx*Dz*sum(sum(sum(PP)));
% max, for plotting purposes
maxP=max(max(max(PP)));
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [0, 5], [0, 5], [0, 5], 'b-', 'LineWidth', 3 );
% draw speherical cloud in 3D space
for ip = (1:10)
p(ip)=patch(isosurface( XX, YY, ZZ, PP, ip*maxP/10 ));
isonormals(XX,YY,ZZ,PP, p(ip))
set(p(ip), 'FaceColor', 'red', 'FaceAlpha', ip/20, 'EdgeColor', 'none');
end
% set view angle and lighting
daspect([1 1 1])
view(3)
camlight; lighting phong
fprintf('Caption: Depiction of Normal pdf (red) centered on d1=d2=d3 line\n');
% gdama06_03
% Likelihood surface for mean/variance
clear all;
fprintf('gdama06_03\n');
% random data
N = 100;
d = random('Normal',2.5,1.5,N,1);
% m1-axis, with m1=mean
Nm = 51;
mmin = 1.5;
mmax = 5.0;
Dm = (mmax-mmin)/(Nm-1);
m = mmin + Dm*[0:Nm-1]';
% s-axis, with s=sqrt(variance)
Ns = 51;
smin = 1;
smax = 2;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% (m,s) grid
[X,Y]=meshgrid( m, s);
% likelihood surface
L=zeros(Nm,Ns);
% tabulate likelihood surface
% Normal P = (1/sqrt(2pi)) * (1/s) * exp( -0.5 (d-m)^2 / s^2 )
% L = -0.5log(2*pi) - log(s) - ( -0.5 (d-m)^2 / s^2 )
for i=(1:Nm)
for j=(1:Ns)
mp=X(i,j);
sp=Y(i,j);
L(i,j) = -N*log(2*pi)/2 - N*log(sp) - 0.5*sum((d-mp).^2)/(sp^2);
end
end
% find (i,j) of point of maximum
[tmp, itmp] = max(L);
[Lmax, Lj] = max(tmp);
Li=itmp(Lj);
% maximum likelihood point
mbest = X(Li,Lj);
sbest = Y(Li,Lj);
fprintf('mean %f sigma %f\n',mbest,sbest);
% ninimum for plotting purposes
Lmin = min(min(L));
% 3D plot
figure(1);
clf;
colormap('jet');
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
pmmin=mmin;
pmmax=mmax;
psmin=smin;
psmax=smax;
pLmin = -500;
pLmax = -100;
axis( [pmmin, pmmax, psmin, psmax, pLmin, pLmax]');
% clip parts of grid by setting L to NaN
for i=(1:Nm)
for j=(1:Ns)
if( L(i,j) < pLmin )
L(i,j)=NaN(1);
end
end
end
% draw the mesh
mesh(X,Y,L);
% sides of a 3D box
pxmin=pmmin+0.01; pxmax=pmmax-0.01;
pymin=psmin+0.01; pymax=psmax-0.01;
pzmin=pLmin+1; pzmax=pLmax-1;
% improvise 3D box on the plot
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 );
DL = (pLmax-pLmin)/20;
plot3( [mbest, mbest], [sbest, sbest], [Lmax-DL, Lmax+DL], 'r-', 'LineWidth', 4 );
plot3( [mbest-0.1, mbest+0.1], [sbest, sbest], [Lmax, Lmax], 'r-', 'LineWidth', 4 );
plot3( [mbest, mbest], [sbest-0.1, sbest+0.1], [Lmax, Lmax], 'r-', 'LineWidth', 4 );
% set view angle
view(3);
xlabel('mean');
ylabel('sigma');
zlabel('L');
colorbar();
fprintf('Caption; Depiction of likelihood surfaceL(mean,sigma)\n');
fprintf('with estimated model (red) superimposed\n');
% gdama06_04
% examples of probability density function p(d1,d2)
% one is peaked, the other had a ridge
clear all;
fprintf('gdama05_04\n');
% d1-axis
Nd1 = 51;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% d2-axis
Nd2 = 51;
d2min = 0;
d2max = 5.0;
Dd2 = (d2max-d2min)/(Nd2-1);
d2 = d2min + Dd2*[0:Nd2-1]';
% (d1,d2) grid
[X,Y]=meshgrid( d1, d2);
% distribution 1, has a peak
P1=zeros(Nd1,Nd2);
dbar = [2.5, 2.5]';
sd1 = 0.5;
sd2 = 0.5;
C = diag( [sd1, sd2]' );
CI = inv(C);
DC = det(C);
norm = (1/(2*pi)) * (1/sqrt(DC));
for i=(1:Nd1)
for j=(1:Nd2)
d =[X(i,j), Y(i,j)]' - dbar;
P1(i,j) = norm*exp( -0.5 * d'*CI* d );
end
end
% for test purposes
% A = Dd1*Dd2*sum(sum(P1)
% distribution 2, ridge
% Note distribution not normalizable
P2=zeros(Nd1,Nd2);
sda = 0.5;
d0a = [2.5, 2.5]';
for i=(1:Nd1)
for j=(1:Nd2)
dr =[X(i,j)-d0a(1), Y(i,j)-d0a(2)]';
r2= (dr'*[1,1]')^2;
P2(i,j) = 0.2*exp(-0.5*r2/sda);
end
end
% plot
figure();
clf;
hold on;
colormap('jet');
% improvise 3D box
pxmin=d1min; pxmax=d1max;
pymin=d2min; pymax=d2max;
pzmin=0; pzmax=0.4;
subplot(1,2,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [pxmin-0.01, pxmax+0.01, pymin-0.01, pymax+0.01, pzmin-0.01, pzmax+0.01 ]' );
% plot distribution P1
mesh(X,Y,P1);
% 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 );
% view angle
view(3);
subplot(1,2,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [pxmin-0.01, pxmax+0.01, pymin-0.01, pymax+0.01, pzmin-0.01, pzmax+0.01 ]' );
% plot distribution P2
mesh(X,Y,P2);
% 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 );
% view angle
view(3);
fprintf('Caption: (left) Pdf p(x1,x2) with a peak. (right) Pdf\n');
fprintf('p(x1,x2) with a ridge.\n');
% gdama06_05
% examples of probability density function p(m1,m2)
% uncorrelated, equal variance, one with small
% variance, the other with large variance
clear all;
fprintf('gdama06_05\n');
% m1-axis
Nm1 = 51;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% m2-axis
Nm2 = 51;
m2min = 0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = m2min + Dm2*[0:Nm2-1]';
% setup for distribution 1, small variance
P1=zeros(Nm1,Nm2);
mbar1 = [2.5, 2.5]';
sd1 = 0.5;
C1 = diag( [sd1^2, sd1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
% setup for distributions 2, large variance
P2=zeros(Nm1,Nm2);
mbar2 = [2.5, 2.5]';
sd2 = 1;
C2 = diag( [sd2^2, sd2^2]' );
CI2 = inv(C2);
DC2 = det(C2);
% make distributions
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
norm2 = (1/(2*pi)) * (1/sqrt(DC2));
for i=(1:Nm1)
for j=(1:Nm2)
x1 =[m1(i), m2(j)]' - mbar1;
x2 =[m1(i), m2(j)]' - mbar2;
P1(i,j) = norm1*exp( -0.5 * x1'*CI1* x1 );
P2(i,j) = norm2*exp( -0.5 * x2'*CI2* x2 );
end
end
% for test purposes
% A1 = Dd1*Dd2*sum(sum(P1));
% A2 = Dd1*Dd2*sum(sum(P2));
% draw distributions
gda_draw(' ', P1, ' ', ' ', P2 );
fprintf('Caption: (Left) Pdf p(m1,m2) with small and equal sigmas. (right)\n');
fprintf('(Right) Pdf p(m1,m2) with large and equal sigmas. (right)\n');
% gdama06_06
% examples of probability density function p(m1,m2)
% unequal variance, uncorrelated
clear all;
fprintf('gdama06_06\n');
% m1-axis
Nm1 = 51;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% m2-axis
Nm2 = 51;
m2min = 0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = m2min + Dm2*[0:Nm2-1]';
% setup for distribution 1, C11<C22
P1=zeros(Nm1,Nm2);
mbar1 = [2.5, 2.5]';
sd1 = 0.5;
sd2 = 1.0;
C1 = diag( [sd1^2, sd2^2]' );
CI1 = inv(C1);
DC1 = det(C1);
% normalization
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% build distribution
for i=(1:Nm1)
for j=(1:Nm2)
x1 =[m1(i), m2(j)]' - mbar1;
P1(i,j) = norm1*exp( -0.5 * x1'*CI1* x1 );
end
end
% for test purposes
% A1 = Dd1*Dd2*sum(sum(P1));
gda_draw(' ', P1 );
fprintf('Caption: Pdf p(m1,m2) with unequal sigmas and zero covarince.\n');
% gdama06_07
%
% examples of probability density function p(m1,m2)
% one is a ridge, the other is Normal but mimics a ridge
clear all;
fprintf('gdama06_07\n');
% m1-axis
Nm1 = 51;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% m2-axis
Nm2 = 51;
m2min = 0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = m2min + Dm2*[0:Nm2-1]';
% setup for distribution 1, Normal
P1=zeros(Nm1,Nm2);
mbar1 = [2.5, 2.5]';
sd12 = 2;
cov = (2-0.1);
C1 = [ [sd12, cov]', [cov, sd12]' ];
CI1 = inv(C1);
DC1 = det(C1);
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% calculate distribution
for i=(1:Nm1)
for j=(1:Nm2)
x1 =[m1(i), m2(j)]' - mbar1;
P1(i,j) = norm1*exp( -0.5 * x1'*CI1* x1 );
end
end
% setup for distribution 2, ridge
% Note distribution not normalizable
P2=zeros(Nm1,Nm2);
sda2 = 0.5;
d0a = [2.5, 2.5]';
for i=(1:Nm1)
for j=(1:Nm2)
dr =[m1(i)-d0a(1), m2(j)-d0a(2)]';
r2= (dr'*[1,-1]')^2;
P2(i,j) = 0.2*exp(-0.5*r2/sda2);
end
end
% plot distributions
gda_draw(' ', P2, ' ', ' ', P1 );
fprintf('Caption (left) Pdf with a true ridge. (right) Pdf with an elongated\n')
fprintf('peak approximating a ridge\n');
% gdama06_08
%
% example of probability density function p(m1,m2)
% implements inequality constraint m1<m2
clear all;
fprintf('gdama06_08');
% m1-axis
Nm1 = 51;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% m2-axis
Nm2 = 51;
m2min = 0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = m2min + Dm2*[0:Nm2-1]';
% compute distribution 1
P1=zeros(Nm1,Nm2);
for i=(1:Nm1)
for j=(1:Nm2)
P1(i,j) = (m1(i)<=m2(j));
end
end
% for test purposes
% A1 = Dd1*Dd2*sum(sum(P1));
% plot distribution
gda_draw(' ', P1 );
fprintf('Caption: Pdf p(m1,m2) implementing the inequality constraint m1<m2\n');
% gdama06_09
% relative entropy of onenormsl pdf p(m)
% with respect to another q(m)
clear all;
fprintf('gdama06_09\n');
% m-axis
N=101;
mmin = -25;
mmax = 25;
Dm = (mmax-mmin)/(N-1);
m = mmin + Dm*[0:N-1]';
% p(m) = pA(m)
mbarp = 0;
sigmamp = 1;
p = (1 /(sqrt(2*pi)*sigmamp))*exp(-((m-mbarp).^2)/(2*sigmamp*sigmamp));
normp = Dm*sum(p);
% q(m) = pN(m)
mbarq = 0;
sigmamq = 5;
q = (1 /(sqrt(2*pi)*sigmamq))*exp(-((m-mbarq).^2)/(2*sigmamq*sigmamq));
normq = Dm*sum(q);
% plot
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, max(p)] );
plot( m, p, 'r-', 'LineWidth', 3);
plot( m, q, 'g-', 'LineWidth', 3);
xlabel('m');
ylabel('p and q');
fprintf('Caption: Two Normal pdfs p(m) (red) and q(m) (green)\n');
% relative entropy S
r = (sigmamp^2)/(sigmamq^2);
Sanalytic = ((mbarp-mbarq)^2)/(2*sigmamq*sigmamq) + 0.5*(r-1-log(r));
Snumeric = Dm*sum( p .* log( p ./ q ) );
fprintf('S = %f (analytic) and %f (numeric)\n', Sanalytic, Snumeric );
% relative entropy as a function of sigmamp
Nv=100;
sigmampv = sigmamq*[1:Nv]'/Nv;
rv = (sigmampv.^2)/(sigmamq^2);
Sv = ((mbarp-mbarq)^2)./(2*sigmamq*sigmamq) + 0.5*(rv-1-log(rv));
% plot
figure(2);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, sigmamq, 0, max(Sv)] );
plot( sigmampv, Sv, 'k-', 'LineWidth', 3);
plot( sigmamp, Sanalytic, 'ro', 'LineWidth', 3);
xlabel('sigmamp');
ylabel('S');
fprintf('Caption: Relative entropy of p(m) and q(m) as the\n');
fprintf('as the sigma of p id varied. One point on curve (red circle)\n');
fprintf('determined by an analytic method\n');
% gdama06_10
%
% examples of probability density function p(d1,m1)
clear all;
fprintf('gdama06_10\n');
% d1-axis
Nd1 = 51;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% m1-axis
Nm1 = 51;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% setup for distribution
P1=zeros(Nd1,Nm1);
mbar1 = [2.5, 2.5]';
sd1 = 0.5;
sd2 = 1;
C1 = diag( [sd1^2, sd2^2]' );
CI1 = inv(C1);
DC1 = det(C1);
% normalization
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% tabulate distribution
for i=(1:Nd1)
for j=(1:Nm1)
x1 =[m1(i), m1(j)]' - mbar1;
P1(i,j) = norm1*exp( -0.5 * x1'*CI1* x1 );
end
end
% for test purposes
% A1 = Dd1*Dd2*sum(sum(P1));
% plot distribution
gda_draw(' ', P1 );
fprintf('Caption: Probability density function p(m,d)');
% gdama06_11
%
% parameteric function passing thru distribution p(d1,m1)
clear all;
% d1-axis
Nd1 = 101;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% setup for distribution
P1=zeros(Nd1,Nm1);
d1bar = 2.25;
m1bar = 2.08;
bar = [d1bar, m1bar]';
sd1 = 0.5;
sm1 = 1;
C1 = diag( [sd1^2, sm1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% tablulate distribution
for i=(1:Nm1)
for j=(1:Nd1)
x1 =[d1(i), m1(j)]' - bar;
P1(i,j) = norm1*exp( -0.5 * x1'*CI1*x1 );
end
end
% s-axis for parametric curve
Ns = 51;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% parameric curve
dp = 1+s-2*(s/smax).^2;
mp = s;
% P on parameteric curve
ipd = 1+floor((dp-d1min)/Dd1);
ipm = 1+floor((mp-m1min)/Dm1);
% insure indices in range
i=find(ipd<1);
ipd(i)=1;
i=find(ipd>Nd1);
ipd(i)=Nd1;
i=find(ipm<1);
ipm(i)=1;
i=find(ipm>Nm1);
ipm(i)=Nm1;
Ps=zeros(Ns,1);
% evaluate P at indices
for i = (1:Ns)
Ps(i) = P1(ipd(i), ipm(i));
end
% maximum along curve
[Pmax, ismax]=max(Ps);
d1smax = dp(ismax);
m1smax = mp(ismax);
% plot p(d1,m1) with parametric curve crossing it
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P1 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( mp, dp, 'w-', 'LineWidth', 3 );
plot( m1smax, d1smax, 'ko', 'LineWidth', 4);
plot( [m1min, m1smax], [d1smax, d1smax], 'w:', 'LineWidth', 2);
plot( [m1smax, m1smax], [d1max, d1smax], 'w:', 'LineWidth', 2);
xlabel('m');
ylabel('d');
fprintf('Caption: Probability density function p(d,m) (colors) with parametric curve\n');
fprintf('f(d,m)=0 (white) superimposed. The point of mximum probability\n');
fprintf('along the curve is shown with a black circle\n');
% plot distribution as a function of arclength 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)');
fprintf('Caption: Probability as a function of arc-length along the parametric curve.\n');
fprintf('The probability has a single maximum\n');
% gdama06_12
%
% parameteric function passing thru distribution p(d1,m2)
clear all;
fprintf('gdama06_12\n');
% d1-axis
Nd1 = 101;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% setup for distribution
P1=zeros(Nd1,Nm1);
d1bar = 2.25;
m1bar = 2.08;
bar = [d1bar, m1bar]';
sd1 = 1;
sm1 = 0.2;
C1 = diag( [sd1^2, sm1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% tabulate distribution
for i=(1:Nm1)
for j=(1:Nd1)
x1 =[d1(i), m1(j)]' - bar;
P1(i,j) = norm1*exp( -0.5 * x1'*CI1*x1 );
end
end
% s-axis for parametric curve
Ns = 51;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% parameric curve
dp = 1+s-2*(s/smax).^2;
mp = s;
% P on parameteric curve
ipd = 1+floor((dp-d1min)/Dd1);
ipm = 1+floor((mp-m1min)/Dm1);
% insure indices in range
i=find(ipd<1);
ipd(i)=1;
i=find(ipd>Nd1);
ipd(i)=Nd1;
i=find(ipm<1);
ipm(i)=1;
i=find(ipm>Nm1);
ipm(i)=Nm1;
Ps=zeros(Ns,1);
% evaluate P at indices
for i = [1:Ns]
Ps(i) = P1(ipd(i), ipm(i));
end
% maximum along curve
[Pmax, ismax]=max(Ps);
d1smax = dp(ismax);
m1smax = mp(ismax);
% plot distribution with parameteric curve superimposed
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P1 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( mp, dp, 'w-', 'LineWidth', 3 );
plot( m1smax, d1smax, 'ko', 'LineWidth', 4);
plot( [m1min, m1smax], [d1smax, d1smax], 'w:', 'LineWidth', 2);
plot( [m1smax, m1smax], [d1max, d1smax], 'w:', 'LineWidth', 2);
xlabel('m');
ylabel('d');
fprintf('Caption: Probability density function p(d,m) (colors) with parametric curve\n');
fprintf('f(d,m)=0 (white) superimposed. In this case the prior model is very certain\n');
fprintf('The point of mximum probability along the curve is shown with a black circle.\n');
% plot distribution as a function of arc-length along curve
figure(2);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [smin, smax, 0, 1] );
axis xy;
plot( s, Ps, 'b-', 'Linewidth', 3 );
xlabel('s');
ylabel('p(s)');
fprintf('Caption: Probability as a function of arc-length along the parametric curve.\n');
fprintf('The probability has a single maximum\n');
% gdama06_13
%
% parameteric function passing thru distribution p(d1,m1)
clear all;
fprintf('gdama06_13\n');
% d1-axis
Nd1 = 101;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% setup for distribution
P1=zeros(Nd1,Nm1);
d1bar = 2.25;
m1bar = 2.08;
bar = [d1bar, m1bar]';
sd1 = 0.2;
sm1 = 1.0;
C1 = diag( [sd1^2, sm1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
% tabulate distribution
for i=[1:Nm1]
for j=[1:Nd1]
x1 =[d1(i), m1(j)]' - bar;
P1(i,j) = norm1*exp( -0.5 * x1'*CI1*x1 );
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
dp = 1+s-2*(s/smax).^2;
mp = s;
% P on parameteric curve
ipd = 1+floor((dp-d1min)/Dd1);
ipm = 1+floor((mp-m1min)/Dm1);
% insure indices in range
i=find(ipd<1);
ipd(i)=1;
i=find(ipd>Nd1);
ipd(i)=Nd1;
i=find(ipm<1);
ipm(i)=1;
i=find(ipm>Nm1);
ipm(i)=Nm1;
Ps=zeros(Ns,1);
% evaluate P at indices
for i = (1:Ns)
Ps(i) = P1(ipd(i), ipm(i));
end
% maximum along curve
[Pmax, ismax]=max(Ps);
d1smax = dp(ismax);
m1smax = mp(ismax);
% plot distribution with parameteric curve superimposed
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P1 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( mp, dp, 'w-', 'LineWidth', 3 );
plot( m1smax, d1smax, 'ko', 'LineWidth', 4);
plot( [m1min, m1smax], [d1smax, d1smax], 'w:', 'LineWidth', 2);
plot( [m1smax, m1smax], [d1max, d1smax], 'w:', 'LineWidth', 2);
xlabel('m');
ylabel('d');
fprintf('Caption: Probability density function p(d,m) (colors) with parametric curve\n');
fprintf('f(d,m)=0 (white) superimposed. In this case the prior model is very uncertain\n');
fprintf('The point of mximum probability along the curve is shown with a black circle.\n');
% plot distribution as a function of arc-length along the curve
figure(2);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [smin, smax, 0, 1] );
axis xy;
plot( s, Ps, 'b-', 'Linewidth', 3 );
xlabel('s');
ylabel('p(s)');
fprintf('Caption: Probability as a function of arc-length along the parametric curve.\n');
fprintf('The probability has a single maximum\n');
% gdama06_14
%
% prior distribution pA(d1,m1) interacting with inexact
% theory pg(d1,m1) to produce total distribution pT(d1,m1)
clear all;
fprintf('gdama06_14\n');
% d1-axis
Nd1 = 101;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% setup for distribution P1=PA
P1=zeros(Nd1,Nm1);
d1bar = 2.25;
m1bar = 2.08;
bar = [d1bar, m1bar]';
sd1 = 0.5;
sm1 = 1;
C1 = diag( [sd1^2, sm1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
% note not normaalized, so max is unity
for i=(1:Nm1)
for j=(1:Nd1)
x1 =[d1(i), m1(j)]' - bar;
P1(i,j) = exp( -0.5 * x1'*CI1*x1 );
end
end
% s-axis for parametric curve s
Ns = 51;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% Pg follows parameric curve
dp = 1+s-2*(s/smax).^2;
mp = s;
% P2=Pg distribution; not normalizable
P2 = zeros(Nd1,Nm1);
sg = 0.35;
sg2 = sg^2;
for i=(1:Nm1)
for j=(1:Nd1)
r2 = (d1(i)-dp).^2 + (m1(j)-mp).^2;
r2min=min(r2);
P2(i,j) = exp( -r2min/sg2 );
end
end
% P3=PT is product of PA and Pg
P3 = P1.*P2;
[tmp, itmp] = max(P3);
[P3max, Pj] = max(tmp);
Pi=itmp(Pj);
Pmaxm1 = m1min+Dm1*(Pj-1);
Pmaxd1 = d1min+Dd1*(Pi-1);
figure(1);
clf;
% plot PA
subplot(1,3,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P1 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
xlabel('m');
ylabel('d');
% Plot Pg
subplot(1,3,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P2 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
xlabel('m');
ylabel('d');
% plot PT
subplot(1,3,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P3 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( Pmaxm1, Pmaxd1, 'ko', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( [Pmaxm1, Pmaxm1], [d1max, d1max-0.1], 'k-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [Pmaxd1, Pmaxd1], 'k-', 'LineWidth', 3 );
xlabel('m');
ylabel('d');
fprintf('Caption: (Left) The prior pdf pA combining observation and prior information\n');
fprintf('(Middle) the pdf of the theory Pg (right) the combined pdf pT\n');
% gdama05_15
%
% Two different cases of:
%
% prior distribution pA(d1,m1) interacting with inexact
% theory pg(d1,m1) to produce total distribution pT(d1,m1)
%
% once case with very exact theory, other with very inexact one
clear all;
fprintf('gdama05_15\n');
% d1-axis
Nd1 = 101;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% Case 1: nearly inexact theory
P1=zeros(Nd1,Nm1);
d1bar = 2.25;
m1bar = 2.08;
bar = [d1bar, m1bar]';
sd1 = 0.5;
sm1 = 1;
C1 = diag( [sd1^2, sm1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
% note not normalized, so max is unity
for i=(1:Nm1)
for j=(1:Nd1)
x1 =[d1(i), m1(j)]' - bar;
P1(i,j) = exp( -0.5 * x1'*CI1*x1 );
end
end
% s-axis for parametric curve
Ns = 51;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% parameric curve
dp = 1+s-2*(s/smax).^2;
mp = s;
% Parametric distribution; not normalizable
P2 = zeros(Nd1,Nm1);
sg = 0.1;
sg2 = sg^2;
for i=(1:Nm1)
for j=(1:Nd1)
r2 = (d1(i)-dp).^2 + (m1(j)-mp).^2;
r2min=min(r2);
P2(i,j) = exp( -r2min/sg2 );
end
end
P3 = P1.*P2;
[tmp, itmp] = max(P3);
[P3max, Pj] = max(tmp);
Pi=itmp(Pj);
Pmaxm1 = m1min+Dm1*(Pj-1);
Pmaxd1 = d1min+Dd1*(Pi-1);
figure(2);
clf;
% PA
subplot(1,3,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P1 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
xlabel('m');
ylabel('d');
% Pg, inexact
subplot(1,3,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P2 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
xlabel('m');
ylabel('d');
% PT (for inexact Pg)
subplot(1,3,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P3 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( Pmaxm1, Pmaxd1, 'ko', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( [Pmaxm1, Pmaxm1], [d1max, d1max-0.1], 'k-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [Pmaxd1, Pmaxd1], 'k-', 'LineWidth', 3 );
xlabel('m');
ylabel('d');
fprintf('Caption: Prior pdf pA(d1,m1) (left) interacting with\n');
fprintf('inexact theory pg(d1,m1) (middle) to produce total pdf pT(d1,m1)\n');
fprintf('(right). In this case, the theory is very exact\n');
% Case 2: exact (narrow) theory
P1=zeros(Nd1,Nm1);
d1bar = 2.25;
m1bar = 2.08;
bar = [d1bar, m1bar]';
sd1 = 0.5;
sm1 = 1;
C1 = diag( [sd1^2, sm1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
% note not normalized, so max is unity
for i=(1:Nm1)
for j=(1:Nd1)
x1 =[d1(i), m1(j)]' - bar;
P1(i,j) = exp( -0.5 * x1'*CI1*x1 );
end
end
% s-axis for parametric curve
Ns = 51;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% parameric curve
dp = 1+s-2*(s/smax).^2;
mp = s;
% Parametric distribution; not normalizable
P2 = zeros(Nd1,Nm1);
sg = 2;
sg2 = sg^2;
for i=(1:Nm1)
for j=(1:Nd1)
r2 = (d1(i)-dp).^2 + (m1(j)-mp).^2;
r2min=min(r2);
P2(i,j) = exp( -r2min/sg2 );
end
end
% PT
P3 = P1.*P2;
[tmp, itmp] = max(P3);
[P3max, Pj] = max(tmp);
Pi=itmp(Pj);
Pmaxm1 = m1min+Dm1*(Pj-1);
Pmaxd1 = d1min+Dd1*(Pi-1);
figure(1);
set( gcf, 'pos', [10, 10, 800, 300] );
clf;
% Plot PA, exact
subplot(1,3,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P1 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
xlabel('m');
ylabel('d');
% Plot Pg (exact)
subplot(1,3,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P2 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
xlabel('m');
ylabel('d');
% PT (exact)
subplot(1,3,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P3 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( Pmaxm1, Pmaxd1, 'ko', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( [Pmaxm1, Pmaxm1], [d1max, d1max-0.1], 'k-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [Pmaxd1, Pmaxd1], 'k-', 'LineWidth', 3 );
xlabel('m');
ylabel('d');
fprintf('Caption: Prior pdf pA(d1,m1) (left) interacting with\n');
fprintf('inexact theory pg(d1,m1) (middle) to produce total pdf pT(d1,m1)\n');
fprintf('(right). In this case, the theory is very inexact\n');
% gdama06_16
%
% total distribution pT(m1,d1) and the projected
% distribution p(m) = (integral) pT(m1,d1) d(d1)
% highlighting the possibility that the maximum
% liklihood points of the two distributions are
% at two different values of m1
clear all;
fprintf('gdama06_16\n');
% d1-axis
Nd1 = 101;
d1min = 0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = d1min + Dd1*[0:Nd1-1]';
% m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = m1min + Dm1*[0:Nm1-1]';
% distribution P1 = pA
P1=zeros(Nd1,Nm1);
d1bar = 2.25;
m1bar = 2.08;
bar = [d1bar, m1bar]';
sd1 = 0.5;
sm1 = 1;
C1 = diag( [sd1^2, sm1^2]' );
CI1 = inv(C1);
DC1 = det(C1);
% note not normaalized, so max is unity
for i=(1:Nm1)
for j=(1:Nd1)
x1 =[d1(i), m1(j)]' - bar;
P1(i,j) = exp( -0.5 * x1'*CI1*x1 );
end
end
% s-axis for parametric curve
Ns = 51;
smin = 0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = smin + Ds*[0:Ns-1]';
% parameric curve
dp = 1+s-2*(s/smax).^2;
mp = s;
% P2 = Pg
P2 = zeros(Nd1,Nm1);
sg = 0.35;
sg2 = sg^2;
for i=(1:Nm1)
for j=(1:Nd1)
r2 = (d1(i)-dp).^2 + (m1(j)-mp).^2;
r2min=min(r2);
P2(i,j) = exp( -r2min/sg2 );
end
end
% P3 = pT with pT=pA pg
P3 = P1.*P2;
[tmp, itmp] = max(P3);
[P3max, Pj] = max(tmp);
Pi=itmp(Pj);
Pmaxm1 = m1min+Dm1*(Pj-1);
Pmaxd1 = d1min+Dd1*(Pi-1);
% find maximum likelihood point
Pm = sum(P3);
[Pmmax, Pmi] = max(Pm);
Pmm = m1min+Dm1*(Pmi-1);
norm=Dm1*sum(Pm);
Pm = Pm/norm;
figure();
clf;
% plot pT(d1,m1)
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [d1min, d1max, m1min, m1max] );
axis ij;
imagesc( [d1min, d1max], [m1min, m1max], P3 );
plot( mp, dp, 'w:', 'LineWidth', 2 );
plot( m1bar, d1bar, 'wo', 'LineWidth', 3 );
plot( Pmaxm1, Pmaxd1, 'ko', 'LineWidth', 3 );
plot( [m1bar, m1bar], [d1max, d1max-0.1], 'w-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [d1bar, d1bar], 'w-', 'LineWidth', 3 );
plot( [Pmaxm1, Pmaxm1], [d1max, d1max-0.1], 'k-', 'LineWidth', 3 );
plot( [m1min, m1min+0.1], [Pmaxd1, Pmaxd1], 'k-', 'LineWidth', 3 );
xlabel('m');
ylabel('d');
fprintf('Caption: Total pdf pT(d,m) with prior value (white circle), maximum\n');
fprintf('likelihood point (black cirle), and theory (dashed line) superimposed.\n');
% plot p(m1)
figure();
clf;
subplot(2,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [m1min, m1max, 0, 1] );
axis xy;
plot( m1, Pm', 'b-', 'LineWidth', 3 );
plot( Pmm, max(Pm), 'ko', 'LineWidth', 3 );
xlabel('m');
ylabel('p(m)');
fprintf('Caption: Projected pdf p(m)=(integral) pT(d,m) dd (blue curve) and\n');
fprintf('maximum likelihood point (black cirle).\n');
% gdama06_17
% sample chi-squared analysis
% distribution of errors P (Phi), E and L
% for a problem of a band-limited timeseries
% with redundant noisy observations of the
% time series and with a priori smoothness
% information
clear all;
fprintf('gdama06_17\n');
% time t-axis
M=101;
Dt=1;
t = Dt*[0:M-1]';
% make a wiggly curve m(t) by starting out with random noise
% and bandpass filtering it around a narrow range of angular
% frquencies centered on w0. The bandpass filtering is
% accomplished by the gda_chebyshevfilt() function, which
% though not mentioned in the text (sorry about that) is pretty
% easy to use and has many other applications. Its arguments are
% output_timeseries = gda_chebyshevfilt( input_timeseries, Dt, f1, f2 )
% where Dt is the sampling of the timeseries (say in s) and
% (f1,f2) are the (low,high) size of the frequency band (say in Hz)
n = 10;
A = 2.0;
w0 = n*pi/t(end);
Dw = w0/10;
mtrue = gda_chebyshevfilt( random('Normal',0,1,M,1), Dt, (w0-Dw)/(2*pi), (w0+Dw)/(2*pi) );
mtrue = A*mtrue/max(abs(mtrue));
% components of data kernel, all normalize to produce
% data of about the same value
% data kernel is just 3 independent observations of same point
N = 3*M;
G = [eye(M,M); eye(M,M); eye(M,M)];
dtrue = G*mtrue;
sigma_d = 0.1;
dobs = dtrue + random('Normal',0,sigma_d,N,1);
sqrtcovdi = eye(N,N)/sigma_d;
% prior information: second derivative close to zero
H = (1/Dt^2)*toeplitz( [1; zeros(M-3,1)], [1, -2, 1, zeros(1,M-3) ] );
h = zeros(M-2,1);
K = M-2;
mddtrue = H*mtrue;
sigma_mdd = (w0^2)*A/sqrt(2);
sigma_mdd2 = std(mddtrue); % for test purposes
sqrtcovhi = eye(K,K)/sigma_mdd;
% overall least squares equation and its solution
F = [ sqrtcovdi*G; sqrtcovhi*H ];
f = [ sqrtcovdi*dobs; sqrtcovhi*h ];
FTFinv = inv(F'*F);
mest = FTFinv*(F'*f);
mddest = H*mest;
sigma_mdd3 = std(mddest); % for test purposes
% plot m
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, t(end), -1.1*A, 1.1*A] );
plot( t, dobs(1:M), 'b.', 'LineWidth', 2 );
plot( t, dobs(M+1:2*M), 'b.', 'LineWidth', 2 );
plot( t, dobs(2*M+1:3*M), 'b.', 'LineWidth', 2 );
plot( t, mtrue, 'k-', 'LineWidth', 4 );
plot( t, mest, 'r-', 'LineWidth', 2 );
xlabel('t (s)');
ylabel('m');
fprintf('Caption: Bandlimited time series, m(t). True (black),\n');
fprintf('estimated (red), redundsnt observations (blue)\n');
% chi-squared values
e = sqrtcovdi*(G*mest-dobs); % normalized prediction error
l = sqrtcovhi*(H*mest-h); % normalized prior error
E = e'*e; % Eobs
L = l'*l; % Lobs
P = E+L; % PHIobs
vP = N+K-M; % degrees of freeedom in PHI
vE = vP*N/(N+K); % approx degrees of freeedom in E
vL = vP*K/(N+K); % approx degrees of freeedom in L
% Chi-squared test of the Null Hypothesis
fprintf('Analytic 95%% confidence tests\n');
if ( P>(vP-2*sqrt(2*vP)) && P<(vP+2*sqrt(2*vP)) )
fprintf('P %.2f vP %.2f null hypothesis cannot be rejected\n', P, vP);
else
fprintf('P %.2f vP %.2f null hypothesis can be rejected\n', P, vP);
end
if ( E>(vE-2*sqrt(2*vE)) && E<(vE+2*sqrt(2*vE)) )
fprintf('E %.2f vE %.2f null hypothesis cannot be rejected\n', E, vE);
else
fprintf('E %.2f vE %.2f null hypothesis can be rejected\n', E, vE);
end
if ( L>(vL-2*sqrt(2*vL)) && L<(vL+2*sqrt(2*vL)) )
fprintf('L %.2f vL %.2f null hypothesis cannot be rejected\n', L, vL);
else
fprintf('L %.2f vL %.2f null hypothesis can be rejected\n', L, vL);
end
fprintf('\n');
% Now do it lots of times to create empirical pdf's
% p(PHI), p(E) and p(L)
Nreal = 10000;
Pvec = zeros(Nreal,1);
Evec = zeros(Nreal,1);
Lvec = zeros(Nreal,1);
for ireal = (1:Nreal)
dobs = dtrue + random('Normal',0,sigma_d,N,1);
f = [ sqrtcovdi*dobs; sqrtcovhi*h ];
mest = FTFinv*(F'*f);
e = sqrtcovdi*(G*mest-dobs);
l = sqrtcovhi*(H*mest-h);
Evec(ireal) = e'*e;
Lvec(ireal) = l'*l;
Pvec(ireal) = Evec(ireal)+Lvec(ireal);
end
Evec = sort(Evec);
Lvec = sort(Lvec);
Pvec = sort(Pvec);
% empirical pdf's via histograms
Nbins = 201;
binmin = 0;
binmax = 2*vP;
Dbin = (binmax-binmin)/(Nbins-1);
bins = binmin + Dbin*[0:Nbins]';
Phist = hist( Pvec, bins );
Ehist = hist( Evec, bins );
Lhist = hist( Lvec, bins );
% plot histigram
figure();
clf;
% p(PHI)
subplot(3,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
top = 1.1*max(Phist);
axis( [binmin, binmax, 0, top ] );
plot( bins,Phist, 'k-', 'LineWidth', 3 );
plot( [vP, vP], [0, top/5], 'r-', 'LineWidth', 2 );
plot( [vP-2*sqrt(2*vP), vP-2*sqrt(2*vP)], [0, top/10], 'b-', 'LineWidth', 4 );
plot( [vP+2*sqrt(2*vP), vP+2*sqrt(2*vP)], [0, top/10], 'b-', 'LineWidth', 4 );
xlabel('P');
ylabel('counts');
% p(E)
subplot(3,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
top = 1.1*max(Ehist);
axis( [binmin, binmax, 0, top ] );
plot( bins,Ehist, 'k-', 'LineWidth', 3 );
plot( [vE, vE], [0, top/5], 'r-', 'LineWidth', 2 );
plot( [vE-2*sqrt(2*vE), vE-2*sqrt(2*vE)], [0, top/10], 'b-', 'LineWidth', 4 );
plot( [vE+2*sqrt(2*vE), vE+2*sqrt(2*vE)], [0, top/10], 'b-', 'LineWidth', 4 );
xlabel('E');
ylabel('counts');
% p(L)
subplot(3,1,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
top = 1.1*max(Lhist);
axis( [binmin, binmax, 0, top ] );
plot( bins,Lhist, 'k-', 'LineWidth', 3 );
plot( [vL, vL], [0, top/10], 'r-', 'LineWidth', 2 );
plot( [vL-2*sqrt(2*vL), vL-2*sqrt(2*vL)], [0, top/10], 'b-', 'LineWidth', 4 );
plot( [vL+2*sqrt(2*vL), vL+2*sqrt(2*vL)], [0, top/10], 'b-', 'LineWidth', 4 );
xlabel('L');
ylabel('counts');
fprintf('Caption: Empirical pdfs for (top) generalized error Phi,\n');
fprintf('(middle) prediction error E and (bottom) error in prior information L\n');
fprintf('Also shown are mean z(red) and 95 percent confidence bounds (blue)\n');
% empirical confidence limits based on the empirical pdf's
fprintf('Numerical 95%% confidence tests\n');
i=find(P<Pvec,1);
p = i/Nreal;
if ( p>0.025 && p<0.975 )
fprintf('P %.2f vP %.2f null hypothesis cannot be rejected\n', P, vP);
else
fprintf('P %.2f vP %.2f null hypothesis can be rejected\n', P, vP);
end
i=find(E<Evec,1);
p = i/Nreal;
if ( p>0.025 && p<0.975 )
fprintf('E %.2f vE %.2f null hypothesis cannot be rejected\n', E, vE);
else
fprintf('E %.2f vE %.2f null hypothesis can be rejected\n', E, vE);
end
i=find(L<Lvec,1);
p = i/Nreal;
if ( p>0.025 && p<0.975 )
fprintf('L %.2f vL %.2f null hypothesis cannot be rejected\n', L, vL);
else
fprintf('L %.2f vL %.2f null hypothesis can be rejected\n', L, vL);
end
% gdama06_18
% example F-test to assess difference in fit of two models
% the two models are linear and cubic fit of the same d(z)
% dataset
clear all;
fprintf('gdama06_18\n');
% auxially variable z
N = 11;
z = [0:N-1]'/(N-1);
dtrue = z - 0.5*z.*z;
% simulated data
sigmad=0.03;
dobs = dtrue + random('Normal', 0, sigmad, N, 1 );
figure();
clf;
% plot the observed data
subplot(2,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 1, -1, 1] );
plot(z,dobs,'ro','LineWidth',3);
% fit 1, straight line
M=2;
G=[ones(N,1), z];
mestA=(G'*G)\(G'*dobs);
dpreA = G*mestA;
% plot the predicted data for the linear fit
plot(z,dpreA,'b-','LineWidth',3);
title('linear fit');
xlabel('z');
ylabel('d(z)');
fprintf('Caption: Linear fit (blue curve) to data (red circles)\n');
% linear error
EA = (dobs-dpreA)'*(dobs-dpreA);
vA = N-M; % degrees of freedom
fprintf('linear error %f, degrees of freedom %d\n', EA, vA);
% plot the observed data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 1, -1, 1] );
plot(z,dobs,'ro','LineWidth',3);
% fit 2, cubic
M=4;
G=[ones(N,1), z, z.*z, z.^3];
mestB=(G'*G)\(G'*dobs);
dpreB = G*mestB;
% plot the predicted data for the cubic fit
plot(z,dpreB,'b-','LineWidth',3);
title('cubic fit');
xlabel('z');
ylabel('d(z)');
fprintf('Caption: Cubic fit (blue curve) to data (red circles)\n');
% error of cubic fit
EB = (dobs-dpreB)'*(dobs-dpreB);
vB = N-M; % degrees of freedom
fprintf('cubic error %f, degrees of freedom %d\n', EB, vB);
Fobs = (EA/vA) / (EB/vB); % F-value
fprintf('1/F %f F %f', 1/Fobs, Fobs);
% probability value associated with F-value
if( Fobs<1 )
Fobs=1/Fobs;
end
P = 1 - (fcdf(Fobs,vA,vB)-fcdf(1/Fobs,vA,vB));
Pleft = fcdf(1/Fobs,vA,vB);
Pright = 1-fcdf(Fobs,vA,vB);
fprintf('P(F<%f) = %f\n', 1/Fobs, Pleft);
fprintf('P(F>%f) = %f\n', Fobs, Pright);
fprintf('P(F<%f or F>%f) = %f\n', 1/Fobs, Fobs, P);
if( (Pleft+Pright)<0.05 )
fprintf('Null Hypothesis can be rejected to 95%% confidence\n');
else
fprintf('Null Hypothesis cannot be rejected to 95%% confidence\n');
end
% gdama06_19
% example of F test for problem with prior information
clear all;
% dense set of model parameters
M1=101; % must be odd
% densely sampled time t variable
Dt1=1;
t1 = Dt1*[0:M1-1]';
% sparse set of model parameters
M2=(M1+1)/2;
% sparsely sampled time t variable
Dt2=2*Dt1;
t2 = Dt2*[0:M2-1]';
% matrix D2S takes dense to sparse model parameters by decimation
% m2 = D2S*m1
D2S = zeros(M2, M1);
j=1;
for i=[1:M2]
D2S(i,j)=1.0;
j=j+2;
end
% matrix S2D takes sparse to dense model parameters by linear interpolation
% m1 = S2D*m2
S2D = zeros(M1, M2);
j=1;
for i=(1:2:M1)
S2D(i,j)=1;
j=j+1;
end
j=1;
for i=(2:2:M1-1)
S2D(i,j)=0.5;
S2D(i,j+1)=0.5;
j=j+1;
end
% make a densely sampled wiggly curve m(t) by starting out with random noise
% and bandpass filtering it around a narrow range of angular
% frquencies centered on w0. The bandpass filtering is
% accomplished by the gda_chebyshevfilt() function, which
% though not mentioned in the text (sorry about that) is pretty
% easy to use and has many other applications. Its arguments are
% output_timeseries = gda_chebyshevfilt( input_timeseries, Dt, f1, f2 )
% where Dt is the sampling of the timeseries (say in s) and
% (f1,f2) are the (low,high) size of the frequency band (say in Hz)
n = 10;
A = 2.0;
w0 = n*pi/t1(end);
Dw = w0/10;
mtrue1 = gda_chebyshevfilt( random('Normal',0,1,M1,1), Dt1, (w0-Dw)/(2*pi), (w0+Dw)/(2*pi) );
mtrue1 = A*mtrue1/max(abs(mtrue1));
% data kernel is just 3 independent observations of same point
N = 3*M1;
G1 = [eye(M1,M1); eye(M1,M1); eye(M1,M1)];
G2 = [S2D; S2D; S2D];
dtrue = G1*mtrue1;
% observed data is noisy
sigma_d = 0.1;
dobs = dtrue + random('Normal',0,sigma_d,N,1);
sqrtcovdi = eye(N,N)/sigma_d;
% prior information: second derivative close to zero
K1 = M1-2;
H1 = (1/Dt1^2)*toeplitz( [1; zeros(M1-3,1)], [1, -2, 1, zeros(1,M1-3) ] );
h1 = zeros(M1-2,1);
K2 = M2-2;
H2 = (1/Dt2^2)*toeplitz( [1; zeros(M2-3,1)], [1, -2, 1, zeros(1,M2-3) ] );
h2 = zeros(M2-2,1);
% prior covariance of model parameters
sigma_mdd = (w0^2)*A/sqrt(2);
sqrtcovhi1 = eye(K1,K1)/sigma_mdd;
sqrtcovhi2 = eye(K2,K2)/sigma_mdd;
% overall least squares equation and its solution
% dense problem
F1 = [ sqrtcovdi*G1; sqrtcovhi1*H1 ];
f1 = [ sqrtcovdi*dobs; sqrtcovhi1*h1 ];
FTFinv1 = inv(F1'*F1);
mest1 = FTFinv1*(F1'*f1);
% sparse problem
F2 = [ sqrtcovdi*G2; sqrtcovhi2*H2 ];
f2 = [ sqrtcovdi*dobs; sqrtcovhi2*h2 ];
FTFinv2 = inv(F2'*F2);
mest2 = FTFinv2*(F2'*f2);
% plot m
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, t1(end), -1.1*A, 1.1*A] );
plot( t1, dobs(1:M1), 'b.', 'LineWidth', 2 );
plot( t1, dobs(M1+1:2*M1), 'b.', 'LineWidth', 2 );
plot( t1, dobs(2*M1+1:3*M1), 'b.', 'LineWidth', 2 );
plot( t1, mtrue1, 'k-', 'LineWidth', 6 );
plot( t1, mest1, 'r-', 'LineWidth', 4 );
plot( t2, mest2, 'g:', 'LineWidth', 2 );
xlabel('t (s)');
ylabel('m');
fprintf('Caption: Bandlimited time series, true (black), estimated\n');
fprintf('sparse solution (red), estimated dense solution (green), data (blue dots)\n');
% error associated with dense problem
e1 = sqrtcovdi*(G1*mest1-dobs);
l1 = sqrtcovhi1*(H1*mest1-h1);
E1 = e1'*e1;
L1 = l1'*l1;
P1 = E1+L1;
vP1 = N+K1-M1;
vE1 = vP1*N/(N+K1);
vL1 = vP1*K1/(N+K1);
% error associated with sparse problem
e2 = sqrtcovdi*(G2*mest2-dobs);
l2 = sqrtcovhi2*(H2*mest2-h2);
E2 = e2'*e2;
L2 = l2'*l2;
P2 = E2+L2;
vP2 = N+K2-M2;
vE2 = vP2*N/(N+K2);
vL2 = vP2*K2/(N+K2);
% F ratio
vA=vP1;
vB=vP2;
PhiA = P1;
PhiB = P2;
Fobs = (PhiA/vA) / (PhiB/vB);
if( Fobs<1 )
Fobs=1/Fobs;
end
fprintf('1/F %f F %f\n', 1/Fobs, Fobs);
Pval = 1 - abs(fcdf(1/Fobs,vA,vB)-fcdf(Fobs,vA,vB));
Pleft = fcdf(1/Fobs,vA,vB);
Pright = 1-fcdf(Fobs,vA,vB);
fprintf('P(F<%f) = %f\n', 1/Fobs, Pleft);
fprintf('P(F>%f) = %f\n', Fobs, Pright);
fprintf('P(F<%f or F>%f) = %f\n', 1/Fobs, Fobs, Pval);
if( (Pleft+Pright)<0.05 )
fprintf('Null Hypothesis can be rejected to 95%% confidence\n');
else
fprintf('Null Hypothesis cannot be rejected to 95%% confidence\n');
end
% now do a whole lot of the same problems
% with different realizations of observational noise
Nreal = 10000;
FFvec = zeros(Nreal,1);
for ireal = (1:Nreal)
% add new noise to data
dobs = dtrue + random('Normal',0,sigma_d,N,1);
% dense problem
f1 = [ sqrtcovdi*dobs; sqrtcovhi1*h1 ];
mest1 = FTFinv1*(F1'*f1);
e1 = sqrtcovdi*(G1*mest1-dobs);
l1 = sqrtcovhi1*(H1*mest1-h1);
PP1 = (e1'*e1+l1'*l1);
% add new noise to data
dobs = dtrue + random('Normal',0,sigma_d,N,1);
% sparse problem
f2 = [ sqrtcovdi*dobs; sqrtcovhi2*h2 ];
mest2 = FTFinv2*(F2'*f2);
e2 = sqrtcovdi*(G2*mest2-dobs);
l2 = sqrtcovhi2*(H2*mest2-h2);
PP2 = (e2'*e2+l2'*l2);
% tabulate F ratio
FFvec(ireal) = (PP1/vP1) / (PP2/vP2);
end
FFvec = sort(FFvec);
FFmean = mean(FFvec);
% histogram and empirical p.d.f.
Nbins = 101;
binmin = 0.5;
binmax = 1.5;
Dbin = (binmax-binmin)/(Nbins-1);
bins = binmin + Dbin*[0:Nbins]';
FFhist = hist( FFvec, bins );
FFhist = FFhist/(Dbin*sum(FFhist));
% plot pdf's
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
top = 1.1*max(FFhist);
axis( [binmin, binmax, 0, top ] );
plot( bins,FFhist, 'k-', 'LineWidth', 3 );
plot( bins,fpdf(bins,vP1,vP2), 'r-', 'LineWidth', 3 );
plot( [Fobs,Fobs]', [0, top/5]', 'k-', 'LineWidth', 4 );
plot( [FFmean,FFmean]', [0, top/5]', 'g-', 'LineWidth', 4 );
plot( [finv(0.5,vP1,vP2),finv(0.5,vP1,vP2)]', [0, top/5]', 'r-', 'LineWidth', 4 );
plot( [finv(0.025,vP1,vP2),finv(0.025,vP1,vP2)]', [0, top/8]', 'b-', 'LineWidth', 4 );
plot( [finv(0.975,vP1,vP2),finv(0.975,vP1,vP2)]', [0, top/8]', 'b-', 'LineWidth', 4 );
xlabel('F');
ylabel('p(F)')
fprintf('Caption: Analytic pdf p(F) (red curve), empirical pdf (black curve)\n');
fprintf('observed F (black bar), empirical mean (green bar), true mean (red bar)\n');
fprintf('95 percent confidence (blue bars)\n');