Recent site activity

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 NIST
y <- 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 Relaxation
coeffsr       <-  1e-12                        # Coefficient of Static Relaxation
terminate     <-  1e-1                        # Initial Termination condition
finalterminate<-  1e-8                        # Final Termination condition
                                              # N = # simulations of the SFA
N <- log10(terminate)-log10(finalterminate)                                     
m <- matrix(0,nrow=N,ncol=6)                  # initialise (=0)matrix of simualtion output
colnames(m) <-c("Terminate","Gens","Beta1","Beta2","Beta3","SSE") # names of rows in m
simulate <- 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 values
ssedata        <- data.frame(beta1,beta2,beta3,sseval)  # assign table to matrix ssedata

sortedssedata  <- 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 generations

while(terminate > finalterminate){
# update row info for simulation matrix
m[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 values
ssedata        <- data.frame(beta1,beta2,beta3,sseval)  # assign table to matrix ssedata

sortedssedata  <- 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 value
oldbest20     <- sortedbest20                 # keep the best 20 from the previous generation

spreadbeta1   <- rangebeta1[2] - rangebeta1[1]
spreadbeta2   <- rangebeta2[2] - rangebeta2[1]
spreadbeta3   <- rangebeta3[2] - rangebeta3[1]

dynrelaxbeta1 <- coeffdr * spreadbeta1        # Dynamic relaxation for parameters
dynrelaxbeta2 <- coeffdr * spreadbeta2
dynrelaxbeta3 <- coeffdr * spreadbeta3


count <- count+1                              # increment generations counter
#                                             # clear vectors & matrices
rm(beta1,beta2,beta3,ssedata,sortedssedata,sortedbest20,sse,sseval)
#                                             # generate new beta values
tempbeta1min <- 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 scale
templogbeta1 <- 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 10
beta1         <- 10^templogbeta1
beta2         <- 10^templogbeta2
beta3         <- 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 definition

sseval           <- sse(beta1,beta2,beta3)    # vector of the sse values
# data.frame(beta1,beta2,beta3,sseval)        # table of parameters & sse values
#                                             # assign table to vector ssedata
ssedata       <- data.frame(beta1,beta2,beta3,sseval)  
ssedata       <- rbind(ssedata,oldbest20)     # concanenate newdata with best 20 from previous generation

sortedssedata <- 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
#sortedbest20
par(mfrow=c(1,3))                             # set three plots per page

par(pch=19)                                   # default plot chars (closed circles)
                                              # these two plot statements create a single scatterplot of the parameters
                                              # with the best 20 in red
plot(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 parameter

plot(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 parameter

plot(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 & sse

plot(sortedbest20,
         main="Best 20 of Final run of SFA",col="Red")  # multiple scatterplots of best20 from final generation
                                              # print matrix m of simulations of SFA
print("      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 model
par(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")

Comments