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

   1 matlab  comment % PROGRAM theta_logistic.m
   2 matlab  comment %    Calculates by simulation the probability that a population
   3 matlab  comment %    following the theta logistic model and starting at Nc will fall
   4 matlab  comment %    below the quasi-extinction threshold Nx at or before time tmax
   5 matlab  blank
   6 matlab  comment % SIMULATION PARAMETERS
   7 matlab  comment % for butterflies (Euphydryas editha bayensis) at Jasper Ridge (population C)
   8 matlab  blank
   9 matlab  code    r=0.3458;                       % intrinsic rate of increase--Butterflies at Jasper Ridge
  10 matlab  code    K=846.017;                      % carrying capacity
  11 matlab  code    theta=1;                    % nonlinearity in density dependence
  12 matlab  code    sigma2=1.1151;          % environmental variance
  13 matlab  code    Nc=94;                          % starting population size
  14 matlab  code    Nx=20;                          % quasi-extinction threshold
  15 matlab  code    tmax=20;                        % time horizon
  16 matlab  code    NumReps=50000;      % number of replicate population trajectories
  17 matlab  blank
  18 matlab  comment % SIMULATION CODE
  19 matlab  blank
  20 matlab  code    sigma=sqrt(sigma2);
  21 matlab  code    randn('state',sum(100*clock));  % seed the random number generator
  22 matlab  blank
  23 matlab  code    N=Nc*ones(1,NumReps);   % all NumRep populations start at Nc
  24 matlab  code    NumExtant=NumReps;      % all populations are initially extant
  25 matlab  code    Extant=[NumExtant];             % vector for number of extant pops. vs. time
  26 matlab  blank
  27 matlab  code    for t=1:tmax,                           % For each future time,
  28 matlab  code            N=N.*exp( r*( 1-(N/K).^theta )...   %   the theta logistic model
  29 matlab  code            + sigma*randn(1,NumExtant) );   %   with random environmental effects.
  30 matlab  code            for i=NumExtant:-1:1,       % Then, looping over all extant populations,
  31 matlab  code                    if N(i)<=Nx,            %   if at or below quasi-extinction threshold,
  32 matlab  code                            N(i)=[];                        %   delete the population.
  33 matlab  code                    end;
  34 matlab  code            end;
  35 matlab  code            NumExtant=length(N);        %   Count remaining extant populations
  36 matlab  code            Extant=[Extant NumExtant];  %   and store the result.
  37 matlab  code    end;
  38 matlab  blank
  39 matlab  comment % OUTPUT CODE
  40 matlab  comment % ComputeS quasi-extinction probability as the fraction of replicate
  41 matlab  comment % populations that have hit the threshold by each future time,
  42 matlab  comment % and plotS quasi-extinction probability vs. time
  43 matlab  blank
  44 matlab  code    ProbExtinct=(NumReps-Extant)/NumReps;
  45 matlab  code    plot([0:tmax],ProbExtinct)
  46 matlab  code    xlabel('Years into the future');
  47 matlab  code    ylabel('Cumulative probability of quasi-extinction');
  48 matlab  code    axis([0 tmax 0 1]);
  49 matlab  blank
  50 matlab  comment % Integrate solution exactly %
  51 matlab  comment % Options=[];
  52 matlab  comment % [T,true] = ode45(@logistic,[0,20],Nc,Options,r,K,theta);
  53 matlab  comment % subplot(1,2,2)
  54 matlab  comment % plot([1:tmax],P,'r.-',T,true,'g.-')
  55 matlab  blank
  56 matlab  comment  ... This is a seriously old-school comment.