Implementing a markov chain for a probability dependent random walk
1 次查看(过去 30 天)
显示 更早的评论
Consider a particle hopping on a one-dimensional lattice with lattice spacing s. The particle starts at the lattice site x0. The particle jumps left or right depending on a position dependent probability distribution. Concretely, I have been able to determine the following recursion relation for the particle motion
For a given N we have .
Here is the probability that the particle reaches position after a total of N jumps. Φ is the normal cumulative distribution function, with mean μ and standard deviation σ. is the Kronecker delta symbol.
To avoid confusion: in my implementation I renamed position to velocity.
My implementation seems to be working fine, but the one-indexing and array-indexing is throwing me off.
If I run an animation with the particle moves to the right towards the mean, and vice versa for . However, if I look at the output matrix called "distribution" I notice that all the elements are zero above column 158 and below column 55. Of course this might be due to rounding in matlab, but I became uncertain if I somehow messed up the indices in my recursion relation.
Unfortunately, my question then reduces to the mundane: are my implemented indices correct in the following loop?
for N = 2:Nmax
for n = n0-(N-1) :2: n0+(N-1)
distribution(N,n)=(1-kronDel(n0-(N-1),n))*distribution(N-1,n-1)*(1-normcdf((V(n-1)-mean)/sigma))...
+(1-kronDel(n0+(N-1),n))*distribution(N-1,n+1)*normcdf((V(n+1)-mean)/sigma);
end
end
Entire code:
%% Settings
clc
close all
clear all
set(0,'defaulttextinterpreter','latex')
set(0,'defaultAxesTickLabelInterpreter','latex');
set(0,'defaultLegendInterpreter','latex');
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultAxesFontSize',15)
N=100;
mean = 400;
sigma = 30;
s = 5;
V0 = 370;
[distribution, Velocity]=f(N,mean,sigma,s,V0);
figure(1)
for i = 1:N
plot(Velocity,distribution(i,:))
title(i)
xlabel('Velocity [m/s]')
ylabel('Probability')
pause(0.1)
drawnow
end
function [distribution_up_down_markov, Velocity] = f(Nmax,mean,sigma,s,V0)
%N=1 (0) case.
%Possible velocities
V = V0-s*Nmax:s:V0+s*Nmax;
%V(Nmax+1)=V0;
distribution = zeros(Nmax,length(V));
distribution(1,Nmax+1)=1;
n0 = Nmax + 1;
for N = 2:Nmax
for n = n0-(N-1) :2: n0+(N-1)
distribution(N,n)=(1-kronDel(n0-(N-1),n))*distribution(N-1,n-1)*(1-normcdf((V(n-1)-mean)/sigma))...
+(1-kronDel(n0+(N-1),n))*distribution(N-1,n+1)*normcdf((V(n+1)-mean)/sigma);
end
end
distribution_up_down_markov=distribution;
Velocity = V;
end
function d=kronDel(j,k)
d=j==k;
end
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!