weasel.matrix.solo <- function(target.length=5, alphabet.length=27,mut.prob=.04,mut.effect=FALSE){

if(mut.effect) mut.prob <- alphabet.length*mut.prob/(alphabet.length-1)

M <- matrix(integer((target.length+1)*(target.length+1)),nrow=target.length+1)

for(k in 1:(target.length)){

verschlechtert <- integer(2*target.length+2)

verbessert <- verschlechtert

verschlechtert[1:k + target.length +1] <- dbinom(0:(k-1),k-1,mut.prob*(alphabet.length-1)/alphabet.length)

verbessert[1:(target.length+2-k) + target.length + 1] <- dbinom(0:(target.length+1-k),(target.length+1-k),mut.prob/alphabet.length)

for(kk in 1:(target.length+1)){

M[k,kk] <-sum( verbessert[(0:target.length) + target.length+3-k]*verschlechtert[0:target.length +target.length+3-kk])

}

}

M[target.length+1,target.length+1] <- 1

return(M)

}

This creates a transition matrix for Dawkins's algorithm with a single child.

`target.length`

is the length of the target phrase, `alphabet.length`

is self-explaining, `mut.prob`

is the probability with which a letter is mutated: if `mut.effect==FALSE`

, such a mutation has no effect with probability 1/N, i.e., it's a change to the same letter :-)

weasel.matrix.cum <- function(M, pop.size = 10){

target.length <- dim(M)[1]-1

M1 <- M

for(k in 1:(target.length)){

M1[k,] <- cumsum(M[k,])

}

M1 <- M1^pop.size

for(k in 1:(target.length))

M1[k,] <- M1[k,] - c(0,M1[k,-target.length-1])

return(M1)

}

creates a transition matrix for a whole population of size `pop.size`

, given a transition matrix `M`

for a single-child process.The last thing which is needed is this little snippet:

`erw.exakt <- function(M,start.vektor){`

M.size <- dim(M)[1]

M <- M[-M.size, -M.size]

M <- solve(diag(M.size-1)-M)

M <- apply(M,1,sum)

return(sum(start.vektor[1:(M.size-1)] * M))

}

Here, the expected value of generations for a transition matrix

`M`

is calculated, using a initial distribution for `start.vektor`

(mostly Bin(target.length,mut.prob) distributed - or just the first unit vector).If I ever prettify this code, I'll take Tom English on his generous offer to put it somewhere on the net.

## No comments:

## Post a Comment