quasi-extinction.nb (Mathematica 5.1 for Students

Quasi-extinction risk is defined as the probability of a population falling below a critical density. It is an important tool for assessing the future viabilty of threatened animal and plant populations. Quasi-extinction is traditionally computed using a Monte Carlo simulation of a discrete stochastic population growth model. For the continuous geometric Brownian motion model there is an analytical formula for computing quasi-extinction risk. The formula was originally derived by Ginzburg, et al. (1982; Risk Analysis 2:385--399).

In order to compute the quasi-extinction probability, I must be able to compute the cumulative distribution function for the normal distribution function. The normal distribution with a mean of zero and a standard deviation of 1 is:

In[2]:=

$General :: spell1 : Possible spelling error: new symbol name \"NormalD\" is similar to existing symbol \"Normal\". More…$

Out[2]=

In[3]:=

[Graphics:HTMLFiles/index_5.gif]

Out[3]=

The cumulative distribution is computed by integrating the normal distribution function from -infinity to y.

In[4]:=

$NormalDCDF[y_] := Integrate[NormalD, {x, -Infinity, y}]$

In[5]:=

[Graphics:HTMLFiles/index_9.gif]

Out[5]=

The formula below is for quasi-extinction defined as the first crossing time. The following parameters are defined:
    r is the mean rate of growth,
    sigma is the standard deviation of the growth rate,
    theta is the critical threshold which is expressed as a fraction of the original abundance, and
    thorizon is  the time horizon.

In[6]:=

QuasiExtinction[r_, sigma_, theta_, thorizon_] := NormalDCDF[(-theta - r thorizon)/Sqrt[thorizon sigma^2]] + Exp[(-2 r theta)/sigma^2] NormalDCDF[(-theta + r thorizon)/Sqrt[thorizon sigma^2]]

In[7]:=

Out[7]=

What is the affect of the variation in the growth rate on the quasi-extinction risk? Intuitively, it is expected that as the population growth rate becomes more variable then the quasi-extinction probability should also increase.

In[8]:=

$Plot[QuasiExtinction[0, sigma, 0.5, 10], {sigma, 0.01, 1.}, AxesLabel  {σ_r, p}]$

[Graphics:HTMLFiles/index_15.gif]

Out[8]=

Using the Plot3D function the mapping of theta and sigma to quasi-extinction risk can be explored Note this computation takes some time to compute. This possibly could be optimized by using Mathematica's built-in functions for probability distributions in the Statistics library.

In[10]:=

$Plot3D[QuasiExtinction[0, sigma, theta, 10], {theta, 0.01, 0.5}, {sigma, 0.01, 0.5}, AxesLabel {σ_r, θ, p}]$

[Graphics:HTMLFiles/index_18.gif]

Out[10]=

Created by Mathematica (July 5, 2005)