Roff D.A. Modeling Evolution: An Introduction to Numerical Methods

Подождите немного. Документ загружается.

2.7.3 Plotting the ﬁtness function

To plot r, our ﬁtness measure, as a function of x we must solve equation (2.35) or

(2.36) for each value of x. We shall ﬁrst examine methods to estimate r without

ourselves doing the integration (i.e., using equation (2.35)) and then methods

using equation (2.36).

2.7.3.1 Using the program to do the integration

R CODE (Figure 2.5):

We ﬁrst consider numerical integration to estimate r. The strategy is to iterate

over a range of x and for each value calculate r.

1. INTEGRAND calculatesða

þ b

xÞe

ða

þb

xþrÞage

2. INTEGRAL calls the R function integrate to obtain the integral, passing to it

INTEGRAND.

3. The function RCALC uses the R function uniroot to ﬁnd the value of r which

satisﬁes the characteristic equation (2.35) for a given x.

4. The main program creates a single column matrix of x values and uses the R

function apply to use RCALC on each x value (row of matrix x) to get the

requisite value of r, stored in the vector r.

5. The vectors x and r are then used in plot to display the relationship between r

(ﬁtness) and x (body size). In sequence, for a given x, RCALC calls uniroot

which calls INTEGRAL which in turn calls INTEGRAND.

rm(list¼ls()) # Remove all objects from memory

# Function to output function to be integrated

INTEGRAND <- function(age,x,r)

{

Af <-0;Bf<- 4*4 ; As <- 1 ;Bs <- 0.5 # parameter values

return ((AfþBf*x)*exp((AsþBs*xþr)*age)) # return function

}

# Function to integrate characteristic equation and return 1-its

value

INTEGRAL <- function(r,x)

{ 1-integrate(INTEGRAND,1,Inf,x,r)$value} # 1-Characteristic

equation

# Function to find r given x

RCALC <- function(x){uniroot(INTEGRAL, interval¼c(1e-7,10),

x)$root}

# MAIN PROGRAM

x <- matrix(seq(0.5,3, length¼100)) # x ranging from 0.5 to 3

length¼100

r <-

apply(x,1,RCALC) # Calculate r for given x

plot(x,r,type¼’l’, xlab¼’Body size, x’, ylab¼’Fitness,

r’,las¼1,lwd¼4)

OUTPUT: Plot shown in Figure 2.5.

88 MODELING EVOLUTION

MATLAB CODE: See Section 2.18.10.

2.7.3.2 User supplied solution to the integral

Using the integrated function shown in equation (2.36) is simpler. However,

because many functions cannot be integrated analytically, it is less general.

R CODE (Figure 2.5):

rm(list¼ls()) # Remove all objects from memory

# Function to evaluate equation (2.36)

FUNC <- function(r,x)

{

Af <-0;Bf<- 4*4 ; As <- 1 ;Bs <- 0.5 # Parameter values

S <- exp(-(rþAsþBs*x))*(AfþBf*x)/(AsþBs*xþr)# RHS of equation

return(1-S) # Subtract 1

}

# Function to find r given x using uniroot

RCALC <- function(x){uniroot( FUNC, interval¼c(1e-07,10),x)$root}

# MAIN PROGRAM – same as previous

x <- matrix(seq(0.5,3, length¼100)) # X ranging from 0.5 to 3

r <- apply(x,1,RCALC) # Calculate r for given x and store

plot(x, r, type¼’l’, xlab¼“Size,x”, ylab¼“Fitness, r”,las¼1,lwd¼4)

MATLAB CODE: See Section 2.18.10.

2.7.4 Finding the maximum using the calculus

Our ﬁtness measure is no longer on one side of the equation, and the equation

cannot be simpliﬁed to make this so. The equation can be differentiated using

implicit differentiation (see Appendix 2). For convenience we ﬁrst take logs to

convert the function into one that is additive:

Body size, x

0.5

0.35

0.40

0.45

0.50

0.55

1.0 1.5 2.0 2.5 3.0

Fitness, r

Figure 2.5 Scenario 5: Fitness versus body size.

FISHERIAN OPTIMALITY MODELS 89

0 ¼ ln ða

þ b

xÞða

þ b

x þ rÞlnða

þ b

x þ rÞ

¼ T1  T2  T3

ð2:37Þ

Now taking each term separately we have

T1 :

þ b

T2 :

¼ b

T3 :

þ b

x þ r

þ b

þ r

ð2:38Þ

Gathering the terms together and rearranging we arrive at

1 þ

þ b

x þ r



þ b

 b



þ b

x þ r

ð2:39Þ

Now

¼ 0 when the right-hand side equals zero, provided the term in parenth-

eses on the left-hand side does not also equal zero, which will not generally be the

case. Thus we can rearrange the right-hand side to make r a function of the

parameters and x:



ða

þ b



 b



 b



 a

ð2:40Þ

where r* is the maximum value of r, obtained at x*. To ﬁnd x* we substitute

equation (2.40) into equation (2.37) to obtain:

0 ¼ ln ða

þ b



Þða

þ b



þ r



Þ ln ða

þ b



þ r



0 ¼ ln ða

þ b



Þ a

þ b



ða

þ b



 b



 b



 a

ln a

þ b



ða

þ b



 b



 b



 a

ð2:41Þ

which can be solved numerically.

R CODE:

I here illustrate three functions that can be used to locate the required value of x .

2.7.4.1 Use of uniroot

The ﬁrst is uniroot. Note that the limits are set fairly close to the required value.

If the limits are set too far apart (e.g., uniroot(f¼RFUNC,interval¼c

(1.2,3))$root ) the function may fail, R returning the error message:

Error in uniroot(f ¼ RFUNC, interval ¼ c(1.2, 3)) :

f.upper ¼ f(upper) is NA

In addition: Warning message:

In log(As þ Bs * x þ r) : NaNs produced

90 MODELING EVOLUTION

The above failure emphasizes the importance of making a preliminary plot to

determine the approximate value of x at which r is a maximum:

rm(list¼ls()) # Remove all objects from memory

# Function to calculate value of equation (2.41) for a given value of x

RFUNC <- function(x)

{

Af <- 0;Bf<-4*4;As<-1;Bs<-0.5 # Set parameter va-

lues

r <- Bs*(AfþBf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As # r from eqn (2.40)

return(log(AfþBf*x)-(AsþBs*xþr)-log(AsþBs*xþr))

}

# MAIN PROGRAM

uniroot(f¼RFUNC,interval¼c(1.2,1.8))$root

OUTPUT:

[1] 1.389974

2.7.4.2 Use of nlm

An alternate approach is to use nlm by taking the absolute value in FUNC, in which

case the minimum must be zero.

rm(list¼ls()) # Remove all objects from memory

# Function to calculate value of equation (2.41) for a given value of x

RFUNC <- function(x)

{

Af<-0;Bf<-4*4;As<-1;Bs<-0.5 # Set parameter values

r <- Bs*(AfþBf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As#rfromeqn(2.40)

return(abs(log(AfþBf*x)-(AsþBs*xþr)-log(AsþBs*xþr)))

}

# MAIN PROGRAM

nlm(RFUNC, p¼1.2)$estimate

Because x is not constrained, this method is not very satisfactory and although the

correct answer is found warning messages are generated.

OUTPUT:

[1] 1.389943

Warning messages:

1: In log(As þ Bs * x þ r) : NaNs produced

2: In nlm(RFUNC, p ¼1.2) : NA/Inf replaced by maximum positive value

2.7.4.3 Use of optimize

The third method is to use the routine optimize, which can ﬁnd either a

minimum or maximum and allows one to enter upper and lower limits for the

search. As in the case of nlm we use the absolute value of the function.

FISHERIAN OPTIMALITY MODELS 91

rm(list¼ls()) # Remove all objects from memory

# Function to calculate value of equation (2.41) for a given value of x

RFUNC <- function(x)

{

Af<-0;Bf<-4*4;As<-1;Bs<-0.5 # Set parameter values

r <-Bs*(AfþBf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As # r fromeqn (2.40)

return(abs(log(AfþBf*x)-(AsþBs*xþr)-log(AsþBs*xþr)))

}

optimize(f ¼RFUNC, interval ¼c(1.2,1.8),m aximu m¼FALSE)$minimum

OUTPUT:

[1] 1.389956

As with nlm, if the limits are too broad a warning message may be generated:

optimize(f ¼ RFUNC, interval ¼ c(1.2,3),maximum ¼ FALSE)$minimum

OUTPUT:

[1] 1.389951

Warning messages:

1: In log(As þ Bs * x þ r) : NaNs produced

2: In optimize(f ¼ RFUNC, interval ¼ c(1.3, 3), maximum ¼ FALSE) :

NA/Inf replaced by maximum positive value

Note that the four outputs are close but not exactly the same.

MATLAB Code: See Section 2.18.11.

2.7.5 Finding the maximum using a numerical approach

We shall consider two approaches: ﬁrst, using the solved integral equation (2.36)

and second, using numerical integration to solve the original model (equation

[2.35]). The latter method is the most general but will also be the most time

consuming. Even if the integral can be solved, it is good practice to use the latter

method as a check on the correctness of the integration.

2.7.5.1 Using the integrated function (equation [2.36])

The main program uses the R function optimize, giving it RCALC which uses

uniroot to ﬁnd the root of the function speciﬁed in FUNC.

rm(list¼ls()) # Remove all objects from memory

# Function to get integral value

FUNC <- function(r,x)

{

Af <-0;Bf<- 4*4 ; As <- 1 ;Bs <- 0.5 # parameter values

S <- exp(-(rþAsþBs*x))*(AfþBf*x)/(AsþBs*xþr) # Function value

return(1-S)

}

92 MODELING EVOLUTION

# Function to find r given x

RCALC <- function(x){uniroot( FUNC, interval¼c(1e-07,10),x)

$root}

# MAIN PROGRAM

# Use optimize which allows us to constrain the search

# Tell optimize that we want the maximum

optimize(f ¼ RCALC, interval ¼ c(.5,3),maximum ¼ TRUE)$maximum

OUTPUT:

[1] 1.389946

2.7.5.2 Using numerical integration of the function.

The main program uses optimize, giving it RCALC which uses uniroot to ﬁnd

the root of the function speciﬁed in FUNC, which itself calls INTEGRAND to obtain

the integral.

rm(list¼ls()) # Remove all objects from memory

# Function to do numerical integration of eqn (2.34)

INTEGRAND <- function(age,x,r)

{

Af <-0;Bf<- 4*4 ; As <- 1 ;Bs <- 0.5 # parameter values

return ((AfþBf*x)*exp((AsþBs*xþ r)*age)) # return function

}

# Function to integrate characteristic equation

FUNC <- function(r,x){1-integrate(INTEGRAND,1,Inf,x,r)$value}

# Function to find r given x

RCALC <- function(x){uniroot( FUNC, interval¼c(1e-0 7,10),x)$root}

# MAIN PROGRAM

optimize(f¼RCALC, interval ¼c(1.2,1.8),maximum ¼TRUE)$maximum

OUTPUT:

[1] 1.389934

MATLAB CODE: See Section 2.18.12.

2.8 Scenario 6: Stochastic variation in parameters

Within any generation there is likely to be stochastic variation in parameter

values. The scenario considered here is that in which one or more parameters

vary within generations according to some probability distribution with the

same distribution among generations. For variation in a single parameter, y,

ﬁtness for a given value of x is given by

PðyÞf ðy; xÞdy,whereP(y)isthe

probability density function and f( y,x) is the ﬁtness function for a given y

and x. A possible model would be a normal density function, so that the

probability of a value of y is given by PðyÞ¼

ﬃﬃﬃﬃ



ym



,wherem(y)isthe

mean parameter value and s

is its variance. Such a function could be

FISHERIAN OPTIMALITY MODELS 93

problematic for traits such as survival that can only vary between zero and

one, because the probability function must either be truncated or values less

than zero set to zero and values greater than one set to one. The problem

with adopting such a model is that the probability distributi on could show

unlikely rises at zero and one. Considerable care should be taken in the

selection of a probability function and in general it is ad visable to use several

functions or one that can take a wide variety of shapes to ensure that the

results are not uniquely dependent on the function chosen. The analysis of

temporal variation in a single parameter is dealt with in the following two

scenarios, the methods discussed in these scenarios are readily adapted to the

present. Here I shall consider the somewhat more complex case, but that

which is most likely to be appropriate for stochastic va riation, in which

several parameters vary. If there is variation in several parameters we need

to take into account possible correlations between parameter values. This

becomes a particular concern when the parameters are themselves traits,

where genetic and phenotypic correlations may be highly likely. For two

parameters the general form is a simple extension of the one trait case,

namely

Pðy

; y

Þf ðy

; y

; xÞdy

. Extension to more parameters is obvious

but solving such equations can prove difﬁcult. A brute force approach will

always work but could be time consuming if the model is very detailed. To

illustrate the general approach using brute force I shall consider the ca se of

variation in two uncorrelated parameters. The case considered here is con-

ceptually the same as in Scenario 1, except that va riation is assumed for

the parameters of the survival function.

2.8.1 General assumptions

1. The organism is semelparous.

2. Fecundity, F, increases with body size, x.

3. Survival, S, decreases with body size, x.

4. Two of the parameters in the functions describing the above two traits are

variable within generations but the same distribution occurs among genera-

tions.

5. Fitness, W, is a function of fecundity and survival.

2.8.2 Mathematical Assumptions

1. As before, fecundity increases linearly with body size:

F ¼ a

þ b

x ð2:42Þ

where a

and b

are constants.

94 MODELING EVOLUTION

2. Survival decreases linearly with body size:

S ¼ a

 b

x ð2:43Þ

where both a

and b

, vary within generations and are speciﬁed by independent

probability density functions. We shall assume the simplest probability function,

namely a uniform function, which is deﬁned as P (y) ¼ c, where c is a constant

determined from the minimum, y

min

, and maximum, y

max

, values of the parame-

ter y by the cumulative probability function:

max

min

cy ¼ 1

½y

max

min

¼ 1

cðy

max

 y

min

Þ¼1

c ¼ 1=ðy

max

 y

min

ð2:44Þ

In the present case the two functions are deﬁned as P(a

) ¼ c

and P(b

) ¼ c

3. As noted above, the appropriate measure of ﬁtness is the average of the ﬁtness

values:

W ¼

max

min

max

min

Pða

ÞPðb

Þða

þ b

xÞða

 b

xÞda

ð2:45Þ

Because a

and b

are independent one might be led to think that the mean ﬁtness

will indeed be equal to the ﬁtness using the mean values. This would indeed be

true if the variable parameters were a

and b

(i.e., the mean of a

þb

x is simply

equal to m

þ m

x, where m stands for the mean). However, survival must vary

within the range 0–1 and hence there are combinations in equation (2.45) that are

not permissible (and hence this equation as written is not strictly correct). For any

given x value the set of pairs of a

and b

must be restricted to those that ensure

that survival does not fall outside 0–1. For combinations in which survival is less

than zero, survival is set to zero and for combinations in which survival is greater

than 1, survival is set to 1.

2.8.3 Plotting the ﬁtness function

The basic approach is to generate two vectors containing values of a

and b

generated from the appropriate probability distributions (or multivariate distribu-

tion if the parameters are correlated). The uniform distribution is generated in

Rbyrunif and in MATLAB by rand. It is generally useful to set the seed prior

(set.seed(n) in R and rand(‘twister’, n) in MATLAB, where n is set by the

user) to the use of any random number generator as this allows the replication of

runs, which can be very useful in debugging. After generation of the two values

FISHERIAN OPTIMALITY MODELS 95

the vector of survival probabilities is calculated and values checked that they fall

within the acceptable range (0–1). In addition to calculating mean ﬁtness as a

function of x (body size), the following code also calculate ﬁtness using the mean

parameter values. For each parameter vectors of 1,000 random values are gener-

ated. These vectors are then used in calculating the ﬁtness for a given body size x.

The survival vector Surv is 1,000 units long. For simplicity and without loss of

signiﬁcant run time a loop is used to iterate over values of body size.

R CODE (Figure 2.6):

rm(list¼ls()) # Remove all objects from memory

Af <-2;Bf<- 2 # Invariant parameter values

Amin <- 0.3; Amax <- 1 # Min and max values of

Bmin <- 0; Bmax <- 0.2 # Min and max values of

Amean <- (AmaxþAmin)/2 # Mean value of

Bmean <- (BmaxþBmin)/2 # Mean value of

# Calculate n parameter combinations

n <- 1000 # Number of values of a

and b

to generate

# We are assuming a uniform distribution of values

set.seed(10) # Set the random number seed

# Generate n random numbers from Bmin to Bmax

Bs <- runif(n, min¼Bmin, max¼Bmax)

# Generate n random numbers from Amin to Amax

As <- runif(n, min¼Amin, max¼Amax)

x <- seq(from¼0, to¼6, length¼100) # Bodysizes from 0 to 6

W <- matrix(0,100,2) # Matrix to take fitness values

for ( i in 1: 100) # Iterate over x values

{

Surv <-AsBs*x[i] # Vector of survivals

# Check that no survival < 0. If so then set to zero

Surv[Surv<0] <-0

# Check that no survival > 1. If so then set to 1

Surv[Surv>1] <-1

# Column 1 contains fitness for variable parameters

W[i,1] <- mean((AfþBf*x[i])*Surv)

# Col 2 contains fitness using mean parameter values

W[i,2] <- (AfþBf*x[i])*(Amean-Bmean*x[i])

}

# Plot fitness¼W vs x for both columns on same graph

plot(x,W[,1], xlab¼’Body size, x’, ylab¼’Fitness, W’,las¼1,

lwd¼4) # Dots

lines(x,W[,2], lwd¼4) # Line

MATLAB CODE: See Section 2.18.13.

Notice that the two parameter vectors are outside of the loop over x.Thisisimpor-

tant for two reasons: ﬁrst, because only one pair of vectors is calculated, it speeds up

96 MODELING EVOLUTION

the program and second, the ﬁtnesses for the different values of x are compared

against the same set of parameters. It is quite evident from the two plots in Figure 2.6

that the optimum body size cannot be calculated using the parameter means.

2.8.4 Finding the maximum using the calculus

In some cases it may be possible to integrate the ﬁtness function and then solve

the resultant function either analytically or numerically. In this particular case

this is actually rather difﬁcult because the integration limits are actually functions

of x and there may be discontinuities because of setting survival equal to 1 if the

survival function exceeds 1. Integration of equation (2.45) as written actually leads

to the prediction that the optimal value of x is the same as that using the mean

parameter values, which is clearly incorrect (Figure 2.6). Here I shall use the

numerical integration routines in R and MATLAB. The multidimensional integra-

tion routine in R is called adapt and must be loaded as a separate package (copy

the ﬁle from the web source http://cran.r-project.org/web/packages/adapt/index.

html into the library folder of R). It was necessary to adjust some of the default

values in the routines adapt and nlm to get something approaching convergence

to a consistent value. The integration routine adapt works by generating a grid of

As and Bs values set by Amin, Amax, Bmin, and Bmax and using these to numeri-

cally evaluate the integral, where the function is deﬁned in the function INTEGRAND.

Body size, x

1.5

2.0

2.5

123456

Fitness, W

Figure 2.6 Scenario 6: Two measures of ﬁtness versus body size. Solid curve shows results

of using mean values (incorrect) and beaded lines shows correct curve.

FISHERIAN OPTIMALITY MODELS 97