Appendix B: Software
273
plot(p2(jump:end,1),d)
subplot(4,1,4)
% histogram of increments
[height,pos]=hist(d,100);
bar(pos,log(height))
The reader is encouraged to explore different control points and observe how the
resulting charts change.
B.5 Metropolis-Hastings Algorithm
The function metropolis3.m encodes the Metropolis-Hastings algorithm, discussed in Sect. 10.7. As input, the function receives the function f from which the
random numbers are sampled, the look-around factor β, the number nmax of random
numbers to return and a starting value x0. It returns an array x with random number
that are drawn from the supplied function f. The function directly implements the
algorithm discussed in Sect. 10.7.
% metropolis3.m
function x=metropolis3(f,beta,nmax,x0)
x=zeros(1,nmax);
x(1)=x0;
fx=f(x0);
for k=1:nmax-1
y=x(k)+beta*(2*rand-1); % new candidate
fy=f(y);
alpha=fy/fx;
% check if new is better
if (alpha>1)
% accept, if better
x(k+1)=y;
fx=fy;
else
% here alpha is smaller than unity
u=rand;
% get random number
if (alpha>u)
% compare with random number u
x(k+1)=y;
% also accept, if larger than u
fx=f(y);
else
x(k+1)=x(k);
% else re-use old value
end
end
end
end
The metropolis function is used to calculate some of the path integrals, discussed in
the next section.