CODEWITHSUNDEEP
rmathematical-optimization
User186739
User186739

Reputation: 3

Multi-parameter optimization in R

I'm trying to estimate parameters that will maximize the likelihood of a certain event. My objective function looks like that:

event_prob = function(p1, p2) {
  x = ((1-p1-p2)^4)^67 *
    ((1-p1-p2)^3*p2)^5 *
    ((1-p1-p2)^3*p1)^2 *
    ((1-p1-p2)^2*p1*p2)^3 *
    ((1-p1-p2)^2*p1^2) *
    ((1-p1-p2)*p1^2*p2)^2 *
    (p1^3*p2) *
    (p1^4)
  return(x)
}

In this case, I'm looking for p1 and p2 [0,1] that will maximize this function. I tried using optim() in the following manner:

aaa = optim(c(0,0),event_prob)

but I'm getting an error "Error in fn(par, ...) : argument "p2" is missing, with no default". Am I using optim() wrong? Or is there a different function (package?) I should be using for multi-parameter optimization?

Upvotes: 0

Views: 2231

Answers (2)

dmuir
dmuir

Reputation: 4431

This problem can in fact be solved analytically.

The objective function simplifies to

F(p1,p2) = (1-p1-p2)^299 * p1^19 * p2^11

which is to be maximised over the region

C = { (p1,p2) | 0<=p1, 0<=p2, p1+p2<=1 }

Note that F is 0 if p1=0 or p2 =0 or p1+p2 = 1, while if none of those are true then F is positive. Thus the maximum of F occurs in the interior of C

Taking the log

f(p1,p2) = 299*log(1-p1-p2) + 19*log(p1) + 11*log(p2)

In fact it is as easy to solve the more general problem: maximise f over C where

f( p1,..pN) = b*log( 1-p1-..-pn) + Sum{ a[j]*log(p[j])}
where b and each a[j] is positive and
C = { (p1,..pN) | 0<pj, j=1..N and p1+p2+..pN<1 }

The critical point occurs where all the partial derivatives of f are zero, which is at

-b/(1-p1-..-pn) + a[j]/p[j] = 0 j=1..N

which can be written as

b*p[j] + a[j]*(p1+..p[N]) = a[j] j=1..N

or

M*p = a
where M = b*I + a*Ones', and Ones is a vector with each component 1

The inverse of M is

inv(M) = (1/b)*(I - a*Ones'/(b + Ones'*a))

Thus the unique critical point is

p^ = inv(M)*a
   = a/(b + Sum{i|a[i]})

Since there is a maximum, and only one critical point, the critical point must be the maximum.

Upvotes: 1

Lennart Oelschl&#228;ger
Lennart Oelschl&#228;ger

Reputation: 518

Based on Erwin Kalvelagen's comment: Redefine your function event_prob:

event_prob = function(p) {
  p1 = p[1]
  p2 = p[2]
  x = ((1-p1-p2)^4)^67 *
    ((1-p1-p2)^3*p2)^5 *
    ((1-p1-p2)^3*p1)^2 *
    ((1-p1-p2)^2*p1*p2)^3 *
    ((1-p1-p2)^2*p1^2) *
    ((1-p1-p2)*p1^2*p2)^2 *
    (p1^3*p2) *
    (p1^4)
  return(x)
}

You may want to set limits to ensure that p1 and p2 fulfill your constraints:

optim(c(0.5,0.5),event_prob,method="L-BFGS-B",lower=0,upper=1)

Upvotes: 0

Related Questions