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

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

Fxt_upper ¼ F_vectors(Index_upper,2);

Fxt_interpolated ¼ Qx*Fxt_upper þ (1-Qx)*Fxt_lower; % Iterpo-

lated value

% Calculate the fitness for this particular egg-clutch size combi-

nation

W ¼ Wmax(Ith_Patch)-sqrt(A(1)*(Egg-Xegg(Ith_Patch))^2 þ A(2)*

(Clutch-Xclutch(Ith_Patch))^2);

W ¼ max(0, W); % Set to zero if negative

W ¼ W þ Psurvival*Fxt_interpolated; % Fitness

W ¼ max(0, W); % Set to zero if negative

end % End of if

W ¼ -W; % Make Fitness negative for fminbnd function

% ———————————————————————————————————————————————————————————

function Temp¼OVER_PATCHES(X, F_vectors, Xcritical, Xmax,

Xinc, Npatch, Psurvival, Wmax, A, Xegg, Xclutch,P1)

% Function to iterate over patches i_e_ over Hosts

Index ¼ 1þ(X-Xcritical)/Xinc;

% Allocate storage of best combinations for each patch

% Columns will contain

% Fitness, Clutch, Egg

Choice_Flag ¼ zeros(2,1); % Store information on number of

choices

Best_in_Patch ¼ zeros(2,3); % Allocate storage for Best Decision

% Iterate over patches

for Ith_Patch ¼ 1:Npatch % Iterate over the two hosts (¼ patches)

% Make a matrix called W_host with the following 3 columns:

% Fitness, Egg size, Clutch size

W_host ¼ zeros(11,3);

for Clutch ¼ 1:11 % Iterate over clutch size

W_host(Clutch,2) ¼ Clutch; % Store clutch size

% Call fminbnd to find best egg size Note that this is finding a

minimum

x ¼ fminbnd(@(x) FITNESS(x,Clutch, X,F_vectors, Xcritical,

Xmax, Xinc, Psurvival, Wmax,A,Xegg,Xclutch,Ith_Patch),0.01,3);

W_host(Clutch,1) ¼ -FITNESS(x,Clutch, X,F_vectors, Xcritical,

Xmax, Xinc, Psurvival, Wmax, A, Xegg, Xclutch,Ith_Patch); % Fitness

W_host(Clutch,3) ¼ x; % Egg size

end % End of clutch size loop

% Get best combination for this host

% Use sort because we need to inspect best two

R ¼ W_host(:,1);

[Sorted_W_host,I] ¼sort(R, ‘descend’); % Sortsinto descending col

Best_in_Patch(Ith_Patch,:) ¼ W_host(I(1),:); % Store best

choice

% Test for several equal optimal choices

418 MODELING EVOLUTION

if(W_host(I(1),1)¼¼W_host(I(2),1))

Choice_Flag(Ith_Patch) ¼ 2;

end % end if

end % Next host

% Overall fitness

W ¼ P1*Best_in_Patch(1,1) þ(1-P1)*Best_in_Patch(2,1);

F_vectors(Index,1)= W ; % Update F(x,t,T)

% Concatenate F(x,t) and the optimal egg and clutch values for both

hosts

% We add to the bottom of the two column matrix F_vectors the follow-

ing

% F_vectors(Index,1), 1 The second entry is simply a dummy variable

% Best_in_Patch(1,2) Best_in_Patch(1,3) % Host 1 Egg size Clutch

size

% Best_in_Patch(2,2) Best_in_Patch(2,3) % Host 2 Egg size Clutch

size

% Choice_Flag(1:2) % Flag for multiple optima

Temp1 ¼ horzcat(F_vectors(Index,1), 1);

Temp2 ¼ horzcat( Best_in_Patch(1,2), Best_in_Patch(1,3));

Temp3 ¼ horzcat( Best_in_Patch(2,2), Best_in_Patch(2,3));

% Add Temp1, Temp2, Temp3 & Choice to bottom of F_vectors and rename

to Temp

Choice ¼ Choice_Flag’ % Convert into column matrix

Temp ¼ vertcat(F_vectors, Temp1, Temp2, Temp3, Choice);

% End of function

%———————————————————————————————————————————————————————————

function Temp¼OVER_STATES(F_vectors, Xcritical, Xmax, Xinc,

Npatch, Psurvival, Max_Index, Wmax, A, Xegg, Xclutch, P1)

% Function to iterate over states of X

% Create matrix for output. Note that we use seven columns

Store ¼ zeros(Max_Index,7);

% Iterate over X ¼ Biomass X(1) is zero so skip

for Index ¼ 2 : Max_Index % Iterate over states of X

% For given X call Over_Patches to determine F(x,t) and best patch

X ¼ (Index-1)*Xinc þ Xcritical;

Temp ¼ OVER_PATCHES(X, F_vectors, Xcritical, Xmax, Xinc,

Npatch, Psurvival,Wmax, A, Xegg, Xclutch, P1);

% Extract components_ Last row-2 is F(x,t) and dummy variable

% Last row-1 is best clutch and egg size for host type 1

% Last row is best clutch and egg size for host type 2

% Last row is flage indicating multiple equal choices

n ¼ size(Temp,1)-4;

F_vectors ¼ Temp(1:n,:); % Extracting F_vectors

% Add the seven output values (omit dummy) to storage

DYNAMIC PROGRAMMING 419

Store(Index,:) ¼ horzcat(Temp(nþ1,1), Temp(nþ2,1:2), Temp

(nþ3,1:2), Temp(nþ4, 1:2));

end % End of X loop

% Add Store values to end of F_vectors for pass back to main program

Temp ¼ horzcat(F_vectors, Store); % Combined by columns

% End of function

%————————————————————————————————————————————————————————————

% MAIN PROGRAM

clear all % Remove all objects from memory

% Initialize parameters

% Create the state variable (X¼reproductive biomass)

Xmax ¼ 10; Xcritical ¼ 0; Xinc ¼ 1;

Max_Index ¼ 1 þ (Xmax-Xcritical)/Xinc;

% Parameter values on the two hosts

Wmax ¼ [10, 20]; % Maximum fitness

A ¼ [100, 1]; % a coefficients

Xegg ¼ [2, 1]; % “optimum” egg size

Xclutch ¼ [5, 10]; % “Optimum” clutch size

P1 ¼ 0.5; % Probability of host 1

Psurvival¼ 0.90; % Survival probability per time increment

Npatch ¼ 2; % Number of patches¼ hosts

Horizon ¼ 10; % Number of time steps

% Set up matrix for fitnesses

% Column 1 is temporary F(x, t)_ Column 2 is F(x, tþ1) Both are zero

F_vectors ¼ zeros(Max_Index,2) ;

% Create matrices for output

FxtT ¼ zeros(Horizon,Max_Index);% F(x,t)

Best_Clutch1 ¼ zeros(Horizon,Max_Index);% Best clutch size

for host 1

Best_Clutch2 ¼ zeros(Horizon,Max_Index);% Best clutch size for

host 2

Best_Egg1 ¼ zeros(Horizon,Max_Index);% Best egg size for

host 1

Best_Egg2 ¼ zeros(Horizon,Max_Index);% Best egg size for

host 2

Choice_H1 ¼ ones(Horizon,Max_Index);% 1 or 2 choices for host

type 1

Choice_H2 ¼ ones(Horizon,Max_Index);% 1 or 2 choices for host

type 2

% Start iterations

Time ¼

Horizon; % Initialize Time

while

( Time > 1)

Time ¼ Time  1; % Decrement Time by 1 unit

% Call OVER_STATES to get best values for this time step

420 MODELING EVOLUTION

Temp ¼ OVER_STATES(F_vectors, Xcritical, Xmax, Xinc, Npatch,

Psurvival, Max_Index, Wmax, A, Xegg, Xclutch, P1);

% Extract F_vectors

TempF ¼ Temp(:,1:2);

% Update F1

for J ¼ 2: Max_Index

F_vectors(J,2) ¼ TempF(J,1);

end % end j loop

% Store results

FxtT(Time,:) ¼ Temp(:,3);

Best_Clutch1(Time,:) ¼ Temp(:,4);

Best_Egg1(Time,:) ¼ Temp(:,5);

Best_Clutch2(Time,:) ¼ Temp(:,6);

Best_Egg2(Time,:) ¼ Temp(:,7);

Choice_H1(Time,:) ¼ Temp(:,8);

Choice_H2(Time,:) ¼ Temp(:,9);

end % End of Time loop

% Output information_ For display add states to last row of matrices

Index ¼ 1:Max_Index;

FxtT(Horizon,:) ¼ (Index-1)*XincþXcritical;

Best_Clutch1(Horizon,:) ¼ (Index-1)*XincþXcritical;

Best_Clutch2(Horizon,:)¼ (Index-1)*XincþXcritical;

Best_Egg1(Horizon,:) ¼ (Index-1)*XincþXcritical;

Best_Egg2(Horizon,:) ¼ (Index-1)*XincþXcritical;

Choice_H1(Horizon,:) ¼ (Index-1)*XincþXcritical;

Choice_H2(Horizon,:) ¼ (Index-1)*XincþXcritical;

% Print Best clutch, best egg choice flag and decision matrix

Best_Clutch1(:,1:Max_Index) % Best clutch on host 1

Best_Clutch2(:,1:Max_Index) % Best clutch on host 2

vpa(Best_Egg1(:,1:Max_Index),3) % Best egg on host 1

vpa(Best_Egg2(:,1:Max_Index),3) % Best egg on host 2

Choice_H1(:,1:Max_Index) % 1 ¼only 1 choice 2 ¼1þ

Choice_H2(:,1:Max_Index) % 1 ¼only 1 choice 2 ¼1þ

vpa(FxtT(:,1:Max_Index),3)% Print Fxt of Decision matrix: 3 sig

places

DYNAMIC PROGRAMMING 421

This page intentionally left blank

Appendix 1

R functions used in this book with equivalent MATLAB codes where applicable.

Except for the ﬁrst four rows, rows are arranged in alphabetical order of the R code

(na ¼ not applicable for this book).

R MATLAB Operation

# Comment % Comment Comment

...................................................................................................................................................

<- = Assigns variable

...................................................................................................................................................

X%*%Y X*Y Matrix multiplication

of X and Y

...................................................................................................................................................

X*Y X.*Y Element by element

multiplication

...................................................................................................................................................

abs(X) abs(X) Absolute value of X

...................................................................................................................................................

adapt dblquad Multidimensional

integration

...................................................................................................................................................

apply use loop construct Applies function to

rows, cols or both

...................................................................................................................................................

as.numeric(X) na Make X numeric

...................................................................................................................................................

barplot(….) bar(…) Make a bar plot

...................................................................................................................................................

Best <- X[Y[1]] Best=X(I(1))

...................................................................................................................................................

cbind horzcat Bind columns

...................................................................................................................................................

ceiling(X) ceil(X) Nearest larger integer

...................................................................................................................................................

contour.(x, y, z.

matrix)

contour(x, y,

zmatrix)

Contour plot

...................................................................................................................................................

contourLines na Set of x,y coordinates

for user-specified

contour lines

...................................................................................................................................................

cor(X,Y) Correlation between X

and Y

...................................................................................................................................................

c(x1, x2, x3…xn) [x1, x2, x3…..xn] Concatenate into a

vector

...................................................................................................................................................

data.frame na Matrix that can take

variables of mixed type

...................................................................................................................................................

dbinom(x=k,

size=3, prob=Pi

(PackSize))

binopdf(k, 3, Pi

(PackSize))

Binomial distribution

...................................................................................................................................................

deriv diff Derivative

...................................................................................................................................................

(continued)

R MATLAB Operation

diag(X) diag(X,0) Diagonal elements of

matrix X

...................................................................................................................................................

dnorm(x=, mean=,

sd= )

normpdf(x, mu, sd) Density at x of the

normal distribution

...................................................................................................................................................

eigen(X) eig Eigen values and other

stats of X

...................................................................................................................................................

exp(x) exp(x) e^x

...................................................................................................................................................

expand.grid(X,Y) meshgrid Matrix of all

combinations of X and Y

...................................................................................................................................................

floor(X) floor(X) Integer part of X

...................................................................................................................................................

for (i in 1:n) for i = 1:n Loop

{ lines of coding} lines of coding

end

...................................................................................................................................................

FUNCTION NAME <−

function(Input

parameters)

function Output

=FUNCTION NAME

(Input parameters)

Function definition

{ Lines of code Last line gives

output parameters

Return output

parameters(s) }

...................................................................................................................................................

X <- FUNCTION NAME

(input

parameters)

X = FUNCTION NAME

(Input parameters)

Function call

...................................................................................................................................................

hist(X) hist(X) Histogram of X

...................................................................................................................................................

if (X > Y) if X>y If construct

{ Lines of coding

Lines of coding end

}

...................................................................................................................................................

if (X > Y) if X > Y If else construct

{ Lines of coding

Lines of coding} end

} else

else{ Lines of

coding

lines of coding

} end

...................................................................................................................................................

Lines of coding}

...................................................................................................................................................

integrate int, quad Integration

...................................................................................................................................................

length(X) length(X) Number of elelments

in X

...................................................................................................................................................

lines(x,y,

optional

parameters)

plot Plot a line on an

existing plot

...................................................................................................................................................

lm(y~x) na Linear regression of y

on x

...................................................................................................................................................

424 MODELING EVOLUTION

R MATLAB Operation

log(X) log(X) Natural log of X

...................................................................................................................................................

matrix(0,n) zeros(n) Column matrix of n

zeros

...................................................................................................................................................

matrix(0,n1,n2) zeros(n1,n2) n1 x n2 matrix of zeros

...................................................................................................................................................

max(X) max(X) Maximum of values in X

...................................................................................................................................................

mean(X) mean(X) Mean of values in X

...................................................................................................................................................

min(X) min(X) Minimum of values in X

...................................................................................................................................................

mvrnorm(n=, mu=,

Sigma= )

mvnrnd(MU,SIGMA) Random multivariate

normal deviates

...................................................................................................................................................

nlm, optimize fminsearch, fminbnd Non-linear

optimization

...................................................................................................................................................

ncol(X) size(X,2) Nos of columns in X

...................................................................................................................................................

nrow(X) size(X,1) Nos of rows in X

...................................................................................................................................................

order(X) [C,I]=max(X) Get indexes of X that

indicate sort

structure

...................................................................................................................................................

outer(x, y, FUNC=) na All combinations of x

and y applying

FUNC=function

...................................................................................................................................................

par(mfcol=c

(1,2))

subplot(2,1,1) Divide graphics page in

two

...................................................................................................................................................

par(mfrow=c

(2,2))

subplot(2,2,2) ivide graphics page

into quadrats

...................................................................................................................................................

persp(x, y, z.

matrix)

surfc(x, y, zmatrix) 3D plot

...................................................................................................................................................

plot(x,y,

optional

parameters)

plot, ezplot, fplot 2D plot

...................................................................................................................................................

pnorm(x, mean=,

sd= )

1−normcdf(X,mu,

sigma)

Integral from x to

infinity of normal

distribution

...................................................................................................................................................

points(x, y) hold on plot(x,y,´:´) Plot points on an

existing plot

...................................................................................................................................................

print(X) X Print X

...................................................................................................................................................

print(c(x1, x2,

x3))

[x1, x2, x3] Print x1, x2, x3

...................................................................................................................................................

qnorm(P, mean=,

sd= )

na normal distribution

function

...................................................................................................................................................

read.table

(file=“Filename.

txt”)

load Filename.txt Read data from text

file

...................................................................................................................................................

(continued)

APPENDIX 1 425

R MATLAB Operation

rep(x, times= ) na Replicate x for times=

...................................................................................................................................................

rbind vertcat Bind rows

...................................................................................................................................................

rm(list=ls()) clear all Clears workspace

...................................................................................................................................................

rnorm(n, mean=,

sd= )

na Random normal deviates

...................................................................................................................................................

rowSums(X) na Row sums of X

...................................................................................................................................................

round(X) round(X) Round X to nearest

integer

...................................................................................................................................................

rpois(n=, lambda=

)

na Random deviate from a

Poisson distribution

...................................................................................................................................................

runif(n, min=X1,

max=X2 )

X1 + (X2−X1)*rand

(n,1)

n uniform random

numbers between X1 and

...................................................................................................................................................

sample(x=, size=,

replace= )

na Random sample from X

...................................................................................................................................................

sd(X) Standard deviation of X

...................................................................................................................................................

setwd(“C:/……”) set directory in gui Set path for output

...................................................................................................................................................

set.seed(x) rand(‘twister’, x) Set seed for random

number

...................................................................................................................................................

seq(from=x1,

to=x2, length=n)

linspace(x1, x2, n) Sequence generator

...................................................................................................................................................

seq(from=x1,

to=x2, 1)

x1:x2 Sequence generator

...................................................................................................................................................

smooth.spline na Fits a cubic smoothing

spline to the supplied

data

...................................................................................................................................................

solve(X) Invert matrix X

...................................................................................................................................................

sqrt(X) sqrt(X) Square-root of X

...................................................................................................................................................

sum symsum Sum X

...................................................................................................................................................

summary(X) na Summary statistics of X

...................................................................................................................................................

table(X) na Tabulate X

...................................................................................................................................................

t(X) na Transpose X

...................................................................................................................................................

t.test(X,

optional

parameters)

na t test

...................................................................................................................................................

uniroot fzero, solve Find root

...................................................................................................................................................

unlist(X) na reduce to a vector of

atomic components of x

...................................................................................................................................................

426 MODELING EVOLUTION

R MATLAB Operation

var(X) na Variance of X

...................................................................................................................................................

while (Condition) while (Condition) While loop

{ Lines of coding Lines of coding

} end

...................................................................................................................................................

write(X,

file=“filename.

txt”)

save Filename.txt Write to text file See

setwd

APPENDIX 1 427