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

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