截断最大似然估计器

Question

我正在尝试估计截断的新广义泊松林德利分布的参数。 在此分布中，我们使用两个参数“

alpha1, beta1

”，当我运行最优函数时，我应该获得估计器“

alpha1, beta1

”的值。 我在 R 中得到这些值，我得到的“alpha1”值是正确的， 但问题是我得到的“beta1”值不正确。

我尝试在 R 中运行它，但 beta1 的估计不正确。

## Data Generating 
##PMF of New Generalized Poisson Lindley distribution

marginal=function(alpha1, beta1,x1){
  ab=(alpha1^2*(1+alpha1+beta1+beta1*x1))/((alpha1+beta1)*(alpha1+1)^(x1+2))
  return(ab)
}

alpha1=0.2;beta1=3;Tr=3

t1=c();
for(i in 1:10000){
  a=0:1000
  pr1=marginal(alpha1,beta1,a)
  y1=sample(a,1,replace=TRUE,prob=pr1);
  t1[i]=y1;
}
t1
##_________________###_____________________##___________
length(t1)
summary(t1)
hist(t1)
# Mean and Variance of Untruncated NGPLD

#Mean of Untruncated NGPLD 
U1=function(alpha1,beta1){
  U1=(alpha1+2*beta1)/(alpha1*(beta1+alpha1))
  return(sum(U1))
}
U1(alpha1,beta1);mean(t1)

# Variance of Untruncated NGPLD
var1=function(alpha1,beta1){
  var1=(2*beta1^2*(1+alpha1)+alpha1^2*(1+alpha1)+alpha1*beta1*(4+3*alpha1))/(alpha1^2*(alpha1+beta1)^2)
  return(sum(var1))
}
var1(alpha1,beta1);var(t1)

#Truncated data
Tx=t1[t1>Tr];Tx
hist(Tx)
length(Tx)
#Truncated NEw genralized poisson lindley Disrtribution ,where
#Truncated NGPLD=(PMF/1-Cf),where (1-Cf) is a Survival function

#survival Function
x=Tr
cf=((alpha1+beta1)*(1+alpha1)^(x+2)-(alpha1+alpha1^2+2*alpha1*beta1+alpha1*beta1*x+beta1))/((alpha1+beta1)*(1+alpha1)^(x+2))
1-cf

#Truncated NEw genralized poisson lindley Disrtribution 
UT <-function(x1,alpha1,beta1,Tr){
  x=Tr
  cf=((alpha1+beta1)*(1+alpha1)^(x+2)-(alpha1+alpha1^2+2*alpha1*beta1+alpha1*beta1*x+beta1))/((alpha1+beta1)*(1+alpha1)^(x+2))
  numerator =alpha1^2 * (1 + alpha1 + beta1 + beta1 * x1)
  denominator = (alpha1 + beta1) * (alpha1 + 1)^(x1 + 2)
  result=numerator/denominator
  result1=result/(1-cf)
  return(result1)
}
UT(4:100,0.2,3,3)
#sum of probilities of Truncated NGPLD is one
sum(UT(4:100,0.2,3,3))
plot(Tx,UT(Tx,alpha1,beta1,Tr))

#Mean and Varince of Truncated NGPLD

#mean of Truncated NGPLD
f1=function(alpha1,beta1,Tr){
  x1=Tr+1:100
  f1=(x1*alpha1^2*(1+alpha1+beta1+beta1*x1))/((alpha1+beta1)*(alpha1+1)^(x1+2))
  return(sum(f1))
}
X=f1(alpha1,beta1,Tr)/(1-cf)
X;mean(Tx)

#varince of Truncated NGPLD 

f3=function(alpha1,beta1,Tr){
  x1=Tr+1:100
  f3=(x1*(x1-1)*alpha1^2*(1+alpha1+beta1+beta1*x1))/((alpha1+beta1)*(alpha1+1)^(x1+2))
  return(sum(f3))
}

vx=f3(alpha1,beta1,Tr)/(1-cf)+f1(alpha1,beta1,Tr)/(1-cf)-(f1(alpha1,beta1,Tr)/(1-cf))^2
vx;var(Tx)



#likelihood function of Truncated NGPLD
lik=function(par,data){
  x1<-data
  alpha1=par[1];beta1=par[2];
  ll=-sum(log(UT(Tx, alpha1, beta1,Tr)))
  return(ll)}

# MLE Of Truncated NGPLD
pp=optim(par=c(alph1=1,beta1=2),fn=lik,data=Tx,
         control = list(parscale=c(alph1=1,beta1=1)))
#estimated parameter 
ap1=pp$par[1];bt1=pp$par[2];ap1;bt1

Answer 1

刚刚更改了您的

log-likelihood

功能以使用

x1

而不是

Tx

，它对我有用

llik <- function(par, data) {
  x1 <- data
  alpha1 <- par[1]
  beta1 <- par[2]
  ll <- -sum(log(UT(x1, alpha1, beta1,Tr)))  # change Tx to x1
  return(ll)
}

# set.seed(1) before data generation
pp <- optim(par = c(alph1=1,beta1=2), fn = llik, data = Tx,
                   control = list(parscale = c(alph1=1,beta1=1)))
#estimated parameter 
ap1 = pp$par[1]; bt1 = pp$par[2]
ap1;bt1
# alph1 
# 0.20076 
# beta1 
# 2.845234 

这是在数据上绘制的 MLE 拟合参数的分布：

Tx = sort(Tx)
hist(Tx, probability = TRUE, main = 'data & MLE fit', col = 'gray')
lines(Tx,UT(Tx,ap1,bt1,Tr), col='red')
legend("topright", legend=c("data", "MLE"), lty = c(NA, 1), border = c("gray", NA), fill = c("gray",NA), col = c(NA, 'red'), x.intersp=c(-1.5,0.5))

这是了解更多信息的好资源

新的广义泊松-林德利分布