

% this program aims at clarifying the concept of prior information and
% "samples" of a prior p.d.f.
%
% we want to generate "likely" depth dependent earth models with 
% varying layer thicknesses and densities

% Exercises:
%
% 1. Change the estimator of dispersion (s) in the log-normal pdf for ro
% and discuss the effects on the resulting samples of the model space!
%
% 2. Calculate for any setup the mean density and the mean layer thickness
% for the model samples. How close are they to the theoretically predicted
% mean values?
%
% 3. What would be the appropriate homogeneous p.d.f. for the density?
% Modify the current log-normal distribution to this pdf and discuss the
% resulting models. 
%
% 4. Somebody assumes that a constant p.d.f. describes the homogeneous pdf.
% Discuss the resulting model samples when both pdf s are constant. 
%
% 5. Define Gaussian pdf with appropriate means and variances for density
% and layer thicknesses. Calculate model samples and check their means.
%
% 6. Verify that the distribution of densities at a certain depth converges
% to (any of your) prior information as the number of model samples increases. 

ntrial=1000;    % number of trial values to test theoretical distributions


% Densities are described by a log normal distribution
nro=100;                        % number of point in ro axis
romax=25;                       % maximum density
ro=romax/nro:romax/nro:romax;   % densitiy ranges from 1-25 g/cm3
roinv=1./ro;                    % lightness
ro0=5;                          % reference density for log normal dsitribution
s=.2;                           % log normal dispersion paramter
pdfro=1./(sqrt(2*pi)*s)*roinv.*exp(-1/(2*s^2).*log10(ro/ro0).^2);   % pdf for ro
pdfro=pdfro/sum(pdfro);                                             % normalization
modelsro=pdfro*0;                                                   % field for test values

% Layer thicknesses described by 1/x pdf (positive parameter)
nld=100;                        % number of point in layer thickness axis
lmax=50;                        % maximum layer
ld=lmax/nl:lmax/nl:lmax;        % densitiy ranges from 1-25 g/cm3
pdfld=1./ld;                    % pdf for layer thickness
pdfld=pdfld/sum(pdfld);         % normalization
modelsld=pdfld*0;               % tield for test values

% Generate randomly ntrial densities and thicknesses and compare with pdf s

for i=1:ntrial,
    
% density
im=get_index_from_pdf(pdfro);
modelsro(im)=modelsro(im)+1;
% thickness
im=get_index_from_pdf(pdfld);
modelsld(im)=modelsld(im)+1;

end

% compare graphically
subplot(121),plot(ro,pdfro,'k-',ro,modelsro/max(modelsro)*max(pdfro),'r-')
xlabel(' Ro (g/cm^3) ')
ylabel(' Prob (ro) ')
title(' Theoretical (black), Samples (red) ')

subplot(122),plot(ld,pdfld,'k-',ld,modelsld/max(modelsld)*max(pdfld),'r-')
xlabel(' D (m) ')
ylabel(' Prob (D) ')
title(' Theoretical (black), Samples (red) ')
disp(' Press any key to continue ')
pause



% Generate random Earth models following the prior probabilities


nplot=5;                        % plot matrix x
mplot=5;                        % plot matrix y
nmodels=nplot*mplot;            % total number of plots
dmax=200;                       % maximum depth for layered models
id=1;                           % layer index initialisation
ii=0;                           % plot index initialisation

disp(' Calculating model samples ... ')

for iplot=1:nplot,
    for jplot=1:mplot,
        
        clear thick                             % clear thickness variable
        clear dens                              % clear density variale
        d=0;                                    % set depth to 0
        ii=ii+1;                                % plot increment
        
while d<dmax,
thick(id) =ld(get_index_from_pdf(pdfld));       % get pointer into model axis for thickness
dens(id)  =ro(get_index_from_pdf(pdfro));       % get pointer into model axis for density
d=d+thick(id);                                  % total thickness
id=id+1;                                        % layer index
end

% Plot each model
subplot(nplot,mplot,ii)
for k=1:id-1,
    if k==1, 
        plot([0 thick(1) thick(1)],[dens(1) dens(1) dens(2)],'-'), hold on
        d1=thick(1);
    elseif k<id-1
        d2=d1+thick(k);
        plot([d1 d2 d2],[dens(k) dens(k) dens(k+1)],'-'), hold on
        d1=d2;
    else
        d2=d1+thick(k);
        plot([d1 d2],[dens(k) dens(k)],'-'), hold on
        d1=d2;        
    end
end
axis([0 dmax 0 romax])
hold off
    
end
end
