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)