# Re: [R] MLE where loglikelihood function is a function of numerical solutions

From: Kristian Lind <kristian.langgaard.lind_at_gmail.com>
Date: Wed, 13 Apr 2011 12:00:03 +0200

Albyn and others,

In order to be more specific I've constructed my program. I know it's long and in some places quite messy. It works until the last part where the log-likelihood function has to be defined and maximized wrt the parameters.

The log-likelihood has the form L = 1/T*sum_t=0^T[log(transdens(v_t, r_t))-log(Jac(r_t,v_t))], where the functions transdens(v_t, r_t) and Jac(r_t,v_t) are defined below.

The problem remains the same, how do I construct the log-likelihood function such that the numerical procedures employed are updated in the maximization procedure?

I was thinking something along the lines of

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
}

Kristian

```----------------------------------Program starts
here------------------------------

# Solving ODEs

parameters <- c(K_vv    <-     0.0047,
K_rv    <-    -0.0268,
K_rr    <-    0.3384,
theta_v <- 50,
theta_r <-    5.68,
Sigma_rv<-    0.0436,
Sigma_rr<-    0.1145,
lambda_v<-    0,
lambda_r<-    -0.0764,
k_v        <-    K_vv*theta_v,
k_r        <-     K_rv*theta_v+K_rr*theta_r,
alpha_r <-     0,
B_rv    <-     (Sigma_rr+Sigma_rv)^2)
```

#declaring state variables i.e. the functions and their intial conditions
state <- c(b_1 = 0,
```            b_2 = 0,
a = 0)
```

#delclaring function and system of equations.
Kristian <- function(t, state, parameters) {

with(as.list(c(state, parameters)),
{

db_1 =
-((K_vv+lambda_v)*b_1+(K_rv+Sigma_rv*lambda_v+Sigma_rr*lambda_r)*b_2+0.5*(b_1)^2+Sigma_rv*b_1*b_2+0.5*((Sigma_rv)^2+(Sigma_rr)^2)*(b_2)^2 )

```        db_2 = -K_rr*b_2+1
da = K_vv*theta_v*b_1+(K_rv*theta_v+K_rr*theta_r)*b_2
list(c(db_1, db_2, da))
```

})
}
# time making a sequence from t to T evaluated at each delta seq(t, T, by =
delta)
times <- seq(0, 10, by = 0.5)

#solving the model using the deSolve function ode
library(deSolve)
outmat <- ode(y = state, times = times, func = Kristian, parms = parameters)
#print(outmat)
#Solving 2 equations in 2 unknowns using nleslv
# simulating data for testing

c2 <- matrix(data = rnorm(20, mean =0, sd =1)+2, nrow = 20, ncol =1) c10 <- matrix(data = rnorm(20, mean =0, sd =1)+5, nrow = 20, ncol =1) besselinput <- matrix(data = 0, nrow =20, ncol = 1) vncChi <- matrix(data = 0, nrow =20, ncol = 1) v <- matrix(data = 0, nrow =20, ncol = 1) r <- matrix(data = 0, nrow =20, ncol = 1) for(i in 1:20)
{

Bo <- function(x, s, outmat)
{

```        f <- rep(NA, length(x))
z <- - outmat[,2:4] %*% c(x[2],x[1],1)
f[1] <- (1-exp(z[4]))/sum(exp(z[1:4])) - s[1]
f[2] <- (1-exp(z[20]))/sum(exp(z)) - s[2]
f
```

}
s <- c(c2[i],c10[i] )
p <- c(50, 5)

library(nleqslv)
ans.nlq <- nleqslv(x=p, fn=Bo, s=s, outmat=outmat, control = list(matxit =100000))
v[i] <- ans.nlq\$x[1]
r[i] <- ans.nlq\$x[2]
#print(ans.nlq\$termcd)
#ans.nlq\$fvec

}
#calculating transition density as a function of parameters
transdens <- function(parameters, x)
{

delta <-1

```    c <- 2*K_vv*(1-exp(-K_vv*delta))^(-1) #2.004704
q <- 2*k_v-1 #0.00959666
f <- rep(NA, 1)
```

#besselinput[2] <- 2*c*(x[2]*x[1]*exp(-K_vv*delta))^(0.5)     f[1]<-
c*exp(c*(x[2]+exp(-K_vv*delta)*x[1]))*(x[2]/(exp(-K_vv*delta)*x[1]))^(q/2)*besselI(2*c*(x[2]*x[1]*exp(-K_vv*delta))^(0.5), q, expon.scaled = FALSE)

f
}
vncChi <- transdens(parameters = parameters, x=c(v[1],v[2])) print(vncChi)
#calculating Determinant of Jacobian as a function of outmat.

Jac <- function(outmat, x2){
f <- rep(NA, 1)

```    y <- - outmat[,2:4] %*% c(x2[1],x2[2],1)
w1 <- outmat[,2]*exp(-outmat[,2:4] %*% c(x2[1],x2[2],1))
w2 <- outmat[,3]*exp(-outmat[,2:4] %*% c(x2[1],x2[2],1))
```
f[1] <- abs((outmat[5,2]*exp(y[4])/(sum(exp(y[1:4])))+
```    (1-exp(y[4]))*sum(w1[1:4])/(sum(exp(y[1:4])))^2)*(outmat[21,3]
*exp(y[21])/(sum(exp(y)))+(1-exp(y[21]))*sum(w2)/
(sum(exp(y)))^2)-(outmat[21,2]*exp(y[21])/(sum(exp(y)))
```
+(1-exp(y[21]))*sum(w1)/(sum(exp(y)))^2)*(outmat[5,3]     *exp(y[4])/(sum(exp(y[1:4])))+(1-exp(y[4]))*sum(w2[1:4])     /(sum(exp(y[1:4])))^2))

f
}

#----------------------------------------------Program works as it is

```supposed to until here----------------------------

```

#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

}

2011/4/10 Albyn Jones <jones_at_reed.edu>

```> to clarify: by  "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.
> "
> I mean the derivatives of LL(psi+eps) are close to the derivatives of
> LL(psi),
> and perhaps you would want the hessian to be close as well.
>
> albyn
>
>
> Quoting Albyn Jones <jones_at_reed.edu>:
>
>  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]]
>>>
>>> ______________________________________________
>>> R-help_at_r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>>
>> ______________________________________________
>> R-help_at_r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>>
> ______________________________________________
> R-help_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

______________________________________________
```
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 13 Apr 2011 - 10:08:42 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 28 Apr 2011 - 11:00:32 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.