The *GenEstim* function presented here uses a very simple genetic algorithm to estimate parameters. The function returns the best estimated set of parameters ($estim), the AIC ($information) at each generation, and the cost of the best model ($bestcost) at each generation.

Results of running the program with a logistic function :

Logis = function(x,p) p[[1]]/(1+p[[2]]*exp(-p[[3]]*x))
RSS = function(par) sum((Logis(X,par)-Y)^2)
aic <- function(yvalues,rss,par)
{
k <- length(par)
n <- length(yvalues)
aic <- 2*k+n*log(rss/n)
return(aic)
}
P <- list(2,10,4)
X <- seq(from=-5,to=5,by=0.1)
Y <- Logis(X,P) + rnorm(length(X),sd=0.1)
plot(X,Y,pch=19,col='grey')
GenEstim <- function(
start.pars,
cost = RSS,
...,
numiter = 1e3,
npop = 1e2)
{
bestcost <- NULL
cur.AIC <- NULL
for(it in 1:numiter)
{
pop <- matrix(0,ncol=length(start.pars),nrow=npop)
for(p in 1:length(start.pars))
{
pop[,p] <- rnorm(npop,start.pars[[p]],sd=1)
pop[1,p] <- start.pars[[p]]
}
Costs <- NULL
for(i in 1:nrow(pop))
{
li <- as.list(pop[i,])
Costs[i] <- cost(li)
}
bestcost <- c(bestcost,min(Costs))
best <-which.min(Costs)
start.pars <- as.list(pop[best,])
cur.AIC <- c(cur.AIC,aic(Y,RSS(start.pars),start.pars))
}
return(list(estim=start.pars,information=cur.AIC,convergence=bestcost))
}
simul <- GenEstim(list(0,0,0))
x.control <- seq(from=-6,to=6,by=0.1)
lines(x.control,Logis(x.control,simul$estim),col='orange',lwd=3)

*Related*

To

**leave a comment** for the author, please follow the link and comment on his blog:

** [R] tricks**.

R-bloggers.com offers

**daily e-mail updates** about

R news and

tutorials on topics such as: visualization (

ggplot2,

Boxplots,

maps,

animation), programming (

RStudio,

Sweave,

LaTeX,

SQL,

Eclipse,

git,

hadoop,

Web Scraping) statistics (

regression,

PCA,

time series,

trading) and more...

If you got this far, why not

__subscribe for updates__ from the site? Choose your flavor:

e-mail,

twitter,

RSS, or

facebook...