CODEWITHSUNDEEP
rnonlinear-optimization
Nick Mars
Nick Mars

Reputation: 117

non-linear optimization in R using optim

I'm a newbie in R!

I would like to find the best gamma distribution parameters to fit my experimental counts data. The optim function's help file says the first argument of the function should be the parameters to be optimized. So I tried :

x = as.matrix(seq(1,20,0.1))
yexp = dgamma(x,2,1)*100 + rnorm(length(x),0,1)
f = function(p,x,yexp) {sum((p[1]*dgamma(x,p[2],scale=p[3]) - yexp)^2)}
mod = optim(c(50,2,1),f(p,x,yexp))

I get the error message :

Error in f(p, x, yexp) : object 'p' not found

Any hint where I'm wrong?

Supplementary question : Is there any other way to fit counts data with standard distribution (gamma, inverse gaussian, etc?)

Upvotes: 2

Views: 2020

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 270348

optim expects its second argument to be a function. Also, the second and third arguments to f are fixed and need to be specified:

optim(c(50, 1, 2), f, x = x, yexp = yexp)

This would also work:

optim(c(50, 1, 2), function(p) f(p, x, yexp))

You could also use nls with default Nelder-Mead algorithm:

nls(yexp ~ a * dgamma(x, sh, scale=sc), start = list(a = 50, sh = 2, sc = 1))

or with plinear in which case no starting value is needed for the first parameter:

nls(c(yexp) ~ dgamma(x, sh, scale=sc), start = list(sh = 2, sc = 1), alg = "plinear")

Upvotes: 4

Related Questions