how to propery specify a gradient function for use in optim() or other optimizer

To check if the gradient is correct, you can compare it with a numeric approximation:

library(numDeriv); grad( function(u) KannistoLik1( c(a=u[1], b=u[2]), .Dx, .Exp ), c(1,1) ); Kannisto.gr(c(a=1,b=1), .Dx, .Exp)

The signs are wrong: the algorithm does not see any improvement when it moves in this direction, and therefore does not move.

You can use some computer algebra system (here, Maxima)to do the computations for you:

display2d: false;f(a,b,x) := a * exp(b*x) / ( 1 + a * exp(b*x) );l(a,b,d,e,x) := - d * log(f(a,b,x)) + e * f(a,b,x);factor(diff(l(exp(a),exp(b),d,e,x),a));factor(diff(l(exp(a),exp(b),d,e,x),b));

I just copy and paste the result into R:

f_gradient <- function(u, .Dx, .Exp, .x.=.5:30.5) {  a <- u[1]  b <- u[1]  x <- .x.  d <- .Dx  e <- .Exp  c(    sum( (e*exp(exp(b)*x+a)-d*exp(exp(b)*x+a)-d)/(exp(exp(b)*x+a)+1)^2 ),    sum( exp(b)*x*(e*exp(exp(b)*x+a)-d*exp(exp(b)*x+a)-d)/(exp(exp(b)*x+a)+1)^2 )  )  }library(numDeriv)grad( function(u) KannistoLik1( c(a=u[1], b=u[2]), .Dx, .Exp ), c(1,1) )f_gradient(c(a=1,b=1), .Dx, .Exp)  # Identical

If you blindly put the gradient in the optimization,there is a numeric instability problem: the solution given is (Inf,Inf)...To prevent it, you can rescale the gradient (a better workaround would be to use a less explosive transformation than the exponential, to ensure that the parameters remain positive).

BFGSab <- optim(  log(c(a = .1, b = .1)),   fn = KannistoLik1,   gr = function(...) f_gradient(...) * 1e-3,   method = "BFGS",  .Dx = .Dx, .Exp = .Exp)exp(BFGSab$par) # Less precise than Nelder-Mead