R Under development (unstable) (2023-12-07 r85661 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(bbmle)
Loading required package: stats4
> 
> ## Simulate data
> 
> set.seed(1)
> x <- 1:5
> y <- 2*x+1
> noise <- rnorm(5, 0, 0.1)
> mydata <- data.frame(x = x, y=y+noise)
> 
> ## Model definition
> 
> model <- function(a, b) with(mydata, a*x+b)
> 
> ## Negative log-likelihood
> 
> nll <- function(par) with(mydata, {
+   a <- par[1]
+   b <- par[2]
+   sum(0.5*((y-model(a,b))/0.1)^2)
+   
+ })
> 
> gr <- function(par) with(mydata, {
+   a <- par[1]
+   b <- par[2]
+   dnllda <- -sum(((y-model(a,b))/0.1)*x/0.1)
+   dnlldb <- -sum(((y-model(a,b))/0.1)*1/0.1)
+   return(c(dnllda, dnlldb))
+ })
> 
> ## optimization
> 
> parnames(nll) <- c("a", "b")
> parnames(gr) <- c("a", "b")
> 
> fit <- mle2(nll, c(a = 1, b=2), gr=gr)
> 
> myprof <- profile(fit)
> myprof_c <- profile(fit,continuation="naive")
> confint(myprof)
      2.5 %   97.5 %
a 1.9712561 2.095215
b 0.7076574 1.118783
> confint(myprof_c)
      2.5 %   97.5 %
a 1.9712561 2.095215
b 0.7076574 1.118783
> 
> fit <- mle2(nll, c(a = 1, b=2), gr=gr, skip.hessian=TRUE)
> myprof2 <- profile(fit,std.err=c(0.1,0.1))
> 
> ## incomplete!
> model2 <- ~a+b*x+c*x^2
> f0 <- deriv(model2,"x",function.arg=c("a","b","c"))
> ## chain rule
> f1 <- function() {
+ ## memoize
+ lastpar <- NULL
+ lastval <- NULL
+ }
> 
> f2 <- function(par) {
+     if (par==lastpar) {
+         return(c(lastval))
+     }
+     lastpar <<- par
+     lastval <<- do.call(f0,par)
+     f1(par)
+ }
> f2.gr <- function(par) {
+     if (par==lastpar) {
+         return(attr(lastval,".grad"))
+     }
+     lastpar <<- par
+     lastval <<- do.call(f0,par)
+     f1.gr(par)
+ }
> parnames(f2) <- parnames(f2.gr) <- c("a","b","c")
> 
> proc.time()
   user  system elapsed 
  12.57    0.14   12.70