Интеграция и ложная сходимость оптимизации в R - PullRequest
Новая
       48

Интеграция и ложная сходимость оптимизации в R

0 голосов
/ 22 мая 2018

Я пытаюсь найти MLE из трех положительных параметров a, mu и theta, а затем значение функции, говоря f1. 

f1<-function(para)
{
 a<-para[1]
 mu<-para[2]
 the<-para[3]
 return(a*mu/the)
}

Шаг 1 Предположим, у нас есть следующая (отрицательная) логарифмическая функция правдоподобия.

loglikelihood

, где

integral

x_ij и t_ij известны 

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

Шаг 2 Давайте исправим истинные значения и сгенерируем данные. 

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

Шаг 3 оптимизация 

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

Одной из возможных причин, приводящих к ложной сходимости, является интеграл в шаге 1 .В функции loglik я использовал glaguerre.quadrature.Однако это не дало правильного результата, потому что интеграл медленно сходится.Я привел пример для поиска некоторого предложения в следующем вопросе

Используйте квадратуру Гаусса-Лагерра для аппроксимации интеграла в R

Здесь я просто предоставляю полноепример.Можно ли использовать какой-либо метод для обработки этого интеграла?

...