AM13 Gaussian with bimodal velocity distribution

A Matlab script for the following example is avalable at sippi_AM13_metropolis_bimodal.m.

The GAUSSIAN and FFTMAa prior types implicitly assume a normal distribution of the model parameter.

It is however possible to change the Gaussian distribution to any shaped distribution, using a normal score transform. Note that when this is done the given semivariogram model for the FFTMA a priori model will not be reproduced. If this is a concern, then the VISIM type a priori model should be used.

The data and forward structures is identical to the one described in the previous example.

%% Load the travel time data set from ARRENAES
clear all;close all
D=load('AM13_data.mat');
options.txt='AM13';

%% SETUP DATA
id=1;
data{id}.d_obs=D.d_obs;
data{id}.d_std=D.d_std;
data{id}.Ct=D.Ct+1; % Covariance describing modeling error

% SETUP THE FORWARD MODEL USED IN INVERSION
forward.forward_function='sippi_forward_traveltime';
forward.sources=D.S;
forward.receivers=D.R;
forward.type='fat';forward.linear=1;forward.freq=0.1;

The desired distribution (the ‘target’ distribution) must be provided as a sample of the target distribution, in the data{id}.d_target distribution.

%% SETUP PRIOR
im=1;
prior{im}.type='FFTMA';
prior{im}.name='Velocity (m/ns)';
prior{im}.m0=0.145;
prior{im}.Va='.0003 Sph(6)';
dx=0.15;
prior{im}.x=[-1:dx:6];
prior{im}.y=[0:dx:13];
prior{im}.cax=[.1 .18];

% SET TARGET
N=1000;
prob_chan=0.5;
dd=.014*2;
d1=randn(1,ceil(N*(1-prob_chan)))*.01+0.145-dd; %0.1125;
d2=randn(1,ceil(N*(prob_chan)))*.01+0.145+dd; %0.155;
d_target=[d1(:);d2(:)];
prior{im}.d_target=d_target;

5 realizations from the corresponding a priori model looks like Figure figure_title compares the distribution from one realization of both prior models considered above.

As for the examples above, the a posteriori distribution can be sampled using e.g.

options.mcmc.nite=500000; % optional, default:nite=30000
options.mcmc.i_sample=500; % optional, default:i_sample=500;
options.mcmc.i_plot=1000; % optional, default:i_plot=50;
options=sippi_metropolis(data,prior,forward,options);

% plot posterior statistics
sippi_plot_posterior(options.txt);