A Review of the Methodologies Used in the Derivation of Formulas for Parametric Survival Functions with Illustrative Numerical Examples

In this review paper formulas for survival functions are derived that take into account risks of deaths in early life including infancy, mid life, a random component of deaths due to accidents and deaths in older ages. The basic ideas used in the derivation of survival function for each of the components just mentioned are risk functions. Given a formula for a risk function, it is possible to derive a formula for the corresponding survival function. By using the theory of competing risks, a formula for survival function that takes into account the risks of deaths in various stages of life expressed as a product of survival functions for the risks of deaths under consideration. For many applications information on the numerical values of parameters in survival functions is not available. Consequently, rationales are developed for assigning plausible values to parameters that take into account personal ideas of an investigator may have for each stage of life. For every assignment of parameter values in the paper, a numerical version of survival functions are plotted in graphs so that an assessment of the plausibility of the chosen parameter values may be made. Also included in the paper is an application of survival functions in an experiment to make an assessment as to whether a small population of chimpanzees, or some other endangered species of animals, will have descendants that make up a surviving population 200 years into the future.


Introduction
In computer simulation experiments designed to study the evolution of age structured animal or human populations with respect to two or more genetic autosomal loci in deep time, such as hundreds or thousands of years before the present, it is almost always the case that no mortality data exists.Consequently, an investigator will need to rely on model parametric survival functions.As is well known, efforts have been made to construct parametric survival functions for nearly the past two centuries.Among the first such efforts was that of Gompertz B. Gompertz 1825.As Gompertz examined human survival data at ages 70 or more, he came to the realization that the risk of death increased exponentially for the older ages.In 1869 Makeham published a paper indicating that a positive constant should be added to Gompertz's observation to take into account deaths by accidents.In 1872, Thiele proposed a formula for the risk of death for the whole life.For further details on the proposals of Gompertz and Makeham, the book by Smith and Keyfitz The work of Thiele is presented In the book and the formula that Thiele proposed is presented on page 72.
In 1982 Mode and Busby published a paper on an eight parameter model of human mortality.In this model attention was focused on three stages of life in humans.The first stage was that for early life for ages x in the interval [0, 10 ) expressed in years.The second stage was for midlife was ages in the interval (10,30], and that for the older ages was the interval (30, k], where k > 0 was the greatest age considered in a model.For humans, this age is often chosen as k = 100.A second paper on formulation of a of parametric survival function was that of Mode and Jacobson (1984) in which age intervals mentioned above were also used.Although the formula for the parametric function used in the two papers for the same three age intervals under consideration were different, the same method of splicing the survival functions together was the same for both papers, which provided a parametric survival function for the entire age interval [0, k].Although both models fit survival functions estimated from many period life tables, they are also mildly flawed, because of small discontinuities at the splice points.The two papers just mentioned were included in the references in the paper by Gage and Mode (1993) along with several other parametric survival functions, which were tested for goodness of fit to data.
It should be emphasized, however, that in this paper the focus of attention will be based on rationals for constructing theoretical parametric survival functions that will be used in computer simulation experiments in evolutionary genetics an other fields in biology when no survival data are available.
In the paper Mode et al. 2013 a survival function related to that of Mode and Busby was used in a paper on a stochastic model of a two sex age structured population with a genetic component.A four parameter parametric survival function that was used to produce the content of the paper was based on partitioning the life span of each individual into two age intervals [0, 30) and [30, r).For this model, age x = 30 was a splice point and there was a noticeable discontinuity at that age point that appeared in age distribution estimated from data generated in Monte Carlo simulation experiments as well as those computed using the embedded deterministic model.Although it seems unlikely that the conclusions of the paper based on computer experiments would change significantly if another survival model was used to remove such discontinuities, in future experiments it would be advantageous to use a parametric survival function that did not have noticeable discontinuities.Moreover, it seems likely that it would also be beneficial to readers of the paper if these annoying discontinuities were removed.
Even though it was realized that there was need for parametric survival functions that did not have any noticeable discontinuities, no sustained efforts were made to formulate parametric survival functions with no discontinues throughout the ages considered in future computer simulation experiments.However, some time in 2015 in a search for APL functions stored in folders on flash memory sticks and other storage devises, by chance an old folder was found containing APL software designed to implement parametric survival functions with no discontinuities.As it turned out the software had been written in connection with research on aging funded by the National Institute of Aging during the late 1980s.By reading the APL code, it was realized that some potentially significant issues were not considered in the software so it was rewritten to make the output more informative as to the significance of the ideas contained in the written code.An inspection of the numerical output of the revised code led to a need to further develop the software to fix potentially flawed ideas.In the remaining sections of this paper, these potential flaws will be investigated and alternative procedures will also be investigated.

A Review of Fundamentals
There is a large literature on survival analysis, but it beyond the scope of this paper to provide review of this literature.There is also a rather large number of books on the subject, and among them are Elandt-Johnson and Johnson (1980) If a reader is interested in a more in depth account of the ideas presented in this section, it is recommended that this book be consulted.To make this paper self contained, in this section attention will be focused on the basic fundamentals and definitions that are necessary to provide a background for the contents of this paper.
Let X denote a random variable representing the life span of any individual in some biological population.Because age is measured in terms of a number in the interval R + = [0, ∞], the range of X is this interval.By definition the distribution function of the random variable X is P for all x ∈ R + , and is the probability that an individual dies in the age interval [0, x).It will be assumed that F (0) = 0 and Any new born individual in the population will be assigned the age x = 0. Let S (x) denote the probability that any individual is alive at age x > 0. Then is called the survival function, and it follows that for all x ∈ R + .The function f (x) is called the probability density function p.d. f .It also follows that for all x ∈ R + .From this equation, it follows that if one constructs or defines a probability density function f (x) that is defined for all x ∈ R + and has the property (2.7) A well known example of a pd f is that of the exponential distribution which has the form for x ∈ R + , where β > 0 is a positive parameter.The distribution function of this distribution is (2.9) and the survival function is for all x ∈ R + .From these equations, it can be seen that A central concept in the formulation of parametric survival functions is the idea of a risk function.As a step toward defining this function, consider the conditional probability and let F (x) denote the distribution of the random variable X and let f (x) denote f (x) the probability density function of the distribution of X, where x > 0 and h > 0. Then (2.12) By definition, the risk function θ (x) is defined by the limit for x > 0. An alternative definition of the risk function is From this definition, it follows that (2.15) By taking the derivative of both sides of this equation, it can be seen that the pd f has the form From the definition of the risk function, it can be seen that if f (x) is the pd f of some distribution that if a formula for f (x) is given, then a formula for the risk function may be derived.Consider, for example, the exponential distribution described above.Then for all x > 0. From this general formula connecting the survival function and the risk function, it can be seen that the survival function of the exponential distribution is as given above.
There are a number of points of view that may be implemented in connection with formulating models of parametric survival functions.The following example from reliable theory, which is used in connection with quality control in manufacturing, will provide a basis for formulating the parametric survival function presented in this paper.Suppose some device has four components, and let the random variables X k for k = 1, 2, 3, 4 taking values in the interval 0 ≤ x < ∞ denote that random life spans of the four components, .Suppose these random variables are independently distributed with survival functions S k (x) for k = 1, 2, 3, 4. Then the survival function of the system is When any component of the system fails, the system fails.Consequently, the life span of the system is given by the random variable Observe that Y > y , if and only if, X k > y for every k = 1, 2, , 3, 4. Therefore, the survival function of the system is for y > 0.
Let θ k (x) denote the risk function for component k = 1, 2, 3, 4.Then the survival function of the system has the form Therefore, from the right most formula in 2.22 it can be seen that (2.23) A procedure that has been widely used in the formulation of parametric survival functions is to define the risk function under consideration and then use expression 2.23 in deriving a symbolic form of a survival function.
The survival function of the Gompertz distribution has the form where α and β are positive parameters.As can be seen from this formula, this function has the property lim y→∞ S G (y) = 0 . (2.25) It is well known that the Gompertz survival function fits the survival function that is estimated in many period life tables of ages y ≥ 30.Throughout this paper, the Gompertz survival function will be used as a basis for formulating parametric survival functions for all ages y > 0.
Because it fits empirical survival functions for ages y ≥ 30, it follows that the function needs to be modified only for ages 0 < y < 30.As will be shown is subsequent sections of this paper, the modifying functions may not converge to zero as y → ∞.But as can be seen from the right most expression in equation 2..22, if the Gompertz survival function is one of the survival functions for the system, then In the following sections of this paper, several examples of the situation just described will be presented.

A Brief History of Proposed Risk Functions for Early Life
The form of the risk function for early life model prosed by Thiele in 1871was where the parameters α 0 > 0 and β o > 0 are positive.As can be seen from this formula µ T (x) → 0 as x → ∞.The integral of this risk function is for x > 0. The survival function in this case is by definition As can be seen from 3.2, it follows that Technically, S T (x) is not a survival function, because S T (x) does not converge to 0 as x → ∞.However, it could be used to adjust the Gompertz survival function for the ages 0 < x < 30 as indicated in 2.24 in section 2. Siler (1970).also used the risk function in 3.1 in his model of animal mortality.
The formula for risk of death in early life in this paper is the modified version that proposed by Thiele and has the form for x > 0. The integral of this risk function is When data are not available, then an investigator needs to explain the rationale used to justify his procedure for choosing numerical values of the parameters .Let p 0 denote the probability that an infant born as age x = 0 survives to age 1 year, and suppose that e −α 0 = p 0 .
(3.8) Then α 0 = − ln p 0 . (3.9) For the case of p 0 = 0.8, ln .8= −0.22314 (3.10) Hence, for this choice of p 0 , the parameter α 0 has the value If an investigator is working with human data and if some period life table are available for a few past centuries, then it would be possible to use this data as a guide for the choice of the numerical value for the parameter p 0 .
The next step is to consider a rational for choosing a value for the parameter β 0 .As can be seen from equation 3.6, the expression (1 − exp ), decreases the value of the parameter α 0 by a fraction for each x > 0.Suppose, for example, that at age 3 years the equation (3.13) and β 0 has the value (3.14) so that the integral of the risk function has the value at x = 3.
Observe that this value is smaller than the value chosen for α 0 = 0.22314.
As can be seen from 3.7 and 3.11, From formula 3.6, it can be seen that H 0 (x) increases for any chosen value of the parameter β 0 , but will eventually converge to the limit α 0 .
In studies of human mortality in computer simulation experiments, ages from 1 to 100 may be under consideration.In such experiments, it has been observed that the value of H 0 (100) is the limit α 0 .For the illustrative numerical example under consideration this value is H 0 (100) = 0.22314 .
(3.17 From equation 2.21, it can be seen that if the survival function S 0 (x) were included in this product, then it may reduce the values of the Gompertz survival more than intended for ages x > 30.If an investigator, thinks this it is valid for the situation under consideration, then the survival function S 0 (x) would be included in the product in 2.21 for all ages x ≥ 1.
For example in such a formulation of a model on evolutionary genetics, the function S 0 (x) could represent a genotype such that individuals of this genotype are more at risk of dying than those of other genotypes so that this survival function would be considered be valid.
An alternative to formulating a survival function is presented in Mode and Busby .A brief history of risk functions just outlined is incomplete.In this connection, it is suggested that an interested reader have a look at the paper of Mode and Jacobson (1984) and the references cited therein.Briefly, the risk of death function used in that paper for early life was related to the probability density function of the Weibull distribution, which was much more complicated than the one in 3.5.In the early planning stages for this paper, a decision was made to use risks functions that have only a few parameters in order to minimize the number of parameters that must be assigned values for computer experimenters connected with evolutionary genetics.

A Brief History of Risk Functions That Have Been Proposed for Midlife
In his proposal for a risk function for the whole life, Thiele (1872) thei offered the risk function of the from for mid-life.As can be seen from this formula is nearly that of the probability density function of a normal distribution with mean c 1 and variance b 2 1 .If the formula below form is chosen at the risk function then it would be a multiple the probability density function of a normal distribution with mean c 1 and variance b 2 1 .Let . By definition the function Φ (x) is the distribution function of the standard normal distribution with mean µ = 0 and variance σ 2 = 1.It can be shown that the integral of the risk function in 4.2 is Therefore, the survival function S 1 (x) has the form From this observation, it follows from 4.4 and 4.5 that Because Φ (x) → 0 as x → −∞, it follows that S 1 (x) may not → 1 as x → 0. For example, from 4.4 and 4.5, it can be seen that Consequently, it seems desirable to search for a distribution such that the range of the random variable X representing age takes on only non-negative values x ≥ 0. Among the distributions with this property is the log-normal distribution, which is defined as follows.Suppose the random variable W has a normal distribution with mean µ and variance σ 2 .Then the random variable X = exp [W] , (4.9) then the random variable X takes on values x such that x > 0. From 4.9 it follows that ln X = W . (4.10) Thus, the random variable X is said to have a log-normal distribution.In what follows, it will be helpful to represent the random variable W in the form W = µ + σZ, where Z has a standard normal distribution with mean 0 and variance 1.
Given this representation of W, the distribution function of X may be easily derived.By definition is the distribution function of the random variable X.But, because the ln x is a monotone increasing function for x > 0, it follows that From this equation, it can be seen from 4.3 that the probability density function of the random variable X has the form From now on the risk function for mid-life will be chosen as µ 1 (x) = α 1 f (x) in 4.13.It can also be seen that from 4.12, it can also be shown that lim which is a desirable property.
To connect the parameters µ and σ with more useful statistical interpretations, it can be shown that the mode of the log-normal density is And if X 1 is a log-normal random variable, then its expectation is Equation 4.15 may be verified by finding the derivative of the log-normal density function in 4.13 and setting it to 0, and then solve that resulting equation for mode in 4.15.It would also be wise to check that m 1 is indeed the maximum of the density function in 4.13.As a first step in showing the µ 1 in 4.16 is the correct formula for the expectation in 4.16, it well be helpful to have a look at the formula for the moment generating function of normal random variable with mean µ and variance σ 2 .It is well known that the moment generating function of the normal distribution is which demonstrates the validity of 4.16.
The integral of the risk function µ 1 (x) = α 1 f (x), where f (x) is the log-normal density is Therefore, the survival function for mid-life has the form but the rate of convergence may be slow.
Given the formulas for the mode m 1 in 4.15 and that in 4.16 for µ 1 , the expectation of a random variable with a log-normal distribution, the next step in the formulation is to derive formulas for computing the parameters µ and σ as functions of m 1 and µ 1 , which in a computer simulation experiment would be assigned values based on some rational used in designing a computer experiment.From 4.15 and 4.16, it follows that Let v 1 and v 2 denote two 2 × 1 vectors defined as follows and Then consider the 2 × 2 matrix From 4.22 and 4, 23, it can be seen that It can be shown that the inverse of the matrix A is Thus, if µ , σ and α 1 are known, then function H 1 (x) in (4, 19) could be computed for chosen values of x > 0. Included in the scientific word processor used to write this paper are symbolic and numerical computation engines.The inverse matrix in 4.29 was calculated using the symbolic computation engine, and the numerical computation engine was used, was used to produce the the numerical computation engine was used to produce the numerical results presented in the rest of this section.
By way of an illustrative example, suppose that in some computer simulation experiment, it was assumed that ages of greatest risk of death in mid-life were x ∈ [15,35].If it assumed that the risk of death is greatest as age 25, then the value of the mode of the log-normal density function could be chosen as µ 1 = 25.From formula 4.32 in order to compute a value of σ > 0, the mean µ 1 of the log-normal distribution should be chosen as a number such that µ 1 > m 1 .Suppose, for example, the value chosen for m 1 = 22.Then, ln 25 = 3. 218 9 and ln 22 = 3. 091.Then ln 25 − ln 22 = 0.127 83, and Given these values of µ and σ, an investigator may wish to graph the log-normal density function in 4.13 for x ∈ [15,35].to provide insights into plausibility of the values the chosen values of µ and σ and their implications for predicting the potential results of the computer simulation experiments being considered.By assigning a value to the promoter α 1 , it would also be possible to compute the values of the survival function for any x > 0.
In the Mode and Sleeman 2000 , the risk function based on the log-normal distribution as indicated in section was suggested as a risk function for mid-life in section in chapter 13 section 2 of the book.But the account of this risk function is much more complete than in the book.Furthermore, one of the formulas presented in this chapter for finding µ and σ may contain an error, but in this section the formula was derived using a different method based on matrix theory and its derivation is more transparent.The review presented in this section is incomplete.For example, the risk function for mid-life suggested by Heligman and Pollard (1980) has not been displayed in this section It was not included, because of its complexity.But, if a reader is interested in an application of this risk function the paper by Mode and Jacobson (1984) may be consulted.

Makeham and Gompertz Risk Functions
The Makeham risk function is included in the formulation to accommodate the risks of accidents throughout the life span of an individual.It is assumed that this risk function is constant and is given by ,where α 2 is a positive constant.In this case, the integral of the risk function is for x ≥ 0. Therefore, for x ≥ 0, the survival function has the simple form (5. 3) It has been suggested that a useful trial value of α 2 , based on period studies of human mortality, is about 0.001.The two-parameter latent risk function, due to Gompertz (19-th Century) deals with risks of deaths at the older ages.Let α 3 and β 3 be positive parameters.Then, it will be assumed that for x ≥ 0 the latent risk function θ 3 (x) has the form (5.4) Observe that, as it should, this risk function increases as age x of an individual increases, and, by assumption, the risk of death increases exponentially with increasing age.The integral of this risk function has the form for x ≥ 0. Therefore, the latent survival function for this component is By applying a general but standard formula that a density is the risk function times the survival function, it can be seen that the probability density function of the Gompertz distribution has the form (5.7) for x ≥ 0. Although this distribution may be derived from intuitively appealing assumptions, it is more difficult to handle from a mathematical point of view than some other distributions that arise in probability and statistics.Nevertheless, because many advanced mathematical functions are now available to research workers in such software packages as MAPLE, MATHEMATICA, and MATLAB, an outline of the mathematics used in analyzing the Gompertz distribution seems appropriate.
Because the parameters α 3 and β 3 do not have obvious statistical interpretations, such as an expectation or variance, it is difficult to assign tentative values to them.Quite often, however, there is some feeling about the modal age of death for those who survive to old age.Let m 3 denote the mode of the Gompertz distribution.Then by using elementary calculus to find the maximum of the density f 3 (x), it can be shown that the equation supplies a connection among the parameters α 3 , β 3 and m 3 .In particular, if m 3 is assigned a value and β 3 is known, then α 3 is determined.But, to find a plausible value of β 3 , more input is needed.To this end, one may also have some idea about plausible values for σ 3 , the standard deviation of the age of death among those who survive to old age.Thus, it would be helpful to express σ 3 in terms of the parameters of the Gompertz distribution.
In this connection, it will be useful to consider the moment generating function of the Gompertz distribution, which is defined by for those values of s for which the integral converges.By using some advanced calculus, it can be shown that the moment generating function of the Gompertz distribution has the form where Γ(•) is the famous gamma function and s ∈ (−ϵ, ∞) with ϵ > 0 near zero.
Recall that the gamma function is defined by the integral (5.11) which converges for all z > 0. As is well known, this function corresponds with the factorial function on the positive integers.In fact, if z = n, a positive integer, then Γ(n) = (n − 1)!.In a word, this function fills the gaps among the integers, and, in particular, it can be shown that Γ( 1 2 ) = √ π.
If the random variable X 3 has a Gompertz distribution, then its expectation is (5.13) and the variance of the distribution is given by the well known formula (5.14) After considerable analysis, it can be shown that the exact formula for the expectation is (5.15)where C ≃ 0.57721 • •• is Euler's constant.Furthermore, the exact formula for the second moment is (5.16) The following two results Γ (1) (1) = −C (5.17) and (5.18) which were used in deriving the above formulas, are given in several books on special functions and advanced calculus.
The symbol Γ (k) (x) is, by definition, the k-th derivative of the Gamma function.
When α 3 > 0 is small, then exp[α 3 ] ≃ 1 and the above infinite series may be neglected.Thus, the approximations and (5.20) hold for small α 3 .Therefore, a formula for the approximate variance of the Gompertz distribution is . (5.21) Equivalently, is an approximate expression connecting the parameter β 3 with the standard deviations σ 3 .However, if an investigator is working with MAPLE, MATHEMATICA, MATLAB or other software packages, programs may be written to compute the variance σ 2 3 more accurately, using the above infinite series.Thus, the accuracy of these approximations may be assessed.Among the classic books that contain information on the above mathematics are those of Artin and Widder.

Introducing Randomness Into the Gompertz Survival Function
In a previous study using Monte Carlo simulation methods, it was found that an age structured stochastic process converged in distribution to a parametrized survival function, see Mode and Sleeman 2012 .In previous work the idea of a stochastic environment has a also been considered see Mode and Jacobson (1987) .In this section.some techniques for incorporating randomness into the Gompertz distribution will be considered.The Gompertz distribution is based on two parameters α and β.Consequently, in this section the focus of attention will be some suggestions for introducing randomness into these parameters.In equation 5.8 there is a formula connecting and parameters α and β as well as the mode m of the Gompertz distribution.
Among those who study and apply Monte Carlo simulation methods there is a well known equation.Let denote the distribution function of the random variable X.It will be assumed that the domain of the random variable X is a continuum of real numbers.Then where U is a random variable on the interval [0, 1].Whenever it is possible to solve this equation to express X as a function of U, then it is possible to compute realizations of the random variable X by computing realizations of the random variable U.
In order to view the mode of the Gompertz distribution as a random variable, it will be necessary to define positive numbers that represent ages.Let a and b be positive numbers representing ages such that a < b , and let V denote a random variable with a uniform distribution on the interval [a, b].Its probability density function is The distribution function of the random variable V is From the comments stated above, it follows that G (V) = U.Thus, By solving this equation for V, it can be seen that From this equation, it becomes transparent that realizations of the random variable V reduces to computing realizations of the random variable U.If U = 0,then V = a, and if U = 1,then V = b.From these observations, it follows that that realizations of the random variable V will always be in the interval [a, b].From 5.8, it follows that if the mode m is a random variable, then so is the parameter α.It is of interest to note that in this case β may be constant or a random variable.
When designing a computer simulation experiment, it is of interest to know the expectation and variance of random variables given in 6.6 and 6.8.It is well known that if the random variable U is uniformly distributed on the interval [0, 1], then its expectation and variance are given by the formulas and From these observations, it follows that the expectation and variance of the random variable V in 6.6 are and For example, if a = 50 and b = 70 as suggested above, then From these observations, it can be seen that the variation around the mean 60.0 will be relatively small in samples of realizations of the random variable V.In the next section, the concepts described in the preceding sections will be examined numerically and illustrated by using graphs.

A Method for Computing Numerical Versions of Parametric Survival Functions
In this numerical version of a survival function, the parameters values for early life, see section 3, were chosen as α 0 = 0.953, which is the probability of an infant surviving one year.The chosen value of the expectation of the exponential distribution was β 0 = 1.361.For midlife the parameter α 1 , which is the probability of surviving midlife, was chosen as α 1 = 1.The values of the mode and mean for the lognormal distribution were, see section 4 for details, were chosen as m 1 = 25 and µ 1 = 33.The value of the Makeham parameter was chosen as α 2 = 0.001.The value of the standard deviation σ 3 for the Gompertz distribution was chosen σ 3 = 13.255 and the value of the mode of the Gompertz was chosen as age m 3 = 75.A detailed account of these parameters is given in section 5.It is also shown in section 5 that the values of parameters determine the values of the Gompertz parameters α 3 and β 3 .The value of the greatest age was chosen as GA = 100.It should be emphasized that the values of the parameters stated above, which were chosen to produce a plausible numerical realization of a survival function for humans that may have long lives, are arbitrary.Any user of the ideas may chose any values of the parameters that are, in his mind, plausible.
Let S 0 (x) denote the function designed to adjust the Gompertz survival function for deaths in early life.Then, S 0 (0) = 1 and S 0 (x) was evaluated numerically for x = 1, 2, 3, • • •, 100.The function S 1 (x) was designed to adjust the Gompertz survival function for midlife deaths was also evaluated numerical for the integers x = 1, 2, 3, •••, 100.Similarly, the survival functions for the Makeham distribution S 2 (x) and Gompertz distribution S 3 (x) were also evaluated numerically for the positive integers listed above.In this section two additional survival functions will be considered and evaluated for x = 1, 2, 3, • • •, 100.Both these survival function satisfy the condition S 0123 (0) = 1 and S 023 (0) = 1.Observe that the survival function S 023 (x) differs from S 0123 (x) because the function S 1 (x) for midlife has been excluded.
If S (x) is the survival function for the random variable X with values in the interval 0 ≤ x ≤ 100, then it is well known that the expected value of X is As can seen from the graph of the survival function S 0123 , it appears that this function over adjusts the Gompertz so it is doubtful whether this function could be applied in a computer simulation experiment, unless the objective of the experiment was to project the evolution of a population with high mortality rates.The expectation of the life of an infant for this survival function is about 41 years,see 7.4 and is much smaller then that of the Gompertz survival function which is about 70 years, see 7.6.The reason why the function S 0123 (x) over adjust the Gompertz survival function is that the numerical value of this function at age 100 is less than one.If it were one, then its affects would be neutral and the numerical value of the Gompertz function would not be changed.
Figure 7.2 contains the graphs of the survival function S 023 and the Gompertz survival function.From these graphs it can be seen that the high rate of infant immortality my be easily seen by viewing the graph of the survival function S 023 in the early years of life.Given this survival function, the expectation of life is about 64 years, see 7.5.It seems more likely that this survival function would be more suitable for using in a computer simulation experiment than the function S 0123 when high mortality rates for infants are under consideration.In experiments not reported in this paper, an attempt was made to use the assumptions that the function S 0 (x) was less than one only for the ages x = 1, 2, • • •, 10 and S 0 (x) = 1 for ages x = 11, 12, • • •, 100.Similarly, the function for mid-life S 1 (x) was less than one only for the ages x = 15, 16, • • •, 30, representing mid-life.For all ages not in the mid-life range, it was assumed that S 1 (x) = 1.According to theory, when these functions were one, the Gompertz function would not be changed in the products 7.1 and 7.2.However, when these products were used at some ages x = 1, 2, • • •, 100, the Gompertz function was not monotone decreasing as the age x increased.Consequently, the adjusted Gompertz function was no longer useful as a survival function for all ages.A decision was made, therefore, to abandon the approach just described.
In previous work, Mode and Busby 1984, as mentioned in a previous section, survival function for infancy, mid-life and older ages were spliced together.When the parameters of this mode were estimated from a life table is was found that there were small bumps at the splice points, but the over all the fit to data was very good.Consequently, this splicing method remains an alternative to the methods described in this paper based on the numerical evaluation of all functions in a product of functions for all ages x = 0, 1, 2, • • •, 100.But, it should be mentioned that if all the parameters of the were assigned plausible numercal values not based on data, the bumps in the splice points may be very noticeable.
In section 6, a method was developed to treat the mode of the Gompertz as a random variable.In an illustrative example, the bounds of variation in the mode of the Gompertz function were 50 and 70.Because the computer implementation of these ideas would entail major a change to existing software, the ideas just mentioned will not be implemented in the software used in this paper.
Natural disasters such as earthquakes, floods and volcanic eruptions are often the cause of deaths in human and animal population in the regions of such disasters.Such disasters are difficult to predict, but it would be interest to design and implement some software based on the idea of random waiting times among disasters.An example of such random waiting times are those among earthquakes in some geological area area.For example in discrete time formulations, it may be assumed that the waiting times T among earthquakes are random variables in some geological area that follow a geometric distribution with the probability density function where 0 < p < 1 and t = 1, 2, • • •.
In any Monte Carlo simulation experiment, it would be straight forward to simulate a sequence T 1 , T 2 , • • • of the random variables that are realizations of the random variable T.In such an experiment, it would be required that some methods for simulating deaths would be needed.A very simple way of simulating such deaths across all ages is to assume that the parameter α 2 in the Makeham distribution is a random variable.To illustrate these ideas, suppose the survival function under consideration is is the Makeham survival function.The value of this parameter in the formulation under consideration was chosen as α 2 = 0.001.If, for example, if this value is increased to α 2 = 0.01, then the survival function in 7.12 would be significantly decreased for all ages x = 0, 1, 2, • • •, 100, and in any in any Monte Carlo simulation experiment the simulated number of death using the survival function in 7.11, would be much greater than if the parameter value α 2 = 0.001 were used.In the next section, an algorithm for simulating deaths across all ages will be described, and will be used if all parameters are constant and not random variables.

Algorithms for Projecting an Age Structured Population With Stochastic Laws of Evolution Forward in Time
The objective of this section is to define as set of algorithms to simulate the stochastic evolution of an age structured population forward in time.It will be assumed that the greatest age any individual may live is 50 years.At time t let X (x, t) denote the number of females in a population for age x = 0, 1, 2, • • •, 50 for t = 0, 1, 2, • • •, N, where N > 1 is the number of years considered in a projection experiment.Let X (0, t) denote the number of female infants in the population at time t.Similarly, let Y (0, t) denote the number of male infants in the population at time t , and let Y (x, t) denote the number of males in the population at time t for ages x = 1, 2, • • •, 50.
Let S f (x) and S m (x) denote, respectively, the survival functions for females and males.By definition S f (0) = S m (0) = 1.Let p f (x) denote the conditional probability that a female of age x − 1, survives to age x for x = 1, 2, • • •, 50.Then, it follows that Observe that p f (1) is the conditional probability that an infant female survives to age 1.Similarly, let p m (x) denote the conditional probability that a male of age x − 1 survives to age x for x = 1, 2, • • •, 50.The evolution of the population will be viewed as stochastic, because it will be assumed that the number of females X (x, t) is a random variable who's realizations may be simulated.
A random variable X is said to have a binomial distribution with parameters N and p if its pdf has the formula In this formula N is a positive integer and p is a probability such that (0 < p < 1) In what follows, it will be convenient to use a kind of shorthand such as to indicate that a random variable X has a binomial distribution with parameters N and p.
The next step in the formulation of the evolutionary model is to take into account the births of infants and their entry into an evolving population.From now on the focus of attention will be the evolution of a population of wild chimpanzees.
Suppose that the ages female chimps may give birth to an infant are x = 15, 16, • • •.40.There are several distributions depending on one or more parameters that may be viewed as candidates for the number of infants a fertile female chimp may give birth during every year of her fertile ages.Let the random variable Y (x) denote the number of offspring a female chimp may have during any year for ages x = 1, 2, • • •, 40.For the sake of simplicity, it will be assumed that the random variable Y (x) has a Poisson distribution with probability density function It is well known that the expectation and variance of B is λ.Let Z be a standard normal random variable with expectation 0 and variance 1.Then it is also well known that the distribution of the random variable when the parameter λ is large, Z has a standard normal distribution with mean 0 and variance 1.From 8, 7 it can be seen that Observe that if λ (x) = 1 for x = 15, 16, • • •, 40.then λ = 26, which is sufficiently large to use the approximation 8.7.The number B in 8, 8 may not be positive.Thus, B will be chosen as Therefore, the distribution of the random variable B is the absolute normal distribution with parameter λ.In general, the number B in 8.9 is not a positive integer.To remedy this situation, for any positive real number X, let [X] denote the greatest integer in X.Thus, if B is chosen as (8.10) then B is non-negative integer.
The next step in the formulation is to define an algorithm for computing the number female and male infant chimps that are added to the population at any time t.Let p f denote the probability that an infant chimp is female, and let p m denote the probability an infant is male.Then, p m = 1 − p f .Let the random variable X (0, t) the number of female infant chimps added to the population at time t, and let B (t) the total number of infants that will become members of the population at time t.Then, it will be assumed that .11)From this result it follows that the number of male infant chimps added to the population at time t is At time t, the number of females in the population by age may be represented by the array denote the array of females by age in the population who are the survivors of the population depicted in 8.13.Then, the survivors in 8.14 are simulated using the formula for x = 0, 1, • • •, 49.Note if there are individuals of age 50 in the population, then they will no longer be viewed as members of an evolving population.
Observe that the age structured population represented by the array in 8.14 does not contain female infants.Therefore, the number X (0, t + 1) of female infants at time t + 1 could be simulated using the algorithm described above and added to the array 8.14 to simulate the array By using the algorithm just described, an array analogous to that in 8.16 could also be simulated for males.In this paper, Monte Carlo simulation methods will be used to simulate a realizations of the stochastic process under consideration, using the algorithms described above.The numerical values of the parameters of the model used in this experiment will be given in the next section.
At any time during a population projection, it is of interest to display an estimates of total population size and the age distribution.At time t, the total size of the female population is Therefore, at any time t, the age distribution is If both females and males are being considered in a population projection experiment, then the same formulas can be used to estimate the total population size for males as well as their age distribution.
When dealing with a stochastic processes evolving in time, it is often of interest to estimate the value of some random variable at time t + 1, given the value of the random variable at time t.For the stochastic processes under consideration, it is well known that the minimum mean square estimate X (x, t + 1) of the random variable X (x, t + 1) is the conditional expectation Observe that this is the expectation of a random variable with a binomial distribution as indicated in 8.3.Estimates of the random function Y (x.t + 1) for males could also be computed using a formula of the type given in 8.19.
The formula in 8.19 may also be written in the form which provides a recursive formula for computing estimates X (x, t) of the random function X (x, t) for t = 1, 2, • • •, given X (x, 0) for all ages x = 0, 1, 2 • ••, 50.A similar recursive equation may be used to compute Y (x, t) as an estimate of the random function Y (x, t) for t = 1, 2, • • • and x = 0, 1, 2, • • •, 50.At each t in the recursive procedure, the function X (0, t) will be assigned the value X (0, t) = λ,the expected total number of births, see 8.6.
The recursive equations for the females and males just described will be referred to as the deterministic model embedded in the stochastic process.For the case of the embedded deterministic model only one trajectory of the process can be computed, Whereas, for the stochastic model in a computer simulation experiment, each trajectory of the process will differ from all other trajectories in a sample of R > 1 realizations of the process.It should also be mentioned that when the embedded deterministic model is used, it will be possible to compute age distributions as well as other sets of values that may be of interest to an investigator.
During the preliminary testing of the ideas underlying stochastic process, it was observed that if it is assumed that λ, the expected total number of offspring born in any year was constant, then the number of births in any realization of the process was nearly constant, i.e., it was not random.To correct this flaw, it was decided to let λ be a random variable defined as follows.Let a denote a positive number less than b.Then it was assumed that λ was a random variable defined by the equation where U is a uniform random variable on the interval [0, 1].At this point, it may be helpful to consult section 6 that contains details on the properties of random variables defined in equation 8.21.For example, the expectation of the random variable λ is In computer simulation experiments the numbers a and b are usually chosen such that the expression on the right is a preassigned positive number.

Assigned Values of Parameters Used in an Illustrative Monte Carlo Simulation Experiment
The parameter values chosen in the simulation experiment considered in this paper were designed to follow the short term evolution a population of animals living in the wild with short lives.It was assumed that the greatest age of any individual in the population as it evolved was 50 years.Differences in the survival on males and female were taken into account by selecting two survival functions with different parameters values.The survival functions for both sexes were chosen using the early life survival function in section 3 and Makeham Gompertz survival functions described in sections 3 and 5.
The survival function, S 023F was chosen for both females and males in the experiment under consideration.The reason for choosing this type of survival function for both females and males is that there is no information from chimp populations in the wild that there is a midlife effect with regarding mortality.as there is in humans, particularly for males in midlife.
.Consequently the midlife component S 1 was omitted.
Observe that to get a numerical version of the survival function S 0 for females, it is necessary to assign values for two parameters in S 0 , the fraction that survive infancy and a mean of an exponential distribution.The fraction that survived infancy was chosen as α 0 = 0, 80 and the mean of the exponential distribution was assigned the value β 0 = 030030.The parameter for the Makeham distribution was assigned the value α 2 = 0.02.Finally to get a numercal version of the Gompertz survival function numerical values must be assigned for the standard deviation σ G and the mode m G .
The numbers chosen for these two parameters were σ G = 6 and m G = 33.All the parameters for the male survival function were the same as those for the females except that the mode of the Gompertz distribution for males was chosen as m G = 30.Numerical values for the female S 023F and male the survival functions S 023M were computed by using the procedures set forth in section 7.
Very little about survivorship in populations of wild Chimpanzees is known.Consequently every numerical version of a survival function must be checked for plausibility.By using the formulas for the expection of life given in section 7, the numerical values for the expectation of life, expressed in terms of years, for females was The difference in these two numbers is due to different modes assigned to the Gompertz survival function.For females the mode of the Gompertz distribution was chosen as m G = 33 and that for males was chosen as m G = 30.These two expectations of life seem to be plausible for wild Chimpanzee populations.
It is often informative to include graph forms of the two survival functions under consideration.Presented in Figure 9.11 is a graph of the survival functions for females and males.
In the simulation experiment under consideration, a realization of the embedded deterministic model as well a sample of realizatioins of the stochastic model were computed.The parameter assignments for the deterministic model were as follows.The number of years of evolution of both the female and male populations was set as 200.The expected number of the total infants born in any year was assigned the value λ = 20.The initial number of females and males in the population was chosen as 300.The probability that an offspring was female was chosen as p f = 100/205, and the probability an offspring was male was the complementary probability p m = 1 − p f .
For the case of the stochastic model, the number of years of evolution of the population was 200, and the number of replications of the experiment was chosen as 100.The initial number for both females and males was 300 and was the same as those used in the experiment with the deterministic model.In the experiment with the stochastic model λ was a random variable, and the numbers in equation 8.21 were chosen as a = 15 and b = 25.It is interesting to note that the expected value of the random variable λ in this case is which is the same value for λ used in the experiment with the deterministic model.An important part of the procedure in setting up a computer simulation experiment is that of specifying the initial age structure for the female and male populations.Let X (0, 0) denote the initial age structure of the female population and let Y (0, 0) denote the initial age structure for the male population at time t = 0. Observe that X (0, 0) and Y (0, 0) are vectors with 1 row and 51 columns, and.as indicated above], the survival function used studying the 200 years of evolution of the female population was S 023F and that for males it was S 023M .The initial vector for the female population was chosen as and similarly that for the males was chosen as Y (0, 0) = 300 × S 023M (9.5)These assignments will work well in experiments with the deterministic model, but when using the stochastic model the vectors of the female and male populations must have values in the set S = (x | x = 0, 1, 2, 3, • • •).This set is often referred to as the set of non-negative integers.Thus if the initial vectors for the age structure of the female and male populations are as X (0, 0) = [X (0, 0)] (9.6) and Y (0, 0) = [Y (0, 0)] , (9.7) where the symbols [X (0, 0)] and [Y (0, 0)] indicate that the greatest integer of each element in the vectors X (0, 0) and Y (0.0) have been selected.The operation defined in equations 9.7 and 9.8 ensures that in the Monte Carlo simulation experiment the elements of all vector belong to the set S of non-negative integers.The software used in the Monte Carlo simulation experiment makes it necessary that the elements of all vectors are in the set S.

Analysis and Interpretation of Monte Carlo Simulated Data
The results presented in this section were based data that was simulated by what is known by Monte Carlo Methods.A review of Monte Carlo Simulation Methods as they apply to mutation and selection formulated in Wright-Fisher models of evolution was published by Mode and Gallop (2008).In this paper, the formula used to compute random numbers in the experiments described in this section is also presented.
By definition the age structure of a population for females and males is the number of individuals in each age classification: 0, 1, 2, • • •, 51.As expected in the experiment with the embedded deteministic model, the age structure of the population converged to a constant for each age classification for both females and males.By a little mode than 50 years into the population projection experiment, convergence to constant age structure had occurred for both sexes.Various terms are used to describe such a population.Among these terms is to say that the projected population has converged to an eqilibrium.Another term that is used that the population has converged to a stable population.
It is of interest to display a graph of the age structure that the deteministic process converged to for females and males after about 50 years into the projection of 200 years as presented in Figure 10.1.The method chosen to statistically summarize the simulated data was to order each column of the array from the smallest number to the largest.In this procedure the array is expressed in terms of the order statistics for each of the 51 ages considered in the experiment.Row 1 of this array contains the minimum age Min for each age classification.And similarly, row 20, 000 contains the maximum age Max for each age classification.Let Q 25 denote the 0.25 quantile, which is the fraction of all columns such that the inequality is valid of the ordered data array.Q 50 is also known as the median of the data.It is easy to see that the quantile Q 75 may also be defined.For the case under consideration Max is defined as row 20, 000 of the ordered data array.It is clear that the notation just described could be extended to the case N, some large number rows and that software could be written to take into account extended the case just outlined Presented in Figure 10.2 are the graphs of the Min, Max and Quantiles for females.Another statistic that was included in the simulated data was total population size for each year of the projection of 200 years.This number was computed for each year of the projection by summing over all ages and sexes.In the experiment with the stochastic model, a sample of total population size was computed for each year the projection, which was replicated 100 times.Presented in Table 10.1 are the summary statistics for this experiment.
From this table it can be seen that there is also an outlier 9744 for total population size.It is also interesting to note that 0.75 of the realizations of the process are less than or equal to 360.It is of compelling interest to attempt to find an explanation as to why there are outliers in the simulated data.It is well known that for the type of age structured model under consideration that the initial age structure of the population has a profound effect on the subsequent evolution of the population that can seen in the experiments with the deterministic and stochastic models.For the case of the deterministic model the expectation of the total number of offspring per year was assigned as constant number.Given this number and a numerical version of the survival function, the population will evolve to a fixed age structure.If the initial age structure differs from the fixed age structure, then unexpected results may be seen in the simulated data such as outliers that may appear in the first years of a projection.The initial age structure used in the experiment with the deterministic model was also used in the experiment with the stochastic model.In the experiment with the stochastic model, outliers were also observed during the initial years of an experiment.
It is of interest to compare, the evolution of total population size that was observed when using the deterministic model with that of the stochastic model.In the early years of the projection using the determinist model, total population size was in the nine thousands but in about fifty years it had converged to 358.912, which is close to the median for the stochastic model displayed in table 10.1.Presented in Figure 10.4 is a comparison of total population size with the Q50 quantile for the stochastic model.As can be seen from Figure 10.4, the trajectory for total population size is close to that of the Q50, the median trajectory for females.
In the experiments with the deterministic and stochastic models, a survival function was computed for females and males in each year of the projection.In the experiment with the stochastic model each projection of 200 years was replicated 100 times and the mean survival function were computed for both females and males In Figure 10.5 a comparison plot of the deterministic survival function in year 100 of the projection with the mean survival function for females computed from the 100 realizations of the stochastic process.Similar graphs for the male population were observed, but they have not been displayed here because they were similar to the graphs for the female population.
When attempting to forecast whether a small population of chimpanzee or that of another species would still exits 200 years into the future, the presence of outliers in a computer simulation experiment is troublesome, because it would not be possible for a small population of chimps to grow to a total population size of over 9, 000 in a few years.Such outliers occur when the initial age structure for females and males differs from that determined by the expected number of offspring each year and the conditional probabilities that individuals of each age survive another year.In the experiment with the deterministic model, the convergence to a constant age structure for females and males occurred between 50 and 100 years into the projection.Given this information, a decision was made to do an experiment with the stochastic model using a slight modification of the age structures for females and males that arose as limits in the experiment with the deterministic model.Presented in Table 10.2 are the summary statistics computed from the simulated data for total Unlike experiment 1 on total population size presented in Table 10.1, where the Max statistic is 9744, the Max statistic in Table 10.2 is only 406 and would not be classified as an outlier when compared with the other statistics in Table 10.2.This second experiment provides evidence that when the initial age structures for females and males as chosen as the greatest integers in the stationary age structure of the deteministic are used as the initial age structures in experiments with the stochastic model, then the stochastic model performs well in the sense that no outliers appeared in the simulated data on total population size.
It is also of interest to display a graph of the of selected quantiles as well as the Min and Max age structures.Figure 10.6 contains the graph of interest for the female population.
The information conveyed in Figure 10.6 may be viewed as a forecast as to whether an existing population of chimps will be the ancestors of a population of chimps in 200 years into the future.It is not possible to state with certainty that an ancestral population will exist 200 years into the future.But, the statistical information in Figure 10.6 will give conservationists ideas about the uncertainties of the existence of a future population.As can be seen from Figure 10.6 the age structure of the quantiles Q 25 , Q 50 and Q 75 are zero only at the higher ages as one would expect.The age structure of the Max is also well above zero except for the higher ages.It is very interesting to note that age structure of the Min is zero beyond age 15.This result suggest that in some realizations of the process the population became extinct.But, on the other hand, the statistical information in Table 10.2 on the total population size provides evidence that the population would not become extinct during 200 years of evolution.The graph for males was very similar to that for females, and has, therefore, been omitted.
It is well known that the class of stochastic models under consideration will converge to what is called a quasi-stationary distribution, given that extinction does not occur.In the experiment just discussed, it is not clear as to whether the data generated by the stochastic mode did converge to the quasi-stationary distribution, because 200 years of evolution is too short a time for convergence to occur.A knowledge of the existence of a quasi-stationary distribution will, however, provide a structure for dealing with the uncertainty that arises in Monte Carlo simulation experiments with stochastic models.
The results just described provide convincing evidence that after all the parameter values have been assigned numerical values in the deterministic and stochastic models, then the first step of an experiment would be that of using the deterministic model to get estimates of the initial age structures for the female and male populations for the stochastic model.In the experiments reported in this section, a few seconds of running time of the deterministic were sufficient to compute 200 years of evolution.The next step in the experiment would that of using the adjusted age structures that were estimated in the experiment with the determinist model as initial age structures in an experiment with the stochastic model.Even though outliers may be observed in the experiment with the deterministic model, it seems very unlikely that outliers would be observed in an experiment with the stochastic model as suggested by observed results of experiment 2. The running time to the completion of the experiment 2 with the stochastic model was in the range of 2 to 4 minutes.These relatively short running times make it feasible to do several exploratory experiments before selecting an experiment with an intention of publishing the results.

For
example, if the mode of the Gompertz distribution is m = 60 , b = 70 and a = 50, then b − a = 20, then the above equation takes the form V = 50 + 20U.(6.7)In this case, in a Monte Carlo simulation experiment, the mode of the Gompertz distribution would be somewhere in the interval [50, 70] for any realization of the process.Similarly, if the standard deviation σ of the Gompertz distribution is viewed as a random variable such that all realizations of σ belong to the interval [c, d], where c and d are positive numbers such that c < d, then σ = c + (d − c) U .(6.8) known as the expectation of the life span for any newly born infant in the population.For the case of the survival function S 0123 this expectation has the value E [X | S 0123 ] = 40.95676571 .(7.4)The expectation of S 023 is E [X | S 023 ] = 64.084919 .(7.5)For the case of the Gompertz survival function, this expection has the value E [X | S 3 ] = 69.592489 .(7.6) Presented in Figure 7.1 is a graph comparing the Gompertz survival function with S 0123 for the ages x = 0, 1, 2, • • •, 100.

Figure 7
Figure 7.1.A Comparison of the Graphs of the Survival Function S 0123 and the Gompertz Survival Function For the case of a mode of 50, the expectation of life would be E [X | S 3 , m 3 = 50] = 44.968, (7.9) and the case m 3 = 70, this expectation is E [X | S 3 , m 3 = 70] = 64.621(7.10)Thus, in any Monte Carlo simulation experiment based on these bounds, the expectation of the Gompertz distribution would lie in the interval [44.968, 64.621], which shows that in such a population the number of deaths would vary sig-nificantly among years.
4)with parameter λ (x) > 0 for each age fertile age x.Let the random variable B denote the total number infants the females of a chimp population have during any year of a population projection.Then, if it is assumed that the random variables Y (x) for x = 15, 16, • • •.40 are independent, then the random variable

Figure 9
Figure 9.1.A Graph of the Female and Male Survival Functions

Figure 10
Figure 10.2.Graphs of the MIN, MAX and Quantiles for Females

Figure 10
Figure 10.4.A Comparison of Age Structure According to the Deterministic Model in Year 100 with the Median for the Stochastic Model

Figure 10
Figure 10.5.A Comparison of the Computed Survival Functions for the Deterministic and Stochastic Models

Figure 10
Figure 10.6.Graph of Selected Quantiles and the Min and Max for Females

Table 10
.1.Summary Statistics for Total Population Size 

Table 10 .
2. Summary Statistics for Total Population Size for Experiment 2 with the Stochastic Model