sampleworld

sampleworld.m
Define the transition probabilities
Initialise some things
Do the sampling
Plot the resulting probabilities
Check detailed balance

sampleworld.m

Sampling from the simple state space model as an example of MCMC

From A First Course in Machine Learning Simon Rogers, August 2016 [simon.rogers@glasgow.ac.uk]

clear all; close all;

Define the transition probabilities

Each row defines the probabilities of going from a state

tran = [0 0.5 0 0 0.5;0.5 0 0.5 0 0;0 0.5 0 0.5 0;0 0 0.5 0 0.5;0.5 0 0 0.5 0];

Initialise some things

The initial state

x = 2;
% The number of samples
S = 2000;
% The probabilities
p = [0.1 0.1 0.4 0.2 0.2];
% A vector to store all sampled states
allx = zeros(S,1);

Do the sampling

for s = 1:S
    % Propose new according to transition matrix
    newx = find(rand<=cumsum(tran(x,:)),1);
    if rand<=p(newx)/p(x)
        x = newx;
    end
    allx(s) = x;
end

Plot the resulting probabilities

figure(1);hold off
[a,b] = hist(allx,1:5);
a = a./S;
c = bar([1:5],[a' p']);
set(c(1),'facecolor',[0 0 0])
set(c(2),'facecolor',[0.6 0.6 0.6])
legend('Sampled','True')

Check detailed balance

The number of 1,2 in the sequence should be the same (roughly) as the number of 2,1

sum(allx(1:end-1)==1 & allx(2:end)==2)
sum(allx(1:end-1)==2 & allx(2:end)==1)
% Ditto 2,3 and 3,2
sum(allx(1:end-1)==2 & allx(2:end)==3)
sum(allx(1:end-1)==3 & allx(2:end)==2)