我想對P和I進行隨機步驟模擬,使用randsample
就像下面這個簡單的步驟一樣。使用randsample的隨機樣本?
P=zeros(1,5); I=zeros(1,5)
%省事
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
%使用randsample
Pvec=[a b c d (some value for doing nothing)]*dt;
Pvec=Pvec./sum(Pvec);
s=randsample(1:5,1,'true',Pvec);
這是不正確的。你將如何有效地做到這一點?
這就是我想要做的,但我不認爲這是完全正確的......
更新與競爭羣我和基於P碼這組方程。
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 2012-02-05 05:14:24
@Jonas不,你說得對,它實際上只應該是在每次迭代中評估的單個值。 – HCAI 2012-02-05 08:22:55