PとIの確率的なステップシミュレーションをrandsample
としたいのですが、これは以下のように簡単です。ランダムサンプルを使った確率的サンプル?
P=zeros(1,5); I=zeros(1,5)
%簡単な方法これは正しくないrandsample
Pvec=[a b c d (some value for doing nothing)]*dt;
Pvec=Pvec./sum(Pvec);
s=randsample(1:5,1,'true',Pvec);
使用
for i=1:5
X=rand; dt=0.01;
a=randi(50,1);
b=randi(50,1);
c=randi(50,1);
d=randi(50,1);
if X<=a*dt,
P(i+1)=P(i+1)+1;
elseif X>a*dt && X<=(a+b)*dt
P(i+1)=P(i)-1;
elseif X>(a+b)*dt && X<=(a+b+c)*dt
I(i+1)=I(i)-1;
elseif X>(a+b+c)*dt && X<=(a+b+c+d)*dt
I(i+1)=I(i)+1;
else %do nothing
P(i+1)=P(i);
I(i+1)=P(i);
end
%。どのように効率的にやりますか? IおよびPコードは、この方程式のセットをベースに、競合集団とUPDATE
これは私がやろうとしているものですが、私はそれが非常に適切ではないと思う...
。
theta_P=0.15;delta_P=0.01;alpha_I=0.4;gamma_I=0.01;delta_I=0.005;lambda_I=0.05;
m=100; % # runs
time=10; % # Total time of simulation
dt=0.01; % # Time step
D=6000; T=10/dt;
P=zeros(m,time/dt); I=zeros(m,time/dt);
for i=1:m
for j=1:time/dt
arrivalI=alpha_I+P(i,j)*lambda_I;
lossI=I(i,j)*gamma_I+P(i,j)*I(i,j)*delta_I;
if j<=T
alpha_P=D/T;
else
alpha_P=0;
end
arrivalP=alpha_P+P(i,j)*theta_P;
lossP=P(i,j)*I(i,j)*delta_P;
X=rand;
Pvec=[arrivalI lossI arrivalP lossP]*dt;%
Pvec=Pvec./sum(Pvec);
s=randsample(1:4,1,'true',Pvec);
if s==1
I(i,j+1)=I(i,j)+1;%;
P(i,j+1)=P(i,j);
elseif s==2
I(i,j+1)=I(i,j)-1;%
P(i,j+1)=P(i,j);
elseif s==3
P(i,j+1)=P(i,j)+1;%
I(i,j+1)=I(i,j);
elseif s==4
P(i,j+1)=P(i,j)-1;%;
I(i,j+1)=I(i,j);
else
P(i,j+1)=P(i,j); %check
I(i,j+1)=I(i,j);
end
end
subplot(2,2,1:2)
%
if P(i,j)>5
loglog(abs(P(i,:)),'-r')
%
else
loglog(abs(P(i,:)),'-b')
%
end
hold on
axis([1 1e3 1 1e4])
end
あなたのステートメント '' X <= a * dt'の場合、 'a'は配列です。それは意図的なのでしょうか? – Jonas
@ Jonasあなたが正しいとは限りませんが、実際にはすべての反復で評価される単一の値であると考えられます。 – HCAI