Nothing
library(bbmle)
old_opt <- options(digits=3)
## source("../R/dists.R")
## source("../R/mle.R")
## an attempt to sketch out by hand
## how one would derive an analytic
## gradient function for a formula-specified
## likelihood and use it ...
## chain rule should be:
## deriv(probability distribution)/[prob params] *
## deriv([prob params])/[model params] *
## {OPTIONAL} deriv([model params])/[linear model params]
set.seed(1001)
x <- rbinom(50,size=10,prob=0.4)
suppressWarnings(mle2(x~dbinom(prob=p,size=10),start=list(p=0.3),data=data.frame(x)))
## step 1: construct gradient function for simplest example
f <- sbinom(prob=0.1,size=1)$formula
d1 <- deriv(parse(text=f),"prob",function.arg=TRUE)
## step 2: chain rule step #1
mle2(x~dbinom(prob=plogis(logitp),size=10),start=list(logitp=-1),
data=data.frame(x))
f <- sbinom(prob=NA,size=NA)$formula
## note: plogis is not in derivatives table!!
## will need to extend by text substitution ...
gsub("plogis(\\([^)]+\\))",
"(1+exp(\\1))^(-1)",
"plogis(logitprob)")
f2 <- gsub("plogis(\\([^)]+\\))",
"(1+exp(\\1))^(-1)","plogis(logitp)")
## start with a single parameter (ignore 'size')
fun1 <- deriv(parse(text=f),c("prob"),function.arg=TRUE)
fun2 <- deriv(parse(text=f2),"logitp", function.arg=TRUE)
size <- 10
a1 <- attr(fun2(logitp=0),"gradient")
a2 <- attr(fun1(prob=plogis(0)),"gradient")
## compute gradient by variable and sum
colSums(apply(a1,2,"*",a2))
## rep(a1,length(x))*a2
## eventually we will want to do something tricky to
## 'memoise' results because optim() requires the
## objective function and gradient to be computed
## *separately*. Not worth worrying about this in the
## first pass!
options(old_opt)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.