gen_rws.m

Revision 1 - 6/13/07 at 2:24 pm by e2holmes

Back to revision history for gen_rws.m
This file is part of the project Generate stochastic population processes
function [simrw,rwprob]=gen_rws(truepar,decline_level,rw_type,reps,simlen)
%rw_type = 1 CRW
%rw_type = 2 RW; all process error
%rw_type = 3 all non-process error
%simlen = number of time points
%reps = number of reps
%Declare output vector
cda=[ones(reps,1)];
da=[ones(reps,1)];
cor=[ones(reps,1)];
mu = truepar(1);
sig2 = truepar(2);
sig2np = truepar(3);
simlen=simlen-1;  %need to do this so that routine returns data of length simlen


%seed random number generator
randn('state',sum(100*clock));

%generate time series
switch rw_type;
   case 1
    var_p =randn(reps,simlen+1)*sqrt(sig2);
    var_np =randn(reps,simlen+1)*sqrt(sig2np);
	 simrw =cumprod(exp(mu+var_p)')';  %this does the whole sim without walking through t's
	 simrw = simrw .* exp(var_np); %add the corruption
 case 2 
    var_pure_p =randn(reps,simlen+1)*sqrt(truepar(4));
    simrw =cumprod(exp(mu+var_pure_p)')';
   case 3 
    var_pure_np =randn(reps,simlen+1)*sqrt(truepar(5));
    simrw = exp(var_pure_np); %only corruption
 end
 
%prob at end of t
rwprob.atend = sum(simrw(:,2:end)<(decline_level*simrw(:,1)*ones(1,simlen)),1)/reps;

%prob within t
tmp=cumsum(simrw(:,2:end)<(decline_level*simrw(:,1)*ones(1,simlen)),2);
rwprob.within=sum(tmp>0,1)/reps;
Sculpin 0.2 | xhtml | problems or comments? | report bugs