git.oblomov.eu Git - ohcount/blob - test/expected_dir/matlab1.m/matlab/code

   1 r=0.3458;                       % intrinsic rate of increase--Butterflies at Jasper Ridge
   2 K=846.017;                      % carrying capacity
   3 theta=1;                    % nonlinearity in density dependence
   4 sigma2=1.1151;          % environmental variance
   5 Nc=94;                          % starting population size
   6 Nx=20;                          % quasi-extinction threshold
   7 tmax=20;                        % time horizon
   8 NumReps=50000;      % number of replicate population trajectories
   9 sigma=sqrt(sigma2);
  10 randn('state',sum(100*clock));  % seed the random number generator
  11 N=Nc*ones(1,NumReps);   % all NumRep populations start at Nc
  12 NumExtant=NumReps;      % all populations are initially extant
  13 Extant=[NumExtant];             % vector for number of extant pops. vs. time
  14 for t=1:tmax,                           % For each future time,
  15 N=N.*exp( r*( 1-(N/K).^theta )...   %   the theta logistic model
  16 + sigma*randn(1,NumExtant) );   %   with random environmental effects.
  17 for i=NumExtant:-1:1,       % Then, looping over all extant populations,
  18 if N(i)<=Nx,            %   if at or below quasi-extinction threshold,
  19 N(i)=[];                        %   delete the population.
  20 end;
  21 end;
  22 NumExtant=length(N);        %   Count remaining extant populations
  23 Extant=[Extant NumExtant];  %   and store the result.
  24 end;
  25 ProbExtinct=(NumReps-Extant)/NumReps;
  26 plot([0:tmax],ProbExtinct)
  27 xlabel('Years into the future');
  28 ylabel('Cumulative probability of quasi-extinction');
  29 axis([0 tmax 0 1]);