MLE where loglikelihood function is a function of numerical solutions

Hi there,

I'm trying to solve a ML problem where the likelihood function is a function
of two numerical procedures and I'm having some problems figuring out how to
do this.

The log-likelihood function is of the form L(c,psi) = 1/T sum [log (f(c,
psi)) - log(g(c,psi))], where c is a 2xT matrix of data and psi is the
parameter vector. f(c, psi) is the transition density which can be
approximated. The problem is that in order to approximate this we need to
first numerically solve 3 ODEs. Second, numerically solve 2 non-linear
equations in two unknowns wrt the data. The g(c,psi) function is known, but
dependent on the numerical solutions.
I have solved the ODEs using the deSolve package and the 2 non-linear
equations using the BB package, but the results are dependent on the
parameters.

How can I write a program that will maximise this log-likelihood function,
taking into account that the numerical procedures needs to be updated for
each iteration in the maximization procedure?

Any help will be much appreciated.


Kristian

	[[alternative HTML version deleted]]

Hi Kristian

The obvious approach is to treat it like any other MLE problem:  evaluation
of the log-likelihood is done as often as necessary for the  
optimizer you are using: eg a call to optim(psi,LL,...)  where  
LL(psi) evaluates the log likelihood at psi.  There may be  
computational shortcuts that would work if you knew that LL(psi+eps)  
were well approximated by LL(psi), for the values of eps used to  
evaluate numerical derivatives of LL.  Of course, then you might  
need to write your own custom optimizer.

albyn

Quoting Kristian Lind <kristian.langgaard.lind at gmail.com>:

Hi there,

I'm trying to solve a ML problem where the likelihood function is a function
of two numerical procedures and I'm having some problems figuring out how to
do this.

The log-likelihood function is of the form L(c,psi) = 1/T sum [log (f(c,
psi)) - log(g(c,psi))], where c is a 2xT matrix of data and psi is the
parameter vector. f(c, psi) is the transition density which can be
approximated. The problem is that in order to approximate this we need to
first numerically solve 3 ODEs. Second, numerically solve 2 non-linear
equations in two unknowns wrt the data. The g(c,psi) function is known, but
dependent on the numerical solutions.
I have solved the ODEs using the deSolve package and the 2 non-linear
equations using the BB package, but the results are dependent on the
parameters.

How can I write a program that will maximise this log-likelihood function,
taking into account that the numerical procedures needs to be updated for
each iteration in the maximization procedure?

Any help will be much appreciated.


Kristian

	[[alternative HTML version deleted]]

HI Berend, 

Thank you for your reply. 

......
Finally the likelihood function at the end of your code

#Maximum likelihood estimation using mle package
library(stats4)
#defining loglikelighood function
#T <- length(v)
#minuslogLik <- function(x,x2)
#{    f <- rep(NA, length(x))
#        for(i in 1:T)
#        {
#            f[1] <- -1/T*sum(log(transdens(parameters = parameters, x =
c(v[i],v[i+1])))-log(Jac(outmat=outmat, x2=c(v[i],r[i])))
#        }
#    f
# }

How do the arguments of your function x and x2 influence the calculations in
the likelihood function?
 
What I was thinking was that the x and x2 would be the input for the transdens() and Jac() functions, but I guess that is already taken care of when defining them.  

As written now with argument x and x2 not being used in the body of the
function, there is nothing to optimize.
 
That is correct and that is the problem. The likelihood need to be stated as a function of the parameters, here the vector "parameters". Because in the maximum likelihood estimation we want to maximize the value of the likelihood function by changing the parameter values. The likelihood function is a function of the parameters but only through the functions, for example Kristian() is a function of "parameters" which feeds into Bo(), transdens() and Jac(). Do you have any suggestions how to get around this issue?  

Shouldn't f[1] be f[i] because otherwise the question is why are looping
for( i in 1:T)?
But then returning f as a vector seems wrong here. Shouldn't a likelihood
function return a scalar?

The likelihood function should return a scalar. I think the fix could be to make the function calculate the function value at each "i" and then make it return the mean of all the f[i]s.