AM13 Gaussian, accounting for modeling errors

[A Matlab script for the following example is avalable at examples/case_tomography/sippi_AM13_metropolis_modeling_error.m.]

The use of any of the forward models defined above, will be approximation to solving the perfect forward problem. This leads to a 'modeling' error as demonstrated by [HCM14]. If one has access to an optimal (but perhaps computational inefficient) forward model, and a faster (less accurate) forward model, then a Gaussian model of the modeling error caused by using the approximate, as opposed to the optima, forward model can be estimated using sippi_compute_modelization_forward_error. sippi_compute_modelization_forward_error.

SIPPI allows accounting for such modeling error through the dt and Ct fields for the data structure.

The setup of the data, prior 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 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];

% 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;

In order to compute the modeling with respect to using the 'Born' type forward model, one can define a new forward structure, here forward_full, and estimate a Gaussian model for the modeling error using

% SETUP THE 'OPTIMAL' FORWARD MODEL
forward_full.forward_function='sippi_forward_traveltime';
forward_full.sources=D.S;
forward_full.receivers=D.R;
forward_full.type='Born';forward_full.linear=1;forward_full.freq=0.1;

% COMPUTE MODELING ERROR DUE TO USE OF forward AS OPPOSED TO forward_full
N=100;
[Ct,dt,dd]=sippi_compute_modelization_forward_error(forward_full,forward,prior,data,N);

% ASSIGN MODELING ERROR TO DATA
data{1}.dt=dt{1	};
data{1}.Ct=Ct{1};

Sampling of the posterior can proceed exactly as for the previous example, using

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);