I am trying to find MLEs of three positive parameters a, mu and theta, and then the value of a function, saying f1.
f1<-function(para)
{
a<-para[1]
mu<-para[2]
the<-para[3]
return(a*mu/the)
}
Step 1 Suppose we have the following (negative) log likelihood function.
where
x_ij and t_ij are known
loglik<-function(para, data)
{
n.quad<-64
a<-para[1]
mu<-para[2]
the<-para[3]
k<-length(table(data$group))
rule<-glaguerre.quadrature.rules(n.quad, alpha = 0)[[n.quad]]
int.ing.gl<-function(y, x, t)
{
(y^(a-mu-1)/(y+t)^(x+a))*exp(-the/y)
}
int.f<-function(x, t) glaguerre.quadrature(int.ing.gl, lower = 0, upper =
Inf, x=x, t=t, rule = rule, weighted = F)
v.int.f<-Vectorize(int.f)
int<-v.int.f(data$count, data$time)
loglik.value<-lgamma(a+data$count)-lgamma(a)+mu*log(the)-lgamma(mu)+log(int)
log.sum<-sum(loglik.value)
return(-log.sum)
}
Step 2 Let's fix true values and generate data.
### Set ###
library(tolerance)
library(lbfgs3)
a<-2
mu<-0.01
theta<-480
k<-10
f1(c(a, mu, theta))
[1] 5e-04
##### Data Generation #####
set.seed(k+100+floor(a*100)+floor(theta*1000)+floor(mu*1024))
n<-sample(50:150, k) # sample size for each group
X<-rep(0,sum(n))
# Initiate time vector
t<-rep(0, sum(n))
# Initiate the data set
group<-sample(rep(1:k,n)) # Randomly assign the group index
data.pre<-data.frame(X,t,group)
colnames(data.pre)<-c('count','time','group')
data<-data.pre[order(data.pre$group),] # Arrange by group index
# Generate time variable
mut<-runif(k, 50, 350)
for (i in 1:k)
{
data$time[which(data$group==i)]<-ceiling(r2exp(n[i], rate = mut[i], shift = 1))
}
### Generate count variable: Poisson
## First, Generate beta for each group: beta_i
beta<-rgamma(k, shape = mu, rate = theta)
# Generate lambda for each observation
lambda<-0
for (i in 1:k)
{
l<-rgamma(n[i], shape = a, rate = 1/beta[i])
lambda<-c(lambda,l)
}
lambda<-lambda[-1]
data<-data.frame(data,lambda)
data$count<-rpois(length(data$time), data$lambda*data$time) # Generate count variable
Step 3 optimization
head(data)
count time group
0 400 1
0 39 1
0 407 1
0 291 1
0 210 1
0 241 1
start.value<-c(2, 0.01, 100)
fit<-nlminb(start = start.value, loglik, data=data,
lower = c(0, 0, 0), control = list(trace = T))
fit
$par
[1] 1.674672e-02 1.745698e+02 3.848568e+03
$objective
[1] 359.5767
$convergence
[1] 1
$iterations
[1] 40
$evaluations
function gradient
79 128
$message
[1] "false convergence (8)"
One of the possible reasons leading to false convergence is the integral in the step 1. In the loglik function, I used glaguerre.quadrature. However, it failed to give correct result because the integral is converging slowly.
I gave an example to look for some suggestion in the following question
Use the Gauss-Laguerre quadrature to approximate an integral in R
Here, I just provide a complete example. Is there any method I can use to handle this integral?
