Albyn and others,

Thank you for your replies.
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
 }

Thank you in advance.

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)
# loading nleqslv package
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@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@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@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@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.
>>>
>>>
>>>
>> ______________________________________________
>> R-help@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.
>>
>>
>>
> ______________________________________________
> R-help@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.
>

	[[alternative HTML version deleted]]

