From: Kristian Lind <kristian.langgaard.lind_at_gmail.com>

Date: Wed, 13 Apr 2011 12:00:03 +0200

f

}

#declaring state variables i.e. the functions and their intial conditions

state <- c(b_1 = 0,

#delclaring function and system of equations.

Kristian <- function(t, state, parameters) {

})

}

# time making a sequence from t to T evaluated at each delta seq(t, T, by =

delta)

times <- seq(0, 10, by = 0.5)

}

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)

{

#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)

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

Date: Wed, 13 Apr 2011 12:00:03 +0200

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_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 >>> PLEASE do read the posting guide >>> 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 >> PLEASE do read the posting guide >> 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 > 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]] ______________________________________________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.
*