Untitled

% gdama12_01

% Monte Carlo search

% applied to exemplary curve-fitting problem

% tunable parameters

% variance of data

clear all;

fprintf('gdama12_01\n');

gdama12_01

sd = 0.3;

% x-axis

N=40;

xmin=0;

xmax=1.0;

Dx=(xmax-xmin)/(N-1);

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

% two model parameters

M=2;

% true model parameters

mt = [1.21, 1.54]';

% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);

w0=20.0;

dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2); % true data

% observed data

sd2 = sd^2;

dobs = dtrue + random('Normal',0.0,sd,N,1);

% plot data

figure();

clf;

hold on;

axis( [0, xmax, 0, 4 ] );

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

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

xlabel('x"');

ylabel('d');

fprintf('Caption: True data (black curve) and observed data (red circles)\n');

Caption: True data (black curve) and observed data (red circles)

% Define 2D grid (for plotting purposes only)

L = 101;

Dm = 0.02;

m1min=0.0;

m2min=0.0;

m1a = m1min+Dm*[0:L-1]';

m2a = m2min+Dm*[0:L-1]';

m1max = m1a(L);

m2max = m2a(L);

% tabulate total error E on grid (for plotting purposed only)

% Note m1 varies along rows of E and m2 along columns

E = zeros(L,L);

for j = (1:L)

for k =(1:L)

dpre = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);

e = dobs-dpre;

E(j,k) =e'*e/sd2;

end

Emax = max(E(:));

% Monte Carlo search

% save best solutions, so far

Niter = 10000;

mbest = zeros(M,Niter);

Ebest = zeros(Niter,1);

iterbest = zeros(Niter,1);

Nsaved = 0;

% starting solution

m0 = [0.1;0.1]; % start in upepr-left hand corner of plot

% predicted data

dpre0 = sin(w0*m0(1)*x) + m0(1)*m0(2);

% prediction error

e0 = dobs-dpre0;

E0 = e0'*e0/sd2;

% save this solution first solution, the best so far

Nsaved=Nsaved+1;

mbest(1:M,Nsaved) = m0;

Ebest(Nsaved) = E0;

iterbest(Nsaved) = 0;

for myiter = (2:Niter)

% randomly select trial solution

m1trial = random('Uniform',m1min,m1max);

m2trial = random('Uniform',m2min,m2max);

mtrial = [m1trial; m2trial];

% predict data

dtrial = sin(w0*m1trial*x) + m1trial*m2trial;

% prediction error

etrial = dobs-dtrial;

Etrial = etrial'*etrial/sd2;

% accept trial solution if error decreases

if( Etrial<E0 )

Nsaved=Nsaved+1;

mbest(1:M,Nsaved) = mtrial;

Ebest(Nsaved) = Etrial;

iterbest(Nsaved) = myiter;

% reset best-solution so far to this one

E0=Etrial;

m0 = mtrial;

end

% throw away unused part of arrays

mbest = mbest(1:2,1:Nsaved);

Ebest = Ebest(1:Nsaved);

iterbest = iterbest(1:Nsaved);

% plot error surface

figure();

clf;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

hold on;

axis( [m2min, m2max, m1min, m1max] );

axis ij;

imagesc( [m2min, m2max], [m1min, m1max], E);

plot( mbest(2,:), mbest(1,:), 'w*', 'LineWidth',2);

plot( mbest(2,:), mbest(1,:), 'w-', 'LineWidth',2);

plot( mbest(2,end), mbest(1,end), 'ro', 'LineWidth',2);

plot( mt(2), mt(1), 'g*', 'LineWidth', 2);

colorbar;

xlabel('m2');

ylabel('m1');

fprintf('Caption: Error surface, with true solution (green) and progression of solution');

Caption: Error surface, with true solution (green) and progression of solution

fprintf('(white circles and lines) and final solution\n');

(white circles and lines) and final solution

fprintf('(red circle)\n');

(red circle)

% plot progress

figure();

clf;

hold on;

K = max(iterbest);

axis( [0, K, 0, 4 ] );

plot( iterbest(:), mbest(1,:),'k-','LineWidth',1);

plot( [0,K]',[mt(1),mt(1)]','r:','LineWidth',1);

plot( iterbest(:), mbest(1,:),'k*','LineWidth',3);

xlabel('iteration, i');

ylabel('m1');

fprintf('Caption: Evolution of model parameter 1 with iteration number (black),\n');

Caption: Evolution of model parameter 1 with iteration number (black),

fprintf('with true model parameter depicted by red line\n');

with true model parameter depicted by red line

figure();

clf;

hold on;

axis( [0, K, 0, 4 ] );

plot( iterbest(:), mbest(2,:),'k-','LineWidth',1);

plot( [1,K]',[mt(2),mt(2)]','r:','LineWidth',1);

plot( iterbest(:), mbest(2,:),'k*','LineWidth',3);

xlabel("iteration, i");

ylabel("m2");

fprintf('Caption: Evolution of model parameter 2 with iteration number (black),\n');

Caption: Evolution of model parameter 2 with iteration number (black),

fprintf('with true model parameter depicted by red line\n');

with true model parameter depicted by red line

figure();

clf;

hold on;

axis( [0, K, min(log10(Ebest)), max(log10(Ebest)) ] );

plot( iterbest(:), log10(Ebest(:)), 'k-','LineWidth',1);

plot( iterbest(:), log10(Ebest(:)), 'k*','LineWidth',3);

xlabel('iteration, i');

ylabel('log10(E)');

fprintf('Caption: Evolution of error with iteration number\n');

Caption: Evolution of error with iteration number

fprintf('%d iterations, %d accepted\n', Niter, Nsaved );

10000 iterations, 10 accepted

% gdama12_02

% Simulated Annealing search

% applied to exemplary curve-fitting problem

clear all;

fprintf('gdama12_02');

gdama12_02

% tunable parameters

% variance of data

sd = 0.3;

% neighborhood parameter

sq = 0.3;

% x axis

N=40;

xmin=0.0;

xmax=1.0;

Dx=(xmax-xmin)/(N-1);

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

% two model parameters

M=2;

% true model parameters

mt = [1.21; 1.54];

% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);

w0=20.0;

dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2); % true data

% observed data

sd2 = sd^2;

dobs = dtrue + random('Normal',0.0,sd,N,1);

% plot data

figure();

clf;

hold on;

subplot(1,1,1);

axis( [0, xmax, 0, 4 ] );

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

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

xlabel('x');

ylabel('d');

fprintf('Caption: True data (black curve) and observed data (red circles)\n');

Caption: True data (black curve) and observed data (red circles)

% Define 2D grid (for plotting purposes only)

L = 101;

Dm = 0.02;

m1min=0.0;

m2min=0.0;

m1a = m1min+Dm*[0:L-1]';

m2a = m2min+Dm*[0:L-1]';

m1max = m1a(L);

m2max = m2a(L);

% tabulate total error psi on grid.

% Note m1 varies along rows of E and m2 along columns

E = zeros(L,L);

for j = (1:L)

for k = (1:L)

dpre = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);

e = dobs-dpre;

E(j,k) = e'*e/sd2;

end

Emax = max(E(:));

% Simulted annealing

% save all solutions and their error

Niter = 4000;

msave = zeros(M,Niter);

Esave = zeros(Niter,1);

n = [1:Niter]';

% temperature

Tstart = 100.0;

Tend = 0.001;

T = exp( log(Tstart)+(log(Tend)-log(Tstart))*[0:Niter-1]'/(Niter-1) );

% starting solution

m0 = [0.2; 0.2]; % start in upepr-left hand corner of plot

% predicted data

d0 = sin(w0*m0(1)*x) + m0(1)*m0(2);

% prediction error

e0 = dobs-d0;

E0 = e0'*e0;

% likelihood

logpm0 = -E0/(T(1)*sd2);

% save this solution first solution, the best so far

msave(1:M,1) = m0;

Esave(1) = E0;

Naccept = 0;

for myiter = (2:Niter)

% randomly select trial solution

mp = random('Normal',m0,sq,M,1);

% predict data

dp = sin(w0*mp(1)*x) + mp(1)*mp(2);

% prediction error

ep = dobs-dp;

Ep = ep'*ep;

logpmp = -0.5*Ep/(T(myiter)*sd2);

logalpha = logpmp - logpm0;

% standard metropolis acceptance test

if( logalpha >= 0.0 )

m0 = mp;

E0 = Ep;

logpm0 = logpmp;

Naccept=Naccept+1;

else

alpha = exp(logalpha);

r = random('Uniform',0.0,1.0);

if( alpha > r )

m0 = mp;

E0 = Ep;

logpm0 = logpmp;

Naccept=Naccept+1;

else

mp = m0;

Ep = E0;

logpmp = logpm0;

end

msave(1:M,myiter) = mp;

Esave(myiter) = Ep;

end

% plot error surface

figure();

clf;

hold on;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

axis( [m2min, m2max, m1min, m1max] );

axis ij;

imagesc( [m2min, m2max], [m1min, m1max], E);

colorbar;

xlabel('m2');

ylabel('m1');

plot( msave(2,:), msave(1,:), 'w*', 'LineWidth', 2);

plot( msave(2,:), msave(1,:), 'w-', 'LineWidth', 2);

plot( mt(2), mt(1), 'g*', 'LineWidth', 2);

plot( msave(1,Niter), msave(1,Niter), 'r*', 'LineWidth', 2);

fprintf('Caption: Error surface, with true solution (green), starting solution (red star),\n');

Caption: Error surface, with true solution (green), starting solution (red star),

fprintf('progression of solution (white circles and lines) and final solution\n');

progression of solution (white circles and lines) and final solution

fprintf('(red circle)\n');

(red circle)

% plot progress

figure();

clf;

hold on;

axis( [0, Niter, 0, 4 ] );

plot( n, msave(1,:)','k-', 'LineWidth', 1);

plot( [1,Niter]',[mt(1),mt(1)]','r:','LineWidth', 1);

xlabel('iteration, i');

ylabel('m1');

fprintf('Caption: Evolution of model parameter 1 with iteration number (black),\n');

Caption: Evolution of model parameter 1 with iteration number (black),

fprintf('with true model parameter depicted by red line\n');

with true model parameter depicted by red line

figure();

clf;

hold on;

axis( [0, Niter, 0, 4 ] );

plot( n, msave(2,:)','k-','LineWidth', 1);

plot( [1,Niter]',[mt(2),mt(2)]','r:','LineWidth', 1);

xlabel('iteration, i');

ylabel('m2');

fprintf('Caption: Evolution of model parameter 2 with iteration number (black),\n');

Caption: Evolution of model parameter 2 with iteration number (black),

fprintf('with true model parameter depicted by red line\n');

with true model parameter depicted by red line

figure();

clf;

hold on;

axis( [0, Niter, min(log10(Esave)), max(log10(Esave)) ] );

plot( n, log10(Esave), 'k-', 'LineWidth', 1);

xlabel('iteration, i');

ylabel('log10(E)');

fprintf('Caption: Evolution of error with iteration number\n');

Caption: Evolution of error with iteration number

figure();

clf;

hold on;

axis( [0, Niter, min(log(T)), max(log(T)) ] );

plot( n, log10(T),'k-','LineWidth', 1);

xlabel("iteration, i");

ylabel("log10(T)");

fprintf('Caption: Decrease in temperture with iteration number\n');

Caption: Decrease in temperture with iteration number

fprintf('%d iterations, %d accepts\n', Niter,Naccept );

4000 iterations, 156 accepts

% gdama12_03

% Using the Metropolis-Hastings algorithm to sample the likelihood

% associated with the exemplary curve-fitting problem

clear all;

fprintf('gdama12_03\n');

gdama12_03

% tunable parameters

% variance of data

sd = 0.3;

% variance of prior information

sm = 0.5;

% neighborhood parameter

sn = 0.05;

% x-axis

N=40;

xmin=0;

xmax=1.0;

Dx=(xmax-xmin)/(N-1);

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

% two model parameters

M=2;

% true model parameters

mt = [1.21, 1.54]';

% prior model parameters and their variance

mA = [1.43, 1.50]';

sm2 = sm^2;

% d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);

w0=20.0;

dtrue = sin(w0*mt(1)*x) + mt(1)*mt(2); % true data

% observed data

sd2 = sd^2;

dobs = dtrue + random('Normal',0.0,sd,N,1);

% plot data

figure();

clf;

hold on;

axis( [0, xmax, 0, 4 ] );

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

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

xlabel('x');

ylabel('d');

fprintf('Caption: True data (black) and observed data (red)\n');

Caption: True data (black) and observed data (red)

% Define 2D grid

L = 101;

Dm = 0.02;

m1min=0.0;

m2min=0.0;

m1a = m1min+Dm*[0:L-1]';

m2a = m2min+Dm*[0:L-1]';

m1max = m1a(L);

m2max = m2a(L);

% tabulate total error psi on grid (for plotting purposes only)

% Note m1 varies along rows of psi and m2 along columns

Psi = zeros(L,L);

for j = (1:L)

for k = (1:L)

dpre = sin(w0*m1a(j)*x) + m1a(j)*m2a(k);

e = dobs-dpre;

l = mA - [m1a(j), m2a(k)]';

Psi(j,k) = e'*e/sd2 + l'*l/sm2;

end

Psimax = max(Psi(:));

% Metropolis Hastings

Nburnt=100000; % burn in

Niter=1000000+Nburnt;

m = zeros(M,Niter);

m1c = (m1min + m1max)/2.0;

m2c = (m2min + m2max)/2.0;

mold = [m1c, m2c]';

m(1:M,1) = mold;

dold = sin(w0*mold(1)*x) + mold(1)*mold(2);

e = dobs-dold;

Eold = e'*e;

l = mA-mold;

Lold = l'*l;

logpold = -0.5*Eold/sd2-0.5*Lold/sm2;

for myiter = (1:Niter)

% perturbation

dm = random('Normal',0.0,sn,M,1);

mnew = mold + dm;

% prior error

l = mA-mnew;

Lnew = l'*l;

% prediction

dnew = sin(w0*mnew(1)*x) + mnew(1)*mnew(2);

e = dobs-dnew;

Enew = e'*e;

logpnew = -0.5*Enew/sd2-0.5*Lnew/sm2;

% acceptance test

if( (logpnew-logpold) >= 0.0)

m(1:M,myiter) = mnew;

mold = mnew;

Eold = Enew;

Lold = Lnew;

logpold = logpnew;

else

alpha = exp( logpnew - logpold );

r = random('Uniform',0.0,1.0);

if( alpha > r )

m(1:M,myiter) = mnew;

mold = mnew;

Eold = Enew;

Lold = Lnew;

logpold = logpnew;

else

m(1:M,myiter) = mold;

end

% keep only values past burn in

m = m(1:M,Nburnt+1:Niter);

[i, Niter] = size(m);

% estimate is the mean

mest = mean(m,2);

% sample covariance

Cm = cov(m');

% 95% confidence of the mean

Dm95 = 2.0*sqrt(diag(Cm)/Niter);

% plot error surface

figure();

clf;

hold on;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

axis( [m2min, m2max, m1min, m1max] );

axis ij;

imagesc( [m2min, m2max], [m1min, m1max], Psi);

colorbar;

xlabel('m2');

ylabel('m1');

plot( mest(2), mest(1), 'go', 'LineWidth', 3);

plot( [mest(2)-Dm95(2), mest(2)+Dm95(2)]', [mest(1)-Dm95(1), mest(1)+Dm95(1)]', 'g-', 'LineWidth', 1);

plot( [mest(2),mest(2)]',[mest(1)-Dm95(1),mest(1)+Dm95(1)]', 'g-', 'LineWidth', 1);

plot( mA(2), mA(1), 'y*', 'LineWidth', 3);

fprintf('Caption: Generalized error surface with prior model (yellow circle)\n');

Caption: Generalized error surface with prior model (yellow circle)

fprintf('estimated model (green circle) and its 95 percent confident\n');

estimated model (green circle) and its 95 percent confident

fprintf('bounds (green lines) superimposed. The confidence bounds may\n');

bounds (green lines) superimposed.  The confidence bounds may

fprintf('be negligibly small\n');

be negligibly small

% histogram of model

H=zeros(L,L);

for myiter = (1:Niter)

mym1 = m(1,myiter);

mym2 = m(2,myiter);

i = floor((mym1-m1min)/Dm)+1; % convert m1 to integer array index

j = floor((mym2-m2min)/Dm)+1; % convert m2 to integer array index

if( (i<1) || (i>L) || (j<1) || (j>L) )

continue;

end

H(i,j)=H(i,j)+1;

end

% plot histogram of model

figure();

clf;

hold on;

set(gca,'LineWidth',3);

set(gca,'FontSize',14);

colormap('jet');

axis( [m2min, m2max, m1min, m1max] );

axis ij;

imagesc( [m2min, m2max], [m1min, m1max], H);

colorbar;

xlabel('m2');

ylabel('m1');

fprintf('Caption: Histogram of model ensemble\n');

Caption: Histogram of model ensemble

% histogram of predicted data

NHx=512;

NHv=128;

xH=xmax*[1:NHx]'/(NHx-1);

vmin=0.0;

vmax=4.0;

Dv = (vmax-vmin)/(NHv-1);

H=zeros(NHv,NHx);

n = [1:NHx]';

for myiter = (1:Niter)

dH = sin(w0*m(1,myiter)*xH) + m(1,myiter)*m(2,myiter);

j = floor((dH-vmin)/Dv)+1; % convert data to integer array index

k = find((j>=1)&(j<=NHv)); % exclude indices that are out of bound

jg = j(k); % data value indices

ig = n(k); % corresponding x position indices

[Nk,i]=size(k);

for i=(1:Nk)

H(jg(i),ig(i))=H(jg(i),ig(i))+1; % accumulate histogram

end

% plot data

figure();

clf;

hold on;

colormap('jet');

axis( [0, xmax, vmin, vmax ] );

Hmin=0.0;

Hmax = max(H(:));

caxis([Hmin,Hmax]);

imagesc([0, xmax], [vmin, vmax], H);

plot(x,dobs,'wo','LineWidth', 3);

xlabel('x');

ylabel('d');

colorbar();

fprintf('Caption: Histogram of predicted data with observed\n');

Caption: Histogram of predicted data with observed

fprintf('data (white circles) superimposed');

data (white circles) superimposed

% gdama12_04

% Using the Metropolis-Hastings algorithm to sample the likelihood

% associated with the exemplary Laplace transform problem

clf;

fprintf('gdama12_04\n');

gdama12_04

% z-axis

M=11;

zmin=0;

zmax=1.0;

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

z = Dz*[0:M-1]';

% true model

mtrue = ones(M,1)+sin(2*pi*z/zmax);

% data kernel, exponentially decaying

N=M;

G = zeros(N,M);

c = 6.0*z; % decay rates

for i = (1:M)

G(i,1:M) = exp(-c(i)*z');

end

% prior model parameters and their variance

mA = ones(M,1);

sm = 0.1;

sm2 = sm^2;

% d = G m

dtrue = G*mtrue;

% observed data

sd=0.01;

sd2 = sd^2;

dobs = dtrue + random('Normal',0.0,sd,N,1);

gda_draw( G, 'caption G', mtrue,'caption m', '=', dobs, 'caption dobs');

fprintf('Caption: Graphical depiction of data equation\n');

Caption: Graphical depiction of data equation

% plot data

figure();

clf;

hold on

axis( [0, zmax, 0, max(dobs) ] );

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

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

xlabel('z');

ylabel('d');

fprintf('Caption: True data (black) and observed data (red)\n');

Caption: True data (black) and observed data (red)

% neighborhood parameter

sn = 0.2;

sn = 0.05;

% Metropolis Hastings

Nburnt=1000; % burn in

Niter=1000000+Nburnt;

m = zeros(M,Niter);

mold = random('uniform',0.8,1.0,M,1);

m(1:M,1) = mold;

dold = G*mold;

e = dobs-dold;

Eold = e'*e;

l = mA-mold;

Lold = l'*l;

logpold = -0.5*Eold/sd2-0.5*Lold/sm2;

for myiter = (2:Niter)

% perturbation

Dm = random('Normal',0.0,sn,M,1);

mnew = mold + Dm;

% prior error

l = mA-mnew;

Lnew = l'*l;

% prediction

dnew = G*mnew;

e = dobs-dnew;

Enew = e'*e;

logpnew = -0.5*Enew/sd2-0.5*Lnew/sm2;

% acceptance test

Dlogp = logpnew - logpold;

if( Dlogp >= 0.0)

m(1:M,myiter) = mnew;

mold = mnew;

Eold = Enew;

Lold = Lnew;

logpold = logpnew;

else

alpha = exp( Dlogp );

r = random('Uniform',0.0,1.0);

if( alpha > r )

m(1:M,myiter) = mnew;

mold = mnew;

Eold = Enew;

Lold = Lnew;

logpold = logpnew;

else

m(1:M,myiter) = mold;

end

% keep only values past burn in

m = m(1:M,Nburnt+1:Niter);

[i, Niter] = size(m);

% estimate is the mean, variance

% sample mean

mest = mean(m,2);

% sample covariance

Cm = cov(m');

% two times std error of mean

Dm95 = 2.0*sqrt(diag(Cm)/Niter);

% plot model

figure();

clf;

hold on;

axis( [0, zmax, 0, 3 ] );

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

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

plot( z, mest-Dm95, 'g-', 'LineWidth',1);

plot( z, mest+Dm95, 'g-', 'LineWidth',1);

xlabel('z');

ylabel('m');

fprintf('Caption: True model (black) and mean (estimated) model (red)\n');

Caption: True model (black) and mean (estimated) model (red)

fprintf('and its 95 percent confidence bounds (green) (confidence\n');

and its 95 percent confidence bounds (green) (confidence

fprintf('bounds may be very close to one another)');

bounds may be very close to one another)

% histogram

Nbins = 101;

binmin=0.0;

binmax=3.0;

Dbin = (binmax-binmin)/(Nbins-1);

bins = binmin+Dbin*[0:Nbins-1]';

H = zeros(Nbins,M);

for iter = (1:Niter)

for i = (1:M)

m0 = m(i,iter);

k = floor((m0-binmin)/Dbin)+1;

if( (k<1) || (k>Nbins) )

continue;

end

H(k,i)=H(k,i)+1.0;

end

% plot error surface

figure();

clf;

hold on

colormap('jet');

axis( [0, zmax, binmin, binmax] );

Hroot = sqrt(H);

Hmax = max(Hroot(:));

caxis([0.0,Hmax]);

imagesc([0, zmax], [binmin, binmax], Hroot);

plot( z, mtrue, 'y-', 'LineWidth',3);

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

colorbar();

xlabel('z');

ylabel('m');

fprintf('Caption: Histogram of model, with mean (estimated) model (white)\n');

Caption: Histogram of model, with mean (estimated) model (white)

fprintf('and true model (yellow) superimposed\n');

and true model (yellow) superimposed

count=sum( (m(3,:)>m(2,:)) &(m(3,:)>m(4,:)) );

fprintf('%.6f percent have m2>m1 and m2>m3\n', 100.0*count/Niter );

61.482500 percent have m2>m1 and m2>m3

% gdama12_05

% Using the extended Metropolis-Hastibgs algorith to sample the a

% trans-dimensional Normal pdf with two data types:

% p(m) = [1d Normal with probability beta] + [2d Normal with probability (1-beta)]

clear all;

fprintf('gdama12_05\n');

gdama12_05

beta = 0.2; % probaility of dim 1

% dim-1 pdf is normal with this variance and normalization facto

sm1 = 0.5;

sm12 = sm1^2;

N1 = beta/sqrt(2.0*pi*sm12);

logN1 = log(N1);

% dim-2 pdf is normal with this variance and normalization factor

sm2 = 1.0;

sm22 = sm2^2;

N2 = (1.0-beta)/(2.0*pi*sm22);

logN2 = log(N2);

% Metropolis Hastings

Nburnt=1000; % burn in

Niter=1000000+Nburnt; % length of Markov chain

doswitch=500; % switch dimensions

% stores links in chain

m = zeros(3,Niter); % m[3,.] 1 or 2 signifies dimension

dimold = 2; % starting dimenion

dimnew = dimold;

% assign first link in chain

if( dimold==1 )

mold = random('Uniform',-1.0,1.0);

logpold = logN1-0.5*(mold^2)/sm12;

m(1,1) = mold;

m(1,2) = 0.0;

m(1,3) = dimold;

elseif( dimold==2 )

mold = random('Uniform',-1.0,1.0,2,1);

logpold = logN2-0.5*(mold(1)^2+mold(2)^2)/sm22;

m(1,1) = mold(1);

m(1,2) = mold(2);

m(1,3) = dimold;

end

% random verctors u are normal with these variance and normlizations

su1 = 0.2;

su12 = su1^2;

Nu1 = 1/(sqrt(2.0*pi)*su1);

logNu1 = log(Nu1);

su2 = 0.2;

su22 = su2^2;

Nu2 = 1/(sqrt(2.0*pi)*su2);

logNu2 = log(Nu2);

su3 = 1.0;

su32 = su3^2;

Nu3 = 1/(sqrt(2.0*pi)*su3);

logNu3 = log(Nu3);

% Transdimensional Metropolis Hastings

for myiter = (2:Niter)

% switch dimensions every doswitch iterations

if( mod(myiter,doswitch)==0 )

if( dimnew == 1 )

dimnew = 2;

elseif ( dimnew == 2 )

dimnew = 1;

end

% three random vectors u1, u2, u3 with transformation

% for the dim-1 to dim-2 transformation

% m1' = [ 1 1 0 0 ] [m1] (m1' = m1 + u1)

% m2' = [ 0 0 0 1 ] [u1] (m2' = u3)

% u1' = [ 0 1 0 0 ] [u2] (u1' = u1)

% u2' = [ 0 0 1 0 ] [u3] (u2' = u2)

% which has a Jacobian of unity

u1 = random('Normal',0.0,su1); % g1

u2 = random('Normal',0.0,su2); % g2

u3 = random('Normal',0.0,su3); % g3

if( (dimold==1) && (dimnew==1) )

mnew = mold + u1;

logpnew = logN1-0.5*(mnew^2)/sm12;

logalpha=logpnew-logpold;

elseif ( (dimold==1) && (dimnew==2) )

% pnew = p12' g1' g2'

% pold = p1 g1 g2 g3

% but since g1=g1' and g2=g2'

% so pnew/pold = p12' / (p1 g3)

logg3 = logNu3-0.5*(u3^2)/su32;

mnew = zeros(2,1);

mnew(1)= mold(1) + u1;

mnew(2)= u3;

logpnew = logN2-0.5*(mnew(1)^2+mnew(2)^2)/sm22;

logalpha=logpnew-(logpold+logg3);

elseif( (dimold==2) && (dimnew==2) )

mnew = zeros(2,1);

mnew(1)= mold(1) + u1;

mnew(2)= mold(2) + u2;

logpnew = logN2-0.5*(mnew(1)^2+mnew(2)^2)/sm22;

logalpha=logpnew-logpold;

elseif( (dimold==2) && (dimnew==1) )

mnew = mold(1) - u1;

logpnew = logN1-0.5*(mnew^2)/sm12;

up3 = mold(2);

loggp3 = logNu3-0.5*(up3^2)/su32;

logalpha=(loggp3+logpnew)-logpold;

end

if( logalpha > 0.0)

dimold = dimnew;

mold = mnew;

logpold = logpnew;

else

alpha=exp(logalpha);

r = random('Uniform',0.0,1.0);

if( alpha > r )

dimold = dimnew;

mold = mnew;

logpold = logpnew;

else

dimnew = dimold;

mnew = mold;

logpnew = logpold;

end

if( dimnew == 1 )

m(1,myiter)=mnew;

m(2,myiter)=0.0;

m(3,myiter)=dimnew;

elseif( dimnew == 2 )

m(1,myiter)=mnew(1);

m(2,myiter)=mnew(2);

m(3,myiter)=dimnew;

end

% keep only values past burn in

m = m(1:3,Nburnt+1:Niter);

[i, Niter] = size(m);

% histogram

Nbins = 41;

binmin = -3.0;

binmax = 3.0;

Dbin = (binmax-binmin)/(Nbins-1);

hbin = binmin+Dbin*[0:Nbins-1]';

H1 = zeros(Nbins,1);

H1count=0;

for myiter = (1:Niter)

if( m(3,myiter)==1 )

m1 = m(1,myiter);

k = floor((m1-binmin)/Dbin)+1;

if( (k<1) || (k>Nbins) )

continue;

end

H1(k) = H1(k)+1;

H1count = H1count+1;

end

H2 = zeros(Nbins,Nbins);

H2count=0;

for myiter = (1:Niter)

if( m(3,myiter)==2 )

m1 = m(1,myiter);

m2 = m(2,myiter);

k1 = floor((m1-binmin)/Dbin)+1;

k2 = floor((m2-binmin)/Dbin)+1;

if( (k1<1) || (k1>Nbins) || (k2<1) || (k2>Nbins) )

continue;

end

H2(k1,k2) = H2(k1,k2)+1;

H2count = H2count+1;

end

gda_draw(H1,H2);

fprintf('Caption: Estimated transdimensionl pdf, with\n');

Caption: Estimated transdimensionl pdf, with

fprintf('1-d compoent pdf (left) and 2-d (right)\n');

1-d compoent pdf (left) and 2-d (right)

fprintf("beta: true: %.4f est: %.4f\n", beta, H1count/(H1count+H2count) );

beta: true: 0.2000  est: 0.2028

% true H1 and H2

H1t = zeros(Nbins,1);

for i = (1:Nbins)

H1t(i) = N1*exp(-0.5*(hbin(i)^2)/sm12);

end

H2t = zeros(Nbins,Nbins);

for i = (1:Nbins)

for j = (1:Nbins)

H2t(i,j) = N2*exp(-0.5*(hbin(i)^2+hbin(j)^2)/sm22);

end

gda_draw(H1t,H2t);

fprintf('Caption: True transdimensionl pdf, with\n');

Caption: True transdimensionl pdf, with

fprintf('1-d compoent pdf (left) and 2-d (right)\n');

1-d compoent pdf (left) and 2-d (right)

H1n = (H1/(H1count+H2count))/(Dbin);

H2n = (H2/(H1count+H2count))/(Dbin^2);

figure();

clf;

hold on;

plot(hbin+Dbin/2,H1n,'k-'); % nis-registered by half a bin

plot(hbin,H1t,'r-');

xlabel('m');

ylabel('p1(m)');

fprintf('Caption: 1D component of transdimensionl pdf, true (red), empirical (black)\n');

Caption: 1D component of transdimensionl pdf, true (red), empirical (black)

% gdama12_06

% Using the extended Metropolis-Hastibgs algorith to sample the a

% trans-dimensional pdf that simultageously fits data to a sinusoid

% or a bell corve

% p(m|d) = [sinusoid, 1 model parameter] + [bell curve, 2 model parameters]

clear all;

fprintf('gdama12_06\n');

gdama12_06

fprintf(' (with true dimension of 1)\n');

 (with true dimension of 1)

% tunable parameter:

% true dimension, can be either 1 or 2

dimtrue = 1;

% x-axis

Nx=41;

xmin = 0.0;

xmax = 1.0;

xmid = (xmax-xmin)/2.0;

Dx = (xmax-xmin)/(Nx-1);

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

% model 1, sine wave of amplitude m1

if( dimtrue == 1)

mtrue = 1.0;

dtrue = mtrue*sin(pi*(x-xmin)/(xmax-xmin));

smnt2 = 0.2^2;

dnottrue = mtrue*exp(-0.5*((x-xmid).^2)/smnt2);

elseif( dimtrue == 2) % model 1, gaussian amplitude mp1 and variance mp2

mptrue = [1.0, 0.2^2 ]';

dtrue = mptrue(1)*exp(-0.5*((x-xmid).^2)/mptrue(2));

dnottrue = mptrue(1)*sin(pi*(x-xmin)/(xmax-xmin));

end

% prior variance of data

sd = 0.2;

sd2 = sd^2;

dobs = dtrue + random('Normal',0.0,sd,Nx,1);

% plot data

figure();

clf;

hold on;

axis( [xmin, xmax, 0, 2 ] );

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

plot(x,dnottrue,'g-','LineWidth', 3);

plot(x,dobs,'r.','LineWidth', 2);

xlabel('z');

ylabel('d');

fprintf('Caption: true data curve (black) and alternative data curve (green)\n');

Caption: true data curve (black) and alternative data curve (green)

fprintf('witn observed data (red) superimposed\n');

witn observed data (red) superimposed

% dim=1

% p(m|d) = p(d|m) pA(m) / p(d)

% dim=2

% pp(mp1,mp2|d) = pp(d|mp1,mp2) ppA(m1,m2) / p(d)

% note that p(d) cancels for all ratios, and can be ignored

% note that normalization p(d|m) and pp(d|mp1,mp2) is the same

% and so cancels for all rations and can be ignored

% pA(m) assumed to be normal

% pA(m) = NA exp(-0.5*(m-mA)**2/sA2) with NA = 1/sqrt(2*pi*sA2)

% ppA(mp1,mp2) asumed to be normal

% ppA(mp1,mp2) = NpA exp(-0.5*((mp1-mpA1)**2+(mp2-mpA2)**2)/spA2)

% with NpA = 1/(2*pi*spA2)

% prior information

% dim=1

mA = 0.9;

sA = 1.0;

sA2 = sA^2;

NA = 1/sqrt(2*pi*sA2);

logNA = log(NA);

% dim=2

mpA = [0.9,0.4^2]';

spA = 1.0;

spA2 = spA^2;

NpA = 1/(2*pi*spA2);

logNpA = log(NpA);

% Metropolis Hastings

Nburnt=10000; % burn in

Niter=1000000+Nburnt; % links in Markov chain

doswitch=50; % switch dimensions

% stores links in chain

m = zeros(3,Niter); % m[3,.] 1 or 2 signifies dimension

dimold = 2; % starting dimension

dimnew = dimold;

% assign first link in chain

if( dimold==1 )

mold = random('Uniform',0.5,1.5);

dpre = mold*sin(pi*(x-xmin)/(xmax-xmin));

e = dobs-dpre;

E = e'*e;

l = mold - mA;

L = l^2;

psi = E/sd2 + L/sA2;

logpold = logNA-0.5*psi;

m(1,1) = mold;

m(1,2) = 0.0;

m(1,3) = dimold;

elseif( dimold==2 )

mold = zeros(2,1);

mold(1) = random('Uniform',0.5,1.5);

mold(2) = random('Uniform',0.1,1.0);

dpre = mold(1)*exp(-0.5*((x-xmid).^2)/mold(2));

e = dobs-dpre;

E = e'*e;

l = mold - mpA;

L = l'*l;

psi = E/sd2 + L/spA2;

logpold = logNpA-0.5*psi;

m(1,1) = mold(1);

m(1,2) = mold(2);

m(1,3) = dimold;

end

% random verctors u are normal with these variance and normlizations

su1 = 0.2;

su12 = su1^2;

Nu1 = 1/(sqrt(2.0*pi)*su1);

logNu1 = log(Nu1);

su2 = 0.2;

su22 = su2^2;

Nu2 = 1/(sqrt(2.0*pi)*su2);

logNu2 = log(Nu2);

su3 = 1.0;

su32 = su3^2;

Nu3 = 1/(sqrt(2.0*pi)*su3);

logNu3 = log(Nu3);

% patch to avoid variace going negative

cutoff=0.001;

% Transdimensional Metropolis Hastings

for myiter =(2:Niter)

% switch dimensions every doswitch iterations

if( mod(myiter,doswitch)==0 )

if( dimnew == 1 )

dimnew = 2;

elseif ( dimnew == 2 )

dimnew = 1;

end

% three random vectors u1, u2, u3 with transformation

% for the dim-1 to dim-2 transformation

% m1' = [ 1 1 0 0 ] [m1] (m1' = m1 + u1)

% m2' = [ 0 0 0 1 ] [u1] (m2' = u3)

% u1' = [ 0 1 0 0 ] [u2] (u1' = u1)

% u2' = [ 0 0 1 0 ] [u3] (u2' = u2)

% which has a Jacobian of unity

u1 = random('Normal',0.0,su1); % g1

u2 = random('Normal',0.0,su2); % g2

u3 = random('Normal',0.0,su3); % g3

if( (dimold==1) && (dimnew==1) )

mnew = mold + u1;

dpre = mnew*sin(pi*(x-xmin)/(xmax-xmin));

e = dobs-dpre;

E = e'*e;

l = mold - mA;

L = l^2;

psi = E/sd2 + L/sA2;

logpnew = logNA-0.5*psi;

logalpha = logpnew-logpold;

elseif( (dimold==1) && (dimnew==2) )

% pnew = p12' g1' g2'

% pold = p1 g1 g2 g3

% but since g1=g1' and g2=g2'

% so pnew/pold = p12' / (p1 g3)

logg3 = logNu3 - 0.5*(u3^2)/su32;

mnew = zeros(2,1);

mnew(1)= mold(1) + u1;

mnew(2)= u3;

if( mnew(2)<cutoff )

mnew(2)=cutoff;

end

dpre = mnew(1)*exp(-0.5*((x-xmid).^2)/mnew(2));

e = dobs-dpre;

E = e'*e;

l = mnew - mpA;

L = l'*l;

psi = E/sd2 + L/spA2;

logpnew = logNpA-0.5*psi;

logalpha = logpnew - (logpold + logg3);

elseif( (dimold==2) && (dimnew==2) )

mnew = zeros(2,1);

mnew(1)= mold(1) + u1;

mnew(2)= mold(2) + u2;

if( mnew(2)<cutoff )

mnew(2)=cutoff;

end

dpre = mnew(1)*exp(-0.5*((x-xmid).^2)/mnew(2));

e = dobs-dpre;

E = e'*e;

l = mnew - mpA;

L = l'*l;

psi = E/sd2 + L/spA2;

logpnew = logNpA - 0.5*psi;

logalpha = logpnew - logpold;

elseif( (dimold==2) && (dimnew==1) )

mnew = mold(1) - u1;

up3 = mold(2);

dpre = mnew*sin(pi*(x-xmin)/(xmax-xmin));

e = dobs-dpre;

E = e'*e;

l = mnew - mA;

L = l^2;

psi = E/sd2 + L/sA2;

logpnew = logNA - 0.5*psi;

loggp3 = logNu3 - 0.5*(up3^2)/su32;

logalpha = (loggp3+logpnew)-logpold;

end

if( logalpha > 0.0)

dimold = dimnew;

mold = mnew;

logpold = logpnew;

else

alpha = exp(logalpha);

r = random('Uniform',0.0,1.0);

if( alpha > r )

dimold = dimnew;

mold = mnew;

logpold = logpnew;

else

dimnew = dimold;

mnew = mold;

logpnew = logpold;

end

if( dimnew == 1 )

m(1,myiter)=mnew;

m(2,myiter)=0.0;

m(3,myiter)=dimnew;

elseif( dimnew == 2 )

m(1,myiter)=mnew(1);

m(2,myiter)=mnew(2);

m(3,myiter)=dimnew;

end

% keep only values past burn in

m = m(1:3,Nburnt+1:Niter);

[i, Niter] = size(m);

if( dimtrue == 1)

fprintf('true dim=1, m = %.4f\n', mtrue);

elseif( dimtrue == 2) % model 1, gaussian amplitude mp1 and variance mp2

fprintf('true dim=2, m = %.4f %.4f\n', mptrue(1), mptrue(2) );

end

true dim=1, m = 1.0000

Nbins = 41;

binmin = 0.0;

binmax = 2.0;

Dbin = (binmax-binmin)/(Nbins-1);

hbin = binmin+Dbin*[0:Nbins-1]';

H1 = zeros(Nbins,1);

H1count=0;

m1sum = 0.0;

for myiter = (1:Niter)

if( m(3,myiter)==1 )

m1 = m(1,myiter);

m1sum = m1sum+m1;

k = floor((m1-binmin)/Dbin)+1;

if( (k<1) || (k>Nbins) )

continue;

end

H1(k) = H1(k)+1;

H1count = H1count+1;

end

if( H1count>0 )

m1est = m1sum/H1count;

fprintf("dim=1 mest = %.4f\n", m1est );

else

fprintf("no dim=1 members\n");

end

dim=1 mest = 0.9878

Nbins = 41;

bin1min = 0.0;

bin1max = 2.0;

bin2min = 0.0;

bin2max = 0.1;

Dbin1 = (bin1max-bin1min)/(Nbins-1);

Dbin2 = (bin2max-bin2min)/(Nbins-1);

hbin1 = bin1min+Dbin1*[0:Nbins-1]';

hbin2 = bin2min+Dbin2*[0:Nbins-1]';

H2 = zeros(Nbins,Nbins);

H2count=0;

m1sum = 0.0;

m2sum = 0.0;

for myiter = (1:Niter)

if( m(3,myiter)==2 )

m1 = m(1,myiter);

m2 = m(2,myiter);

m1sum = m1sum + m1;

m2sum = m2sum + m2;

k1 = floor((m1-bin1min)/Dbin1)+1;

k2 = floor((m2-bin2min)/Dbin2)+1;

if( (k1<1) || (k1>Nbins) || (k2<1) || (k2>Nbins) )

continue;

end

H2(k1,k2) = H2(k1,k2)+1;

H2count = H2count+1;

end

if( H2count>0 )

m1est = m1sum/H2count;

m2est = m2sum/H2count;

fprintf('dim=2 mest = %.4f %.4f\n', m1est,m2est );

else

fprintf('no dim=2 members\n');

end

dim=2 mest = 1.0588 0.0582

fprintf('percent of dim=1 %.4f\n', 100.0*H1count/(H1count+H2count));

percent of dim=1 99.5350

gda_draw(H1,H2);

fprintf('Caption: Graphical depiction of transdimensional pdf, with component 1-d\n');

Caption: Graphical depiction of transdimensional pdf, with component 1-d

fprintf('pdf at left and component 2-d pdf at right\n');

pdf at left and component 2-d pdf at right

% gdama12_06a

% Using the extended Metropolis-Hastibgs algorith to sample the a

% trans-dimensional pdf that simultageously fits data to a sinusoid

% or a bell corve

% p(m|d) = [sinusoid, 1 model parameter] + [bell curve, 2 model parameters]

clear all;

fprintf('gdama12_06a (with true dimension of 2)\n');

gdama12_06a (with true dimension of 2)

% tunable parameter:

% true dimension, can be either 1 or 2

dimtrue = 2;

% x-axis

Nx=41;

xmin = 0.0;

xmax = 1.0;

xmid = (xmax-xmin)/2.0;

Dx = (xmax-xmin)/(Nx-1);

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

% model 1, sine wave of amplitude m1

if( dimtrue == 1)

mtrue = 1.0;

dtrue = mtrue*sin(pi*(x-xmin)/(xmax-xmin));

smnt2 = 0.2^2;

dnottrue = mtrue*exp(-0.5*((x-xmid).^2)/smnt2);

elseif( dimtrue == 2) % model 1, gaussian amplitude mp1 and variance mp2

mptrue = [1.0, 0.2^2 ]';

dtrue = mptrue(1)*exp(-0.5*((x-xmid).^2)/mptrue(2));

dnottrue = mptrue(1)*sin(pi*(x-xmin)/(xmax-xmin));

end

% prior variance of data

sd = 0.2;

sd2 = sd^2;

dobs = dtrue + random('Normal',0.0,sd,Nx,1);

% plot data

figure();

clf;

hold on;

axis( [xmin, xmax, 0, 2 ] );

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

plot(x,dnottrue,'g-','LineWidth', 3);

plot(x,dobs,'r.','LineWidth', 2);

xlabel('z');

ylabel('d');

fprintf('Caption: true data curve (black) and alternative data curve (green)\n');

Caption: true data curve (black) and alternative data curve (green)

fprintf('witn observed data (red) superimposed\n');

witn observed data (red) superimposed

% dim=1

% p(m|d) = p(d|m) pA(m) / p(d)

% dim=2

% pp(mp1,mp2|d) = pp(d|mp1,mp2) ppA(m1,m2) / p(d)

% note that p(d) cancels for all ratios, and can be ignored

% note that normalization p(d|m) and pp(d|mp1,mp2) is the same

% and so cancels for all rations and can be ignored

% pA(m) assumed to be normal

% pA(m) = NA exp(-0.5*(m-mA)**2/sA2) with NA = 1/sqrt(2*pi*sA2)

% ppA(mp1,mp2) asumed to be normal

% ppA(mp1,mp2) = NpA exp(-0.5*((mp1-mpA1)**2+(mp2-mpA2)**2)/spA2)

% with NpA = 1/(2*pi*spA2)

% prior information

% dim=1

mA = 0.9;

sA = 1.0;

sA2 = sA^2;

NA = 1/sqrt(2*pi*sA2);

logNA = log(NA);

% dim=2

mpA = [0.9,0.4^2]';

spA = 1.0;

spA2 = spA^2;

NpA = 1/(2*pi*spA2);

logNpA = log(NpA);

% Metropolis Hastings

Nburnt=10000; % burn in

Niter=1000000+Nburnt; % links in Markov chain

doswitch=50; % switch dimensions

% stores links in chain

m = zeros(3,Niter); % m[3,.] 1 or 2 signifies dimension

dimold = 2; % starting dimension

dimnew = dimold;

% assign first link in chain

if( dimold==1 )

mold = random('Uniform',0.5,1.5);

dpre = mold*sin(pi*(x-xmin)/(xmax-xmin));

e = dobs-dpre;

E = e'*e;

l = mold - mA;

L = l^2;

psi = E/sd2 + L/sA2;

logpold = logNA-0.5*psi;

m(1,1) = mold;

m(1,2) = 0.0;

m(1,3) = dimold;

elseif( dimold==2 )

mold = zeros(2,1);

mold(1) = random('Uniform',0.5,1.5);

mold(2) = random('Uniform',0.1,1.0);

dpre = mold(1)*exp(-0.5*((x-xmid).^2)/mold(2));

e = dobs-dpre;

E = e'*e;

l = mold - mpA;

L = l'*l;

psi = E/sd2 + L/spA2;

logpold = logNpA-0.5*psi;

m(1,1) = mold(1);

m(1,2) = mold(2);

m(1,3) = dimold;

end

% random verctors u are normal with these variance and normlizations

su1 = 0.2;

su12 = su1^2;

Nu1 = 1/(sqrt(2.0*pi)*su1);

logNu1 = log(Nu1);

su2 = 0.2;

su22 = su2^2;

Nu2 = 1/(sqrt(2.0*pi)*su2);

logNu2 = log(Nu2);

su3 = 1.0;

su32 = su3^2;

Nu3 = 1/(sqrt(2.0*pi)*su3);

logNu3 = log(Nu3);

% patch to avoid variace going negative

cutoff=0.001;

% Transdimensional Metropolis Hastings

for myiter =(2:Niter)

% switch dimensions every doswitch iterations

if( mod(myiter,doswitch)==0 )

if( dimnew == 1 )

dimnew = 2;

elseif ( dimnew == 2 )

dimnew = 1;

end

% three random vectors u1, u2, u3 with transformation

% for the dim-1 to dim-2 transformation

% m1' = [ 1 1 0 0 ] [m1] (m1' = m1 + u1)

% m2' = [ 0 0 0 1 ] [u1] (m2' = u3)

% u1' = [ 0 1 0 0 ] [u2] (u1' = u1)

% u2' = [ 0 0 1 0 ] [u3] (u2' = u2)

% which has a Jacobian of unity

u1 = random('Normal',0.0,su1); % g1

u2 = random('Normal',0.0,su2); % g2

u3 = random('Normal',0.0,su3); % g3

if( (dimold==1) && (dimnew==1) )

mnew = mold + u1;

dpre = mnew*sin(pi*(x-xmin)/(xmax-xmin));

e = dobs-dpre;

E = e'*e;

l = mold - mA;

L = l^2;

psi = E/sd2 + L/sA2;

logpnew = logNA-0.5*psi;

logalpha = logpnew-logpold;

elseif( (dimold==1) && (dimnew==2) )

% pnew = p12' g1' g2'

% pold = p1 g1 g2 g3

% but since g1=g1' and g2=g2'

% so pnew/pold = p12' / (p1 g3)

logg3 = logNu3 - 0.5*(u3^2)/su32;

mnew = zeros(2,1);

mnew(1)= mold(1) + u1;

mnew(2)= u3;

if( mnew(2)<cutoff )

mnew(2)=cutoff;

end

dpre = mnew(1)*exp(-0.5*((x-xmid).^2)/mnew(2));

e = dobs-dpre;

E = e'*e;

l = mnew - mpA;

L = l'*l;

psi = E/sd2 + L/spA2;

logpnew = logNpA-0.5*psi;

logalpha = logpnew - (logpold + logg3);

elseif( (dimold==2) && (dimnew==2) )

mnew = zeros(2,1);

mnew(1)= mold(1) + u1;

mnew(2)= mold(2) + u2;

if( mnew(2)<cutoff )

mnew(2)=cutoff;

end

dpre = mnew(1)*exp(-0.5*((x-xmid).^2)/mnew(2));

e = dobs-dpre;

E = e'*e;

l = mnew - mpA;

L = l'*l;

psi = E/sd2 + L/spA2;

logpnew = logNpA - 0.5*psi;

logalpha = logpnew - logpold;

elseif( (dimold==2) && (dimnew==1) )

mnew = mold(1) - u1;

up3 = mold(2);

dpre = mnew*sin(pi*(x-xmin)/(xmax-xmin));

e = dobs-dpre;

E = e'*e;

l = mnew - mA;

L = l^2;

psi = E/sd2 + L/sA2;

logpnew = logNA - 0.5*psi;

loggp3 = logNu3 - 0.5*(up3^2)/su32;

logalpha = (loggp3+logpnew)-logpold;

end

if( logalpha > 0.0)

dimold = dimnew;

mold = mnew;

logpold = logpnew;

else

alpha = exp(logalpha);

r = random('Uniform',0.0,1.0);

if( alpha > r )

dimold = dimnew;

mold = mnew;

logpold = logpnew;

else

dimnew = dimold;

mnew = mold;

logpnew = logpold;

end

if( dimnew == 1 )

m(1,myiter)=mnew;

m(2,myiter)=0.0;

m(3,myiter)=dimnew;

elseif( dimnew == 2 )

m(1,myiter)=mnew(1);

m(2,myiter)=mnew(2);

m(3,myiter)=dimnew;

end

% keep only values past burn in

m = m(1:3,Nburnt+1:Niter);

[i, Niter] = size(m);

if( dimtrue == 1)

fprintf('true dim=1, m = %.4f\n', mtrue);

elseif( dimtrue == 2) % model 1, gaussian amplitude mp1 and variance mp2

fprintf('true dim=2, m = %.4f %.4f\n', mptrue(1), mptrue(2) );

end

true dim=2, m = 1.0000 0.0400

Nbins = 41;

binmin = 0.0;

binmax = 2.0;

Dbin = (binmax-binmin)/(Nbins-1);

hbin = binmin+Dbin*[0:Nbins-1]';

H1 = zeros(Nbins,1);

H1count=0;

m1sum = 0.0;

for myiter = (1:Niter)

if( m(3,myiter)==1 )

m1 = m(1,myiter);

m1sum = m1sum+m1;

k = floor((m1-binmin)/Dbin)+1;

if( (k<1) || (k>Nbins) )

continue;

end

H1(k) = H1(k)+1;

H1count = H1count+1;

end

if( H1count>0 )

m1est = m1sum/H1count;

fprintf("dim=1 mest = %.4f\n", m1est );

else

fprintf("no dim=1 members\n");

end

dim=1 mest = 0.7581

Nbins = 41;

bin1min = 0.0;

bin1max = 2.0;

bin2min = 0.0;

bin2max = 0.1;

Dbin1 = (bin1max-bin1min)/(Nbins-1);

Dbin2 = (bin2max-bin2min)/(Nbins-1);

hbin1 = bin1min+Dbin1*[0:Nbins-1]';

hbin2 = bin2min+Dbin2*[0:Nbins-1]';

H2 = zeros(Nbins,Nbins);

H2count=0;

m1sum = 0.0;

m2sum = 0.0;

for myiter = (1:Niter)

if( m(3,myiter)==2 )

m1 = m(1,myiter);

m2 = m(2,myiter);

m1sum = m1sum + m1;

m2sum = m2sum + m2;

k1 = floor((m1-bin1min)/Dbin1)+1;

k2 = floor((m2-bin2min)/Dbin2)+1;

if( (k1<1) || (k1>Nbins) || (k2<1) || (k2>Nbins) )

continue;

end

H2(k1,k2) = H2(k1,k2)+1;

H2count = H2count+1;

end

if( H2count>0 )

m1est = m1sum/H2count;

m2est = m2sum/H2count;

fprintf('dim=2 mest = %.4f %.4f\n', m1est,m2est );

else

fprintf('no dim=2 members\n');

end

dim=2 mest = 0.9467 0.0377

fprintf('percent of dim=1 %.4f\n', 100.0*H1count/(H1count+H2count));

percent of dim=1 1.5000

gda_draw(H1,H2);

fprintf('Caption: Graphical depiction of transdimensional pdf, with component 1-d\n');

Caption: Graphical depiction of transdimensional pdf, with component 1-d

fprintf('pdf at left and component 2-d pdf at right\n');

pdf at left and component 2-d pdf at right

% gdama12_07

% Using the extended Metropolis-Hastings algorith to sample the

% pdf associated with the Laplace Transform problem, in which the

% model is trans-dimensional, with 1, 2, 4 ... layers.

% The scriot uses a fun transformation and its

% inverse for layered media

% that presumes the nummber of layers is a power of two.

% forward tranformation:

% Each old layer is split into N = 2 or 4 or 8 ... new layers,

% with N=Mnew/Mold and with new layers having same value os old,

% plus random variation

% [mp;up] = T [m;u] where m are layer values, u are random variation

% invere transformation:

% Each group of N old Layers are grouped into one new layer

% with new layer having same value as the mean of the old group,

% plus random variation

% [m;u] = Ti [mp;up] where m are layer values, u are random variation

% definitions Mold, Mnew and N=Mnew/Mold, K=Mnew-Mold, J=K/Mold

% the matrix T is 2*Mnew by 2*Mnew and is broken into several blocks

% Block A, in the upper left is Mnew by Mnew, and is broken into 2 sub-blocks

% Block AL on the upper left is Mnew by Mold

% divided into (Mold groups of N rows) and Mold columns

% each group has a column of 1's, moving one place to right

% Block AR on the upper right, is Mnew by K

% it is divided into Mold by Mold sub-blocks that are N by J

% only the diagonal sub-blocks are populated

% the top row of each diagonal sub-block is populated by -1's

% the rest have 1's along their diagonal,starting flush left

% Block B, on the upper right is Mnew by Mnew

% it is the identity matrix

% Block C, on the lower left is Mnew by Mnew and is unpopulated

% Block D, on the lower left is Mnew by Mnew and is an identify matrix

% for no change of dimension, there is only a forward transformation

% with random noise being added to each model parameter

clear all;

fprintf('gdama12_07\n');

gdama12_07

% tunable parameters (beware: they are carefully tuned!)

% sigma of data = sigmad_factor * datamax

sigmad_factor = 0.005;

sigmad_factor = 0.006;

% sigma of prior model

sigmam = 5.0;

sigmam2 = sigmam^2;

% random variation of models

sigmau = 0.01;

% for error checking

maxE = 1.0E-6;

Ecount=0;

% this codebuilds and checks the transformation

dmax = 4; % maximum dimension, with 0<=dim<=dmax

Mmax = 2^dmax; % maximum number of layers

% transformation and its invere (only dimension changes considered here)

% as the first domension is zero, dimensions are misregisterd by 1

% in the arrays. My convention is that the dimension variable always

% is the actual dimension (0<=dim<=dmax) and that one is added in

% the array reference

T = zeros(dmax+1,dmax+1,2*Mmax,2*Mmax);

% Jacobian

Jac = zeros(dmax+1,dmax+1);

% build transformations

for oldpower = (0:dmax)

Mold = 2^oldpower;

for newpower = (oldpower+1:dmax)

Mnew = 2^newpower;

N = floor(Mnew/Mold);

K = Mnew-Mold;

J = floor(K/Mold);

TA = zeros(Mnew,Mnew);

TB = eye(Mnew);

TC = zeros(Mnew,Mnew);

TD = eye(Mnew);

% AL block

for i = (1:Mold)

for j = (1:N)

TA((i-1)*N+j,i)=1.0;

end

% AR block

for i =(1:Mold)

% populate top row with -1's

for k = (1:J)

TA((i-1)*N+1,Mold+(i-1)*J+k)=-1.0;

end

for k = (1:min(N,J))

TA((i-1)*N+k+1,Mold+(i-1)*J+k)=1.0;

end

To = [[TA,TB];[TC,TD]];

Toi = inv(To);

% T is always 2*Mmax by 2*Mmax

% initialized to diagional matrix

T(oldpower+1,newpower+1,1:2*Mmax,1:2*Mmax) = eye(2*Mmax);

T(newpower+1,oldpower+1,1:2*Mmax,1:2*Mmax) = eye(2*Mmax);

% then set upper-left block

T(oldpower+1,newpower+1,1:2*Mnew,1:2*Mnew) = To;

T(newpower+1,oldpower+1,1:2*Mnew,1:2*Mnew) = Toi;

% Jacobian

Jac(newpower+1,oldpower+1) = abs(det(To));

Jac(oldpower+1,newpower+1) = abs(det(Toi));

% check crital part of forward transformation, that it preserve value

m = random('Uniform',0.0,1.0,Mold,1);

u = zeros(K+Mnew,1);

mu = [m;u];

mup = To*mu;

E = 0.0;

for i = (1:Mold)

for j = (1:N)

e = m(i) - mup((i-1)*N+j);

E = E+e^2;

end

if( E>maxE )

Ecount=Ecount+1;

end

% check critcal part of inverse transformation, that it preserve mean value

mp = random('Uniform',0.0,1.0,Mnew,1);

up = zeros(Mnew,1);

mup = [mp;up];

mu = Toi*mup;

E = 0.0;

for i = (1:Mold)

s=0.0;

for j = (1:N)

s = s + mup((i-1)*N+j);

end

e = mu(i) - s/N;

E = E+e^2;

end

if( E>maxE )

Ecount=Ecount+1;

end

% check that the mp-up part of Ti has at least 1 nonzero value per row

% so that the resulting m's have random variation

for i = (1:Mold)

Nonzero=0;

for j = (Mnew+1:2*Mnew)

if (abs(Toi(i,j)) > maxE )

Nonzero=Nonzero+1;

end

if (Nonzero==0)

Ecount=Ecount+1;

end

% no change in dimension

for pof2 =(0:dmax)

To = zeros(2*Mmax,2*Mmax);

M = 2^pof2;

for i = (1:M)

To(i,i) = 1.0;

To(i,i+M) = 1.0;

end

for i = (M+1:2*Mmax)

To(i,i) = 1.0;

end

Jo = abs(det(To));

T(pof2+1,pof2+1,1:2*Mmax,1:2*Mmax) = To;

Jac(pof2+1,pof2+1) = Jo;

end

fprintf('Transformations built for dmax=%d Mmax=%d: number of errors: %d\n', dmax, Mmax, Ecount);

Transformations built for dmax=4 Mmax=16: number of errors: 0

fprintf('\n');

dold = 2;

dnew = 2;

Mold = 2^dold;

Mnew = 2^dnew;

fprintf('Example: dold=%d bnew=%d Mold=%d Mnew=%d\n', dold, dnew, Mold, Mnew);

Example: dold=2 bnew=2 Mold=4 Mnew=4

fprintf('Jacobian: %.4f\n', Jac(dold+1,dnew+1) );

Jacobian: 1.0000

fprintf('Transformation:');

Transformation:

for i = (1:2*Mnew)

line = '';

for j = (1:2*Mnew)

mystr = sprintf('%2.0f ', T(dold+1,dnew+1,i,j) );

line = [line,mystr];

end

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

end

0  0  0  1  0  0  0 
1  0  0  0  1  0  0 
0  1  0  0  0  1  0 
0  0  1  0  0  0  1 
0  0  0  1  0  0  0 
0  0  0  0  1  0  0 
0  0  0  0  0  1  0 
0  0  0  0  0  0  1 

fprintf('\n');

dold = 1;

dnew = 2;

Mold = 2^dold;

Mnew = 2^dnew;

fprintf('Example: dold=%d bnew=%d Mold=%d Mnew=%d\n', dold, dnew, Mold, Mnew);

Example: dold=1 bnew=2 Mold=2 Mnew=4

fprintf('Jacobian: %.4f\n',Jac(dold, dnew));

Jacobian: 0.5000

fprintf('Transformation:');

Transformation:

for i = (1:2*Mnew)

line = '';

for j = (1:2*Mnew)

mystr=sprintf('%2.0f ', T(dold+1,dnew+1,i,j) );

line = [line,mystr];

end

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

end

0 -1  0  1  0  0  0 
0  1  0  0  1  0  0 
1  0 -1  0  0  1  0 
1  0  1  0  0  0  1 
0  0  0  1  0  0  0 
0  0  0  0  1  0  0 
0  0  0  0  0  1  0 
0  0  0  0  0  0  1 

fprintf('\n');

fprintf('Inverse Transformation:\n')

Inverse Transformation:

for i = (1:2*Mnew)

line = '';

for j = (1:2*Mnew)

tmp = T(dnew+1,dold+1,i,j);

if( abs(tmp) < 1.e-6 )

tmp=0;

end

mystr=sprintf('%2.0f ', tmp );

line = [line,mystr];

end

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

end

0  0  0 -0 -0  0  0 
0  0  0  0  0 -0 -0 
-0  0  0  0  0 -0  0  0 
0 -0  0  0  0  0 -0 
0  0  0  1  0  0  0 
0  0  0  0  1  0  0 
0  0  0  0  0  1  0 
0  0  0  0  0  0  1 

fprintf('\n');

fprintf('Transformation * Inverse Transformation\n')

Transformation * Inverse Transformation

P = squeeze(T(dold+1,dnew+1,1:2*Mnew,1:2*Mnew)) * squeeze(T(dnew+1,dold+1,1:2*Mnew,1:2*Mnew));

for i = (1:2*Mnew)

line = '';

for j = (1:2*Mnew)

mystr=sprintf('%2.0f ', P(i,j) );

line = [line,mystr];

end

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

end

0  0  0  0  0  0  0 
1  0  0  0  0  0  0 
0  1  0  0  0  0  0 
0  0  1  0  0  0  0 
0  0  0  1  0  0  0 
0  0  0  0  1  0  0 
0  0  0  0  0  1  0 
0  0  0  0  0  0  1 

fprintf('\n');

% I think this easier to understand if I use labels 's' for small 'b' for big

% then the vectors is [ms, ub] where ms is length Ms and the other [mb, us]

% where mb is length Mb. Also, let K = Mb-Ms.

% Case 1: Small to Big (dold<dnew). The you start with [ms, ub] where you supply ub

% by generating 2*Mb-Ms random numbers. Then you compute [mb, us] via the trasnformation

% [mb, us] = T[dold,dnew,0:2*Mnew,0:2*Mnew] [ms, ub]

% The acceptance function is

% p(mb) Jsb gb(ub_K+1) ... gb(ub_2*Mmax)

% Asb = -------------------------------- ----------------------------

% p(ms) gs(us_1) ... gs(us_K) gs(us_K+1) ... gs(us_2*Mmax)

% where the right hand fraction is unity, as gs(us_i) is literlly the

% same function as gb(ub_i).

% Case 2: Big to Small (dnew<dold). The you start with [mb, us] where you supply us

% by generating Mb random numbers. Then you compute [mb, us] via the trasnformation

% [mb, us] = T[dnew,dold,0:2*Mnew,0:2*Mnew] [ms, ub]

% The acceptance function is

% p(ms) gs(us_1) ... gs(us_K) Jbs gs(us_K+1) ... gs(us_2*Mmax)

% Abs = -------------------------------- ----------------------------

% p(mb) gb(ub_K+1) ... gb(ub_2*Mmax)

% where the right hand fraction is unity. Note that the two A's are exactly

% reciprocals, as Jsb=1/Jbs.

% Case 3: Same dimension (d=dnew=dold). The you start with [m0, u0] where you supply u0

% by generating M=Mb=Ms random numbers. Then you compute [mp, up] via the trasnformation

% [mp, up] = T[d,d,0:2*M,0:2*M] [m0, u0]

% The acceptance function is

% p(mp)

% A = -----

% p(m0)

% as the Jacobian is unity and all the g's cancel..

% prior for number of dimensions; note misregistration by 1

i = [0:dmax]';

cD = 0.2;

PD = exp(-cD*i);

PD = PD/sum(PD);

logPD = log(PD);

% synthetic data

G = zeros(Mmax,Mmax);

cmax = Mmax/20.0;

c = cmax*[0:Mmax-1]'/(Mmax-1);

Dz = 1.0;

zmax=Dz*(Mmax-1);

z = zmax*[0:Mmax-1]'/(Mmax-1);

for i = (1:Mmax)

r = exp(-c(i)*z)';

G(i,1:Mmax) = r;

end

dtrue = 1;

Mtrue = 2^dtrue;

mtrue = [0.5, 1.5]';

mtrueext = [mtrue; zeros(Mmax-Mtrue,1)];

To = squeeze(T(dtrue+1,dmax+1,1:Mmax,1:Mmax));

mint = To*mtrueext;

N = Mmax;

datatrue = G*mint;

datamax = max(abs(datatrue));

sigmad = sigmad_factor*datamax;

sigmad2 = sigmad^2;

logsigmad = log(sigmad);

logtwopi = log(2.0*pi);

dataobs = datatrue + random('Normal',0.0,sigmad,Mmax,1);

gda_draw(G,'caption G',mint,'caption m','=',datatrue,'caption d');

fprintf('Caption: Graphical depiction of data equation\n');

Caption: Graphical depiction of data equation

% prior model

sigmam2 = sigmam^2;

logsigmam = log(sigmam);

mpri = ones(Mtrue,1);

mpriext = [mpri; zeros(Mmax-Mtrue,1)];

% random variation of models

sigmau2 = sigmau^2;

logsigmau = log(sigmau);

% starting model

dold = 0;

Mold = 2^dold;

mold = ones(Mold,1);

moldext = [mold; zeros(Mmax-Mold,1)];

To = squeeze(T(dold+1,dmax+1,1:Mmax,1:Mmax));

moldint = To*moldext;

dataold = G*moldint;

eold = dataobs - dataold;

Eold = eold'*eold;

logpEold = -logtwopi-Mmax*logsigmad-0.5*Eold/sigmad2;

lold = mpriext-moldext;

Lold = (lold'*lold)/(Mmax/Mold); % (Mmax/Mold) accounts for duplicates

logpLold = -logtwopi-Mold*logsigmam-0.5*Lold/sigmam2;

% Transdimensional Metropolis Hastings

Nburnt=100000; % burn in

Niter=1000000+Nburnt; % length of chain

doswitch=50; % switch model types every doswitch iterations

% stores links in chain

msave = zeros(Mmax,Niter);

dsave = zeros(Niter,1);

% save starting model

msave(1:Mmax,1) = moldint; % save interpolated models

dsave(1) = dold; % save dimension

for myiter = (2:Niter)

% switch dimensions every doswitch iterations

if( mod(myiter,doswitch)==0 )

dnew = randi([0,dmax]);

else

dnew = dold;

end

% new model parameters

Mnew=2^dnew;

% new model from random perturbation of old

uold=random('Normal',0.0,sigmau,2*Mmax-Mold,1);

muold = [mold; uold];

To = squeeze(T(dold+1,dnew+1,1:2*Mmax,1:2*Mmax));

Jo = Jac(dold+1,dnew+1);

munew = To*muold;

mnew = munew(1:Mnew);

unew = munew(Mnew+1:2*Mmax);

if( dnew>dold )

loggstuff = log(Jo);

K=Mnew-Mold;

for k = (1:K)

% in numerator, so minus

loggstuff = loggstuff - logtwopi - logsigmau - 0.5*(uold(k)^2)/sigmau2;

end

elseif( dold>dnew )

loggstuff = log(Jo);

K=Mold-Mnew;

for k = (1:K)

% in denominator, so positive

loggstuff = loggstuff + logtwopi + logsigmau + 0.5*(unew(k)^2)/sigmau2;

end

else

loggstuff = 0.0;

end

if( dnew < dmax )

mnewext = [mnew; zeros(Mmax-Mnew,1)];

To = squeeze(T(dnew+1,dmax+1,1:Mmax,1:Mmax));

mnewint = To*mnewext;

else

mnewext = mnew;

mnewint = mnew;

end

datanew = G*mnewint;

enew = dataobs - datanew;

Enew = enew'*enew;

logpEnew = -logtwopi-Mmax*logsigmad-0.5*Enew/sigmad2;

lnew = mpriext-mnewext;

Lnew = (lnew'*lnew)/(Mmax/Mnew); % (Mmax/Mnew) accounts for duplicates

logpLnew = -logtwopi-Mnew*logsigmam-0.5*Lnew/sigmam2;

logalpha = (logpEnew + logpLnew + loggstuff + logPD(dnew+1)) - (logpEold + logpLold + logPD(dold+1));

% error check on reversibility; the loggstuff must be equal and opposite

% (the reveribility condition must be correct for the method to work)

RCHECK=0;

RCHECKITER=151;

if(RCHECK && (myiter==RCHECKITER))

% swap old and new

loggstuffsave = loggstuff;

logalphasave = logalpha;

% save old

doldsave = dold;

Moldsave = Mold;

moldsave = mold;

uoldsave = uold;

muoldsave = muold;

moldextsave = moldext;

moldintsave = moldint;

dataoldsave = dataold;

logpEoldsave = logpEold;

logpLoldsave = logpLold;

% old becomes new

dold = dnew;

Mold = Mnew;

mold = mnew;

uold = unew;

muold = munew;

moldext = mnewext;

moldint = mnewint;

dataold = datanew;

logpEold = logpEnew;

logpLold = logpLnew;

% new becomes old

dnew = doldsave;

Mnew = Moldsave;

mnew = moldsave;

unew = uoldsave;

munew = muoldsave;

mnewext = moldextsave;

mnewint = moldintsave;

datanew = dataoldsave;

logpEnew = logpEoldsave;

logpLnew = logpLoldsave;

To = squeeze(T(dold+1,dnew+1,1:2*Mmax,1:2*Mmax));

Jo = Jac(dold+1,dnew+1);

if( dnew>dold )

loggstuff = log(Jo);

K=Mnew-Mold;

for k = (1:K)

% in numerator, so minus

loggstuff = loggstuff - logtwopi - logsigmau - 0.5*(uold(k)^2)/sigmau2;

end

elseif( dold>dnew )

loggstuff = log(Jo);

K=Mold-Mnew;

for k = (1:K)

% in denominator, so positive

loggstuff = loggstuff + logtwopi + logsigmau + 0.5*(unew(k)^2)/sigmau2;

end

else

loggstuff = 0.0;

end

if( dnew < dmax )

mnewext = [mnew; zeros(Mmax-Mnew,1)];

To = squeeze(T(dnew+1,dmax+1,1:Mmax,1:Mmax));

mnewint = To*mnewext;

else

mnewext = mnew;

mnewint = mnew;

end

datanew = G*mnewint;

enew = dataobs - datanew;

Enew = enew'*enew;

logpEnew = -logtwopi-Mmax*logsigmad-0.5*Enew/sigmad2;

lnew = mpriext-mnewext;

Lnew = (lnew'*lnew)/(Mmax/Mnew); % (Mmax/Mnew) accounts for duplicates

logpLnew = -logtwopi-Mnew*logsigmam-0.5*Lnew/sigmam2;

logalpha = (logpEnew + logpLnew + loggstuff + logPD(dnew+1)) - (logpEold + logpLold + logPD(dold+1));

fprintf('reversibility check 1: should be zero: %e\n',loggstuff+loggstuffsave);

fprintf('reversibility check 2: should be zero: %e\n',logalpha+logalphasave);

% this will only work once, so terminate script

fprintf('Exit Early Due to Error Check Being Activted with RCHECK=1');

return;

end

% standard acceptance test

if( logalpha > 0.0)

dold = dnew;

Mold = Mnew;

mold = mnew;

uold = unew;

muold = munew;

moldint = mnewint;

Eold = Enew;

Lold = Lnew;

logpEold = logpEnew;

logpLold = logpLnew;

else

alpha = exp(logalpha);

r = random('Uniform',0.0,1.0);

if( alpha > r )

dold = dnew;

Mold = Mnew;

mold = mnew;

uold = unew;

muold = munew;

moldint = mnewint;

moldext = mnewint;

Eold = Enew;

Lold = Lnew;

logpEold = logpEnew;

logpLold = logpLnew;

else

dnew = dold;

Mnew = Mold;

mnew = mold;

unew = uold;

munew = muold;

mnewint = moldint;

mnewext = moldext;

Enew = Eold;

Lnew = Lold;

logpEnew = logpEold;

logpLnew = logpLold;

end

msave(1:Mmax,myiter) = moldint; % save interpolated models

dsave(myiter) = dold; % save dimension

end

% delete burn-in

msave = msave(1:Mmax,Nburnt+1:Niter);

dsave = dsave(Nburnt+1:Niter);

[Niter, i] = size(dsave);

% Histogram of model dimensions; note misregisered by 1

Hdim = zeros(dmax+1,1);

for i =(1:Niter)

j = dsave(i);

Hdim(j+1) = Hdim(j+1)+1;

end

Hdim = Hdim/sum(Hdim);

fprintf('Histogram of model dimensions\n');

Histogram of model dimensions

for i=(0:dmax)

fprintf('%d %d\n',i,Hdim(i+1));

end

1.709500e-01
2.811000e-01
2.252000e-01
1.756990e-01
1.470510e-01

fprintf('\n');

fprintf('mean interpolated model\n');

mean interpolated model

mest= mean(msave,2);

for i=(1:Mmax)

fprintf('%d %d\n',i,mest(i));

end

6.053080e-01
6.084120e-01
6.292330e-01
6.300257e-01
7.162999e-01
7.166941e-01
7.203799e-01
7.207781e-01
1.204280e+00
1.204293e+00
1.205650e+00
1.206190e+00
1.217996e+00
1.217819e+00
1.218373e+00
1.218400e+00

fprintf('\n');

% plot histogram of dimensions

n=[0:dmax]';

figure();

clf;

hold on;

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

plot(n,PD,'k-','LineWidth',3);

plot(n,Hdim,'g-','LineWidth',3);

xlabel('dimensions, n');

ylabel('P(n)');

fprintf('Caption: Histogram of dimensions (green)\n');

Caption: Histogram of dimensions (green)

fprintf('together with prior pdf (black)\n')

together with prior pdf (black)

% Histogram of layer probabilities

Nv = 1001;

vmin = 0.0;

vmax = 2.0;

Dv = (vmax-vmin)/(Nv-1)

Dv = 0.0020

Hv = zeros(Nv,Mmax);

for myiter = (1:Niter)

mint = msave(1:Mmax,myiter);

for i = (1:Mmax)

j = floor((mint(i)-vmin)/Dv)+1;

if( (j<1) || (j>Nv) )

continue;

end

Hv(j,i) = Hv(j,i)+1;

end

% plot model probabilities

figure();

clf;

hold on;

colormap(jet);

axis( [0, zmax, vmin, vmax] );

Hmax = max(Hv(:));

caxis( [0.0, Hmax] );

imagesc( [0, zmax], [vmin, vmax], Hv );

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

colorbar();

xlabel('z');

ylabel('m(z)');

fprintf('Caption: Histogram of interpolated moddels, with\n');

Caption: Histogram of interpolated moddels, with

fprintf('mean model superimposed in white\n')

mean model superimposed in white