## Monday, October 12, 2009

### Code I'm not proud of...

I was asked for the code I used to make my little pics. Well, it's not pretty, it's just a few functions in R. For example:
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] <- 1return(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.