### R code for 3-parameter Logistic model

 #pdf("rat42test.pdf")##___________________________________________________________________##          Stochastic Funnel Algorithm - SFA#             by #       Adrian O'Connor (aoconnor46@gmail.com)#___________________________________________________________________##   Using Logistic Growth model# (source: http://www.ncfaculty.net/dogle/fishR/gnrlex/VonBertalanffy/VonBertalanffy.pdf, page 45)#___________________________________________________________________## The SFA is a type of evolutionary algorithm that relies on a # trust region within the landscape of the SSE.  The trust region # is defined by the limits of the best 20 parameter sets based on # minimal SSE. New generations of candidates are created using a # uniform r.v. generator operating within the trust region. These # limits define rectangles that converge to a minimum SSE value as # the generations pass.# The SFA is a (100+20)EA, it creates 100 new points from the# trust region defined by the best 20 of the previous generation.# The resulting 120 points are sorted to get the best 20.# This code runs (simulates) the SFA for terminate criterion changing# from 1e+1 to 1e-9, dividing by 10 for each susequent generation.  # The resulting 11 scatterplots illustrate the convergence of the method.# Coded in R. Code available at http://tinyurl.com/9rtek7v.  #____________________________________________________________________## Initialisation including initial run of SFA#____________________________________________________________________#rm(list=ls())                                  # clear all definitions#x <- c(1:8)                           # data from NISTy <- c(0.1,0.2,1.2,3.1,5.1,6,5.9,6)#numrows <- length(x)                           # extract number of rows in vector x ( and y by default)#__________________________________________________________________________________##         Parameters for SFA search#__________________________________________________________________________________coeffdr       <-  0.4                        # Coefficient of Dynamic Relaxationcoeffsr       <-  1e-12                        # Coefficient of Static Relaxationterminate     <-  1e-1                        # Initial Termination conditionfinalterminate<-  1e-8                        # Final Termination condition                                              # N = # simulations of the SFAN <- log10(terminate)-log10(finalterminate)                                     m <- matrix(0,nrow=N,ncol=6)                  # initialise (=0)matrix of simualtion outputcolnames(m) <-c("Terminate","Gens","Beta1","Beta2","Beta3","SSE") # names of rows in msimulate <- 0                                 # counter for rows in m#__________________________________________________________________________________###  Function called "sse" to calculate the sse value at a point (beta1, beta2, beta3).#  This is a sum of squares whose value depends on beta1, beta2, beta3, the data & the model.  #  This function describes the topology of the sse landscape.#   #sse            <- function(beta1,beta2,beta3){       # begin sse function definition                sse <- 0                       # initialise sse                for(i in 1:numrows)                sse <- sse + (y[i]-beta1 / (1+exp(beta2-beta3 * x[i])))^2                return (sse)                                       }       # end sse function definition#sse(72.462,2.6181,0.06736)                     # Known global minimum of sse = 1168.009(source NIST)##  Generate initial vectors of beta1,beta2 and beta2#  These are distributed uniformly on a log10 scale between 10^-6 and 10^+6#  Actual parameter values extracted using "10^".#beta1          <- 10^runif(100,min=-6,max=+6)beta2          <- 10^runif(100,min=-6,max=+6)beta3          <- 10^runif(100,min=-6,max=+6)#sseval            <- sse(beta1,beta2,beta3)             # vector of the sse values# data.frame(beta1,beta2,beta3,sseval)                  # table of parameters & sse valuesssedata        <- data.frame(beta1,beta2,beta3,sseval)  # assign table to matrix ssedatasortedssedata  <- ssedata[order(ssedata\$sseval),] # sort table on basis of sse values## best 20 from sortedsedata on basis of minimal SSE#sortedbest20   <- sortedssedata[c(1:20),]      # assigned to vector#rangebeta1     <- c(min(sortedbest20[,1]),max(sortedbest20[,1]) )rangebeta2     <- c(min(sortedbest20[,2]),max(sortedbest20[,2]) )rangebeta3     <- c(min(sortedbest20[,3]),max(sortedbest20[,3]) )####___________________________________________________________________# #     Now start convergent search by updating beta value ranges#___________________________________________________________________count         <-0                             # initiate counter for generationswhile(terminate > finalterminate){# update row info for simulation matrixm[simulate,] <- c(terminate,count,sortedbest20[1,1],sortedbest20[1,2],sortedbest20[1,3],sortedbest20[1,4])#beta1          <- 10^runif(100,min=-6,max=+6)beta2          <- 10^runif(100,min=-6,max=+6)beta3          <- 10^runif(100,min=-6,max=+6)#sseval            <- sse(beta1,beta2,beta3)             # vector of the sse values# data.frame(beta1,beta2,beta3,sseval)                  # table of parameters & sse valuesssedata        <- data.frame(beta1,beta2,beta3,sseval)  # assign table to matrix ssedatasortedssedata  <- ssedata[order(ssedata\$sseval),] # sort table on basis of sse values## best 20 from sortedsedata on basis of minimal SSE#sortedbest20   <- sortedssedata[c(1:20),]      # assigned to vector#rangebeta1     <- c(min(sortedbest20[,1]),max(sortedbest20[,1]) )rangebeta2     <- c(min(sortedbest20[,2]),max(sortedbest20[,2]) )rangebeta3     <- c(min(sortedbest20[,3]),max(sortedbest20[,3]) )simulate <-simulate+1                         # update row counter for matrix of simulation output                                              # while loop below works on relative error in the sse. It relies on the sse # eventually being parabolic in the region of the minimum.#while((  sortedbest20[20,4]-sortedbest20[1,4])/sortedbest20[1,4]> terminate         & abs(sortedbest20[20,3]-sortedbest20[1,3])/min(sortedbest20[20,3],sortedbest20[1,3])> terminate       & abs(sortedbest20[20,2]-sortedbest20[1,2])/min(sortedbest20[20,2],sortedbest20[1,2])> terminate       & abs(sortedbest20[20,1]-sortedbest20[1,1])/min(sortedbest20[20,1],sortedbest20[1,1])> terminate) {                                              # while condition is relative error between 1st sse value & 20th sse valueoldbest20     <- sortedbest20                 # keep the best 20 from the previous generationspreadbeta1   <- rangebeta1[2] - rangebeta1[1]spreadbeta2   <- rangebeta2[2] - rangebeta2[1]spreadbeta3   <- rangebeta3[2] - rangebeta3[1]dynrelaxbeta1 <- coeffdr * spreadbeta1        # Dynamic relaxation for parametersdynrelaxbeta2 <- coeffdr * spreadbeta2dynrelaxbeta3 <- coeffdr * spreadbeta3count <- count+1                              # increment generations counter#                                             # clear vectors & matricesrm(beta1,beta2,beta3,ssedata,sortedssedata,sortedbest20,sse,sseval)#                                             # generate new beta valuestempbeta1min <- max(1e-6,rangebeta1[1]-dynrelaxbeta1-coeffsr)tempbeta1max <- min(1e+6,rangebeta1[2]+dynrelaxbeta1+coeffsr)tempbeta2min <- max(1e-6,rangebeta2[1]-dynrelaxbeta2-coeffsr)tempbeta2max <- min(1e+6,rangebeta2[2]+dynrelaxbeta2+coeffsr)tempbeta3min <- max(1e-6,rangebeta3[1]-dynrelaxbeta3-coeffsr)tempbeta3max <- min(1e+6,rangebeta3[2]+dynrelaxbeta3+coeffsr)#                                             # now use logs to generate uniformly on a log scaletemplogbeta1 <- runif(100,min=log10(tempbeta1min),max=log10(tempbeta1max))templogbeta2 <- runif(100,min=log10(tempbeta2min),max=log10(tempbeta2max))templogbeta3 <- runif(100,min=log10(tempbeta3min),max=log10(tempbeta3max))#                                             # extract vectors of beta values using powers of 10beta1         <- 10^templogbeta1beta2         <- 10^templogbeta2beta3         <- 10^templogbeta3#sse           <- function(beta1,beta2,beta3){ # begin sse function definition                sse <- 0                      # initialise sse                for(i in 1:numrows)                sse <- sse + (y[i]-beta1 / (1+exp(beta2-beta3 * x[i])))^2                return (sse)                                       }      # end sse function definitionsseval           <- sse(beta1,beta2,beta3)    # vector of the sse values# data.frame(beta1,beta2,beta3,sseval)        # table of parameters & sse values#                                             # assign table to vector ssedatassedata       <- data.frame(beta1,beta2,beta3,sseval)  ssedata       <- rbind(ssedata,oldbest20)     # concanenate newdata with best 20 from previous generationsortedssedata <- ssedata[order(ssedata\$sseval),]  # sort table of 120 rows on basis of sse values                                              # only the best 20 rows (lowest sse) are now actually used.## best 20 from sortedssedata on basis of minimal SSE#sortedbest20  <- sortedssedata[c(1:20),]      # assigned to vector sortedbest20##                                             # create ranges using min and max parameter values from best20#rangebeta1    <- c(min(sortedbest20[,1]),max(sortedbest20[,1]) )rangebeta2    <- c(min(sortedbest20[,2]),max(sortedbest20[,2]) )rangebeta3    <- c(min(sortedbest20[,3]),max(sortedbest20[,3]) )#par(ask=TRUE)                                 # ask to print plots#}                                             # end of while loop for each terminate criterion#sortedbest20par(mfrow=c(1,3))                             # set three plots per pagepar(pch=19)                                   # default plot chars (closed circles)                                              # these two plot statements create a single scatterplot of the parameters                                              # with the best 20 in redplot(sortedssedata[,1],sortedssedata[,2],xlab="Beta1",ylab="Beta2",    xlim=c(min(sortedssedata[,1]),max(sortedssedata[,1])),    ylim=c(min(sortedssedata[,2]),max(sortedssedata[,2])),    main=paste("Terminate Criterion = ",as.character(terminate))    )#par(new=TRUE,pch=19)                          # overprint existing plot (new=TRUE) with red (col=2)filled circles(pch=19)                                              # for the best 20 parameter sets on basis of sse.                                              # Note scales and labels are the same for both plots.#plot(sortedbest20[,1],sortedbest20[,2],xlab="Beta1",ylab="Beta2",    xlim=c(min(sortedssedata[,1]),max(sortedssedata[,1])),    ylim=c(min(sortedssedata[,2]),max(sortedssedata[,2])),    col=2    )par(new=FALSE)                                # reset overprint parameterplot(sortedssedata[,1],sortedssedata[,3],xlab="Beta1",ylab="Beta3",    xlim=c(min(sortedssedata[,1]),max(sortedssedata[,1])),    ylim=c(min(sortedssedata[,3]),max(sortedssedata[,3])),    main=paste("Terminate Criterion = ",as.character(terminate))    )#par(new=TRUE,pch=19)                          # overprint existing plot (new=TRUE) with red (col=2)filled circles(pch=19)                                              # for the best 20 parameter sets on basis of sse.                                              # Note scales and labels are the same for both plots.#plot(sortedbest20[,1],sortedbest20[,3],xlab="Beta1",ylab="Beta3",    xlim=c(min(sortedssedata[,1]),max(sortedssedata[,1])),    ylim=c(min(sortedssedata[,3]),max(sortedssedata[,3])),    col=2    )par(new=FALSE)                                # reset overprint parameterplot(sortedssedata[,2],sortedssedata[,3],xlab="Beta2",ylab="Beta3",    xlim=c(min(sortedssedata[,2]),max(sortedssedata[,2])),    ylim=c(min(sortedssedata[,3]),max(sortedssedata[,3])),    main=paste("Terminate Criterion = ",as.character(terminate))    )#par(new=TRUE,pch=19)                          # overprint existing plot (new=TRUE) with red (col=2)filled circles(pch=19)                                              # for the best 20 parameter sets on basis of sse.                                              # Note scales and labels are the same for both plots.#plot(sortedbest20[,2],sortedbest20[,3],xlab="Beta2",ylab="Beta3",    xlim=c(min(sortedssedata[,2]),max(sortedssedata[,2])),    ylim=c(min(sortedssedata[,3]),max(sortedssedata[,3])),    col=2    )par(new=FALSE)                                # reset overprint parameter#plot(sortedssedata,          main=paste("Terminate Criterion = ",as.character(terminate)))terminate     <-terminate/10                  # update terminate criterion}                                             # end of while loop for limit of terminate criterion#plot(sortedssedata,          main="Final run of SFA")            # multiple scatterplots of final generation of  b1, b2 & sseplot(sortedbest20,         main="Best 20 of Final run of SFA",col="Red")  # multiple scatterplots of best20 from final generation                                              # print matrix m of simulations of SFAprint("      Matrix of simulations of the SFA")m#dev.off()##                                             # Parameter estimation using "nls" #                                             # and starting values from final SFA simulation"# modelout     <-  nls(y~beta1/(1+exp(beta2 - beta3 * x)),                              data=data.frame(x,y),                              start=list(beta1=m[N,3],beta2=m[N,4],beta3=m[N,5]),                              trace=TRUE)                                              # Trace shows starting values and path to final estimates.modelout#                                             # Plot of data and best fit modelpar(new=FALSE,mfrow=c(1,1))plot(x,y,xlim=(c(0,8)),ylim=(c(0,8)),ylab="Response",main="3-parameter logistic model",sub="data & model",col="Red")par(new=TRUE)curve(coef(modelout)[1]/(1+exp(coef(modelout)[2]-coef(modelout)[3] * x)),xlim=(c(0,8)),ylim=(c(0,8)),ylab="Response")