如何使自定义PrWG分布的参数估计结果为正值？

问题描述 #

自定义PrWG分布的对数似然函数后，使用bbmle::mle2进行参数估计时，尽管尝试多组初始值，仍无法得到满足theta>0、lambda>0、gamma>0且0<alpha<1的参数结果，总有至少一个参数估计值为负。

原代码 #

x<-c(36.8, 47.2, 35.6, 36.7, 55.8, 58.7, 42.3, 37.8, 55.4, 45.2, 31.8, 48.3, 45.3, 48.5, 52.8, 45.4, 49.8, 48.2, 54.5, 50.1, 48.4, 44.2, 41.2, 47.2, 39.1, 40.7, 40.3, 41.2, 30.4, 42.8, 38.9, 34, 33.2, 56.8, 52.6, 40.5, 40.6, 45.8, 58.9, 28.7, 37.3, 36.8, 40.2, 58.2, 59.2, 42.8, 46.3, 61.2, 58.4, 38.5, 34.2, 41.3, 42.6, 43.1, 42.3, 54.2, 44.9, 42.8, 47.1, 38.9, 42.8, 29.4, 32.7, 40.1, 33.2, 31.6, 36.2, 33.6, 32.9, 34.5, 33.7, 39.9)
#### PrWG Distribution Log-likelihood ####
LL_PrWG<-function(theta,lambda,gamma,alpha){
  PDF<-(1-alpha)*(theta^2*x^(theta-1)+lambda*gamma*x^(gamma-1)*(1+x^theta))*(1+x^theta)^(-theta-1)*exp(-lambda*x^gamma)/(1-(alpha*(1+x^theta)^(-theta)*exp(-lambda*x^gamma)))^2
  LL<- -sum(log(PDF))
  return(LL)
}
#### PrWG Distribution Optimization ####
library(bbmle) ## calling R package bbmle###
fit<-mle2(LL_PrWG, start=list(theta=0.05,lambda=0.00001,gamma=0.001,alpha=0.0001), data=list(x), method="BFGS")

解决方案：参数变换实现约束 #

mle2默认使用无约束优化算法，直接指定初始值无法强制参数满足约束。通过参数变换将带约束的参数映射到无约束空间，是最可靠的解决方式：

• 对theta、lambda、gamma：使用指数变换par = exp(z)，确保变换后参数恒大于0
• 对alpha：使用logistic变换alpha = plogis(z)（即1/(1+exp(-z))），确保变换后参数在(0,1)区间内

修改后的代码如下：

x<-c(36.8, 47.2, 35.6, 36.7, 55.8, 58.7, 42.3, 37.8, 55.4, 45.2, 31.8, 48.3, 45.3, 48.5, 52.8, 45.4, 49.8, 48.2, 54.5, 50.1, 48.4, 44.2, 41.2, 47.2, 39.1, 40.7, 40.3, 41.2, 30.4, 42.8, 38.9, 34, 33.2, 56.8, 52.6, 40.5, 40.6, 45.8, 58.9, 28.7, 37.3, 36.8, 40.2, 58.2, 59.2, 42.8, 46.3, 61.2, 58.4, 38.5, 34.2, 41.3, 42.6, 43.1, 42.3, 54.2, 44.9, 42.8, 47.1, 38.9, 42.8, 29.4, 32.7, 40.1, 33.2, 31.6, 36.2, 33.6, 32.9, 34.5, 33.7, 39.9)

# 变换后的对数似然函数：优化无约束参数z1-z4
LL_PrWG_transformed <- function(z1, z2, z3, z4) {
  # 将无约束参数转换为原约束参数
  theta <- exp(z1)
  lambda <- exp(z2)
  gamma <- exp(z3)
  alpha <- plogis(z4)
  
  # 原PDF计算逻辑
  PDF <- (1-alpha)*(theta^2*x^(theta-1)+lambda*gamma*x^(gamma-1)*(1+x^theta))*(1+x^theta)^(-theta-1)*exp(-lambda*x^gamma)/(1-(alpha*(1+x^theta)^(-theta)*exp(-lambda*x^gamma)))^2
  
  # 返回负对数似然
  -sum(log(PDF))
}

library(bbmle)
# 初始值对应原初始值的变换结果
start_vals <- list(
  z1 = log(0.05),    # theta=0.05的变换
  z2 = log(0.00001), # lambda=0.00001的变换
  z3 = log(0.001),   # gamma=0.001的变换
  z4 = qlogis(0.0001)# alpha=0.0001的变换
)

# 执行优化
fit_transformed <- mle2(LL_PrWG_transformed, start = start_vals, data = list(x), method = "BFGS")

# 提取并转换回原参数
estimates <- coef(fit_transformed)
original_params <- list(
  theta = exp(estimates["z1"]),
  lambda = exp(estimates["z2"]),
  gamma = exp(estimates["z3"]),
  alpha = plogis(estimates["z4"])
)

print(original_params)

补充说明 #

1. 参数变换的核心是让优化算法在无约束的实数域搜索，避免因参数越界导致的收敛问题
2. 如果变换后仍出现收敛问题，可以尝试调整初始值（比如结合样本统计特征推导，如gamma可参考样本变异系数相关值），或更换优化方法（如method="Nelder-Mead"）
3. 可通过summary(fit_transformed)查看优化的收敛诊断信息，确认结果可靠性

内容的提问来源于stack exchange，提问作者S.A. Osagie