2012-02-03 13 views
1

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

これは私がやろうとしているものですが、私はそれが非常に適切ではないと思う... enter image description here

enter image description here

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 
+0

あなたのステートメント '' X <= a * dt'の場合、 'a'は配列です。それは意図的なのでしょうか? – Jonas

+0

@ Jonasあなたが正しいとは限りませんが、実際にはすべての反復で評価される単一の値であると考えられます。 – HCAI

答えて

0

私はあなたの最初のコードブロックrandsampleを呼び出して、「簡単な方法」を複製することができるとは思いません。

最初のコードブロックは、PIを再帰的に生成します。

一方、randsampleは、集団の置換の有無にかかわらずサンプルを生成する:1:5この場合、

randseq(Bioinformatics Toolboxが必要です)を試してみてください。効率の面では、一般に再帰的操作をベクトル化する方法はありません。

+0

これを見ていただきありがとうございます。基本的には、確定的な方程式を確率的に解くためのものです。あなたはrandsampleまたはrandseqを使ってmatlabでこれを見つけたことがありますか?または実際に%の簡単な方法が好きではないいくつかの他の方法ですか? – HCAI

+0

再帰的な問題の性質上、必要な出力を生成するすべてのメソッドは実質的に「簡単な方法」になります。シミュレーションを高速化する必要がある場合は、[Matlab Coder](http://www.mathworks.com/products/matlab-coder/)テクノロジの一部を確認してください。 – MattLab

関連する問題