2012-02-03 2 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 

%. 어떻게 효율적으로 할 수 있습니까?

이 내가 할 노력하고있어,하지만 난이 상당히 잘 생각하지 않습니다 ... 방정식의 집합을 집단에게 나는 경쟁과 기반 P 코드와 enter image description here

UPDATE.

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