Stochastic sample using randsample?

Question

I'd like to make a stochastic step simulation for P and I with randsample like this simple one below.

P=zeros(1,5); I=zeros(1,5)

%easy way

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

%Using randsample

    Pvec=[a b c d (some value for doing nothing)]*dt;
    Pvec=Pvec./sum(Pvec);
    s=randsample(1:5,1,'true',Pvec);

This isn't correct. How would you do it efficiently?

This is what I'm trying to do, but I don't think it's quite right...

UPDATE with competing populations I and P code Based this set of equations.

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

Answer 1

I don't think you can replicate your first code block "the easy way" with a call to randsample.

The first code block generates P and I recursively.

Whilst, randsample generates samples with or without replacement of the population: 1:5 in this case.

You may want to try randseq (requires the Bioinformatics Toolbox). In terms of efficiency, in general there is no way to vectorize a recursive operation.