1

我想为 P 和我做一个随机步骤模拟，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);

这是不正确的。你将如何有效地做到这一点？

这就是我正在尝试做的，但我认为这不太对... 在此处输入图像描述

更新与竞争种群 I 和 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
4

1 回答 1

0

我认为您无法通过调用randsample.

第一个代码块以P递归I方式生成。

同时，randsample生成有或没有替换总体的样本：1:5在这种情况下。

您可能想尝试randseq（需要生物信息学工具箱）。在效率方面，一般来说没有办法向量化一个递归操作。

于 2012-02-04T00:45:28.820 回答