Implementing a markov chain for a probability dependent random walk

4 次查看(过去 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

回答(1 个)

VBBV
VBBV 2022-10-20
编辑:VBBV 2022-10-20
for n = 1-n0 :2: n0-1
  1 个评论
VBBV
VBBV 2022-10-20
From the given relation of N the inner loop index may be changed to take account of every 2 element in the range

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Mathematics 的更多信息

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by