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 <jo...@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 <jo...@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.l...@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]]
______________________________________________
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.
