MCMC Affine invariant ensemble sampler tutorial

This tutorial demonstrates an example of how to use the affine invariant ensemble MCMC solver in this package for model parameter estimation. This is a lower level usage in contrast to the mcmc.optimisation function. The main differences are:

  • no masking of input measurement data

  • only posterior distirbution is available (no summary statistics derived from posterior distribution)

Let’s say we have a simple monoexponential decay model:

\[S = S0 \times e^{-R_{2}^{*}t}\]

In this model, we have two parameters to be estimated: \(S0\) and \(R_{2}^{*}\).

The first thing is to create a function to generate the forward signal. Here is an example:

Note that the design of the forward function is slightly stricter for the MCMC solver. The output signal s must have a dimension of Nmeas by Nvoxel.

We can simulate the measurements using this function

%% generate some signal based on monoexponential decay
% reproducibility
seed = 5438973; rng(seed); gpurng(seed);


% set up estimation parameters; must be the same as in FWD function
modelParams = {'S0','R2star','noise'};

% define number of voxels and SNR
Nsample = 50;
SNR     = 100;

mask        = ones(1,Nsample)>0;
t           = linspace(0,40e-3,15); 
% GT
S0          = 1 + randn(1,Nsample)*0.3;
R2star      = 30 + 5*randn(1,Nsample);
% forward signal generation
pars.(modelParams{1}) = S0; 
pars.(modelParams{2}) = R2star;
S                     = Example_monoexponential_FWD_GD(pars,t);

% realistic signal with certain SNR
noise   = mean(S0) / SNR;           % estimate noise level
y       = S + noise*randn(size(S)); % add Gaussian noise

Note

Note that ‘y’ must be a [Nmeas x Nvoxels] matrix here.

To estimate \(S0\) and \(R_{2}^{*}\) from y,

  1. Set up the starting point for the estimation

pars0.(modelParams{1}) = 1 + randn(1,Nsample)*0.3;  % S0
pars0.(modelParams{2}) = 30 + 5*randn(1,Nsample);   % R2*
pars0.(modelParams{3}) = ones(1,Nsample)*0.001;     % noise

  1. Set up the model parameters and fitting boundary

% set up fitting algorithm
fitting                     = [];
% define model parameter name and fitting boundary
fitting.modelParams         = modelParams;
fitting.lb                  = [0, 0, 0.001];    % lower bound 
fitting.ub                  = [2, 50, 0.1];     % upper bound

  1. Set up optimisation setting

% Estimation algorithm setting
fitting.iteration    = 1e4;
fitting.burnin       = 0.1;     % 10% iterations
fitting.thinning     = 5;
fitting.StepSize     = 2;
fitting.Nwalker      = 50;

  1. Define the forward function

% define your forward model
modelFWD = @Example_monoexponential_FWD_GD;

  1. Define fitting weights (optional)

% equal weights
weights = [];

  1. Start the optimisation

mcmc_obj    = mcmc;
xPosterior  = mcmc_obj.goodman_weare(y,pars0,weights,fitting,modelFWD,t);

  1. Plot the estimation results

%% plot the estimation results
% compute the mean values of the posterior distribution
% xPosterior.({modelParamss{k}}) is organised in [Nvoxel,Nwalker,Nmcmcsample]
S0_mean     = mean(reshape(xPosterior.S0    ,[Nsample, prod(size(xPosterior.S0,2:3))]),2);
R2star_mean = mean(reshape(xPosterior.R2star,[Nsample, prod(size(xPosterior.R2star,2:3))]),2);

figure;
nexttile;scatter(S0,pars0.(modelParams{1}));hold on; scatter(S0,S0_mean);refline(1);
xlabel('GT'); ylabel('S0')
nexttile;scatter(R2star,pars0.(modelParams{2}));hold on; scatter(R2star,R2star_mean);refline(1)
xlabel('GT'); ylabel('R2*')
legend('Start','fitted')

The full example script can be found in here.