R_using non linear regression with constraints

I am not as expert as John, but I thought it worth pointing out that the 
variable substitution technique gives up one set of constraints for 
another (b=0 in this case). I also find that plots help me see what is 
going on, so here is my reproducible example (note inclusion of library 
calls for completeness). Note that NONE of the optimizers mentioned so far
appear to be finding the true best fit. The fact that myfun() yields 0 
always if t=0 and that condition is within the data given seems likely to 
be part of the problem. I don't know how to resolve this... perhaps John 
will look at it again.

David: Your thinking makes fine sense if you are using a Monte Carlo or 
brute force solution, but the fact that it creates a discontinuity in 
the objective function will confuse any optimizer that uses analytic or 
numerically estimated slopes.

##----------begin
library(minpack.lm)
library(ggplot2)

mydata <- data.frame( x = c( 0, 5, 9, 13, 17, 20 )
                     , y = c( 0, 11, 20, 29, 38, 45 )
                     )

myfun <- function( a, b, r, t ) {
   a * b * ( 1 - exp( -b * r * t ) )
}

objdta <- expand.grid( a = seq( 1000, 3000, by=20 )
                      , b = seq( -0.01, 1, 0.01 )
                      , rowidx = seq.int( nrow( mydata ) )
                      )
objdta[ , c( "y", "t" ) ] <- mydata[ objdta$rowidx
                                    , c( "y", "x" ) ]
objdta$tf <- factor( objdta$t )
objdta$myfun <- with( objdta
                     , myfun( a = a, b = b, r = 2, t = t )
                     )
objdtass <- aggregate( ( objdta$myfun - objdta$y )^2
                      , objdta[ , c( "a", "b" ) ]
                      , FUN = function( x )
                               sum( x, na.rm=TRUE )
                      )
objdtassmin <- objdtass[ which.min( objdtass$x ), ]

myfit <- nlsLM( y ~ myfun( a, b, r=2, t=x )
               , data = mydata
               , start = list( a = 2000
                             , b = 0.05
                             )
               , lower = c( 1000, 0 )
               , upper = c( 3000, 1 )
               )
a <- as.vector( coef( myfit )[ "a" ] )
b <- as.vector( coef( myfit )[ "b" ] )

brks <- c( 500, 1e7, 2e7, 3e7, 4e7 )
ggplot( objdtass, aes( x=a, y=b, z = x, fill=x ) ) +
     geom_tile() +
     geom_contour( breaks= brks ) +
     geom_point( x=a, y=b, colour="red" ) +
     geom_point( x=objdtassmin$a
               , y=objdtassmin$b
               , colour="green" ) +
     scale_fill_continuous( name="SumSq", breaks = brks )
# Green point is brute-force solution
# Red point is optimizer solution for myfun

##############

myfun2 <- function( a, log1ab, r, t ) {
   ab <- 1000 - exp( log1ab )
   ab * ( 1 - exp( -ab/a * r * t ) )
}

objdta$log1ab <- with( objdta, log( 1000 - a * b ) )
objdta$myfun2 <- with( objdta
                      , myfun2( a = a
                              , log1ab = log1ab
                              , r = 2
                              , t = t
                              )
                      )
objdtass2 <- aggregate( ( objdta$myfun2 - objdta$y )^2
                       , objdta[ , c( "a", "b" ) ]
                       , FUN = function( x )
                                if ( all( is.na( x ) ) ) NA
                                else sum( x, na.rm=TRUE )
                       )
objdtass2min <- objdtass2[ which.min( objdtass2$x ), ]

myfit2 <- nlsLM( y ~ myfun2( a, log1ab, r = 2, t = x )
                , data = mydata
                , start = list( a = 2000
                              , log1ab = 4.60517
                              )
                , lower = c( 1000, 0 )
                , upper = c( 3000, 8.0063 )
                )
a2 <- as.vector( coef( myfit2 )[ "a" ] )
b2 <- ( 1000
       - exp( as.vector( coef( myfit2 )[ "log1ab" ] ) )
       ) / a

brks <- c( 500, 1e6, 2e6, 3e6, 4e6 )
ggplot( objdtass2, aes( x=a, y=b, z = x, fill=x ) ) +
     geom_tile() +
     geom_contour( breaks = brks ) +
     geom_point( x=a2, y=b2, colour="red" ) +
     geom_point( x=objdtass2min$a
               , y=objdtass2min$b
               , colour="green" ) +
     scale_fill_continuous( name="SumSq", breaks = brks )
# Green point is brute-force solution
# Red point is optimizer solution for myfun2

##----------end

R_using non linear regression with constraints

Thread (9 messages)