R中复杂多项式（右偏函数）的最大似然参数估计答案

【问题标题】：Maximum Likelihood parameter estimation for complex polynomial (right skew function) in RR中复杂多项式（右偏函数）的最大似然参数估计
【发布时间】：2021-03-01 06:00:38
【问题描述】：

我有一个包含 y.size 和 x.number 的数据集。我正在尝试比较模型和自定义模型的线性回归估计的 AIC 值。我能够成功地运行线性回归的估计。自定义模型估计会产生此错误“优化错误（开始，f，方法 = 方法，黑森州 = TRUE，...）：非有限有限差分值 [2]”我是 ML 模型的新手，所以任何帮助将不胜感激。

y.size <- c(2.69,4.1,8.04,3.1,5.27,5.033333333,3.2,7.25,6.29,4.55,6.1,2.65,3.145,3.775,3.46,5.73,5.31,4.425,3.725,4.32,5,3.09,5.25,5.65,3.48,6.1,10,9.666666667,6.06,5.9,2.665,4.32,3.816666667,3.69,5.8,5,3.72,3.045,4.485,3.642857143,5.5,6.333333333,4.75,6,7.466666667,5.03,5.23,4.85,5.59,5.96,5.33,4.92,4.255555556,6.346666667,4.13,6.33,4,7.35,6.35,4.63,5.13,7.4,4.28,4.233333333,4.3125,6.18,4.3,4.47,4.88,4.5,2.96,2.1,3.7,3.62,5.42,3.8,5.5,3.27,3.36,3.266666667,2.265,3.1,2.51,2.51,4.4,2.64,4.38,4.53,2.29,2.87,3.395,3.26,2.77,3.22,4.31,4.73,4.05,3.48,4.8,4.7,3.05,4.21,5.95,4.39,4.55,4.27,4.955,4.65,3.32,3.48,3.828571429,4.69,4.68,3.76,3.91,4,4.41,4.19,4.733333333,4.32,2.83,3.41,4.42,3.47,3.84,4.39)

x.number <- c(69,62,8,80,13,12,2,22,19,49,840,44,31,56,33,58,91,8,15,86,11,69,12,24,32,27,1,4,26,4,28,33,1516,41,20,58,44,29,58,14,3,3,6,3,26,52,26,29,92,30,18,11,27,19,38,78,57,52,17,45,56,7,37,7,14,13,164,76,82,14,273,122,662,434,126,374,1017,522,374,602,164,5,191,243,134,70,23,130,306,516,414,236,172,164,92,53,50,17,22,27,92,48,30,55,28,296,35,12,350,17,22,53,97,62,92,272,242,170,37,220,452,270,392,314,150,232)

require(bbmle)

linreg <- function(a, b, sigma){
  y.pred <- a + b * x.number
  -sum(dnorm(y.size, mean = y.pred, sd = sigma, log=TRUE ))
}

mle2.linreg.model <- mle(linreg, start = list(a = 5, b = -0.01 , sigma = 1))

summary(mle2.linreg.model)
-logLik(mle2.linreg.model)
AIC(mle2.linreg.model)

skewfun <- function(aa,bb, Tmin, Tmax, sigma){
   y.pred <- aa * x.number * (x.number - Tmin) * ((Tmax - x.number )^(1/bb)) + 2
   -sum(dnorm(y.size, mean = y.pred, sd = sigma, log=TRUE ))
   }

mle2.skewfun.model <- mle(skewfun, start = list(aa = 4/10^27, bb = 1/10 , Tmin = 2 , Tmax = 300, sigma = 0.1))

编辑：在尝试了 Simon Woodward 提供的初始答案后，我尝试了使用新参数估计的正确偏斜函数，但我得到了类似的错误：Error in solve.default(oout$hessian) : Lapack routine dgesv: system is exactly singular: U[2,2] = 0。我自定义拟合函数以了解曲线应如何查找此虚拟数据。当我遇到错误时，我使用这些参数来运行 ML。这是原始数据+曲线的样子。下面是生成它的代码

y.size <- c(2.69,4.1,8.04,3.1,5.27,5.033333333,3.2,7.25,6.29,4.55,6.1,2.65,3.145,3.775,3.46,5.73,5.31,4.425,3.725,4.32,5,3.09,5.25,5.65,3.48,6.1,10,9.666666667,6.06,5.9,2.665,4.32,3.816666667,3.69,5.8,5,3.72,3.045,4.485,3.642857143,5.5,6.333333333,4.75,6,7.466666667,5.03,5.23,4.85,5.59,5.96,5.33,4.92,4.255555556,6.346666667,4.13,6.33,4,7.35,6.35,4.63,5.13,7.4,4.28,4.233333333,4.3125,6.18,4.3,4.47,4.88,4.5,2.96,2.1,3.7,3.62,5.42,3.8,5.5,3.27,3.36,3.266666667,2.265,3.1,2.51,2.51,4.4,2.64,4.38,4.53,2.29,2.87,3.395,3.26,2.77,3.22,4.31,4.73,4.05,3.48,4.8,4.7,3.05,4.21,5.95,4.39,4.55,4.27,4.955,4.65,3.32,3.48,3.828571429,4.69,4.68,3.76,3.91,4,4.41,4.19,4.733333333,4.32,2.83,3.41,4.42,3.47,3.84,4.39)

x.number <- c(69,62,8,80,13,12,2,22,19,49,840,44,31,56,33,58,91,8,15,86,11,69,12,24,32,27,1,4,26,4,28,33,1516,41,20,58,44,29,58,14,3,3,6,3,26,52,26,29,92,30,18,11,27,19,38,78,57,52,17,45,56,7,37,7,14,13,164,76,82,14,273,122,662,434,126,374,1017,522,374,602,164,5,191,243,134,70,23,130,306,516,414,236,172,164,92,53,50,17,22,27,92,48,30,55,28,296,35,12,350,17,22,53,97,62,92,272,242,170,37,220,452,270,392,314,150,232)
df <- data.frame(x.number, y.size)
df <- df[df$x.number < 750,]

aa <- 1.25/10^88 
bb <- 1/30 
Tdata <- df$x.number
Tmin <- 1
Tmax <- 750


y.pred <- (aa * df$x.number * (df$x.number - Tmin) * abs(Tmax - df$x.number) ^ (1/bb)) + 3
raw.data <- df$y.size

min = Tdata - Tmin
max = Tmax - Tdata

df1 <- data.frame(df$x.number, min, max, y.pred, raw.data)
library(tidyr)
df.long <- gather(df1, data.type, data.measurement, y.pred:raw.data)

#ggplot(aes(x = Tdata , y = raw.data), data = df1) + geom_point()
ggplot(aes(x = Tdata , y = y.pred), data = df1) + geom_point()

library(ggplot2)

ggplot(aes(x = df.x.number , y = data.measurement, color = data.type), data = df.long) + geom_point()

【问题讨论】：

您需要确保 skewfun 始终给出答案。您需要防止或捕获 bb
您还需要捕捉并处理 Tmax - x.number 可能为负数的情况。
如何绑定参数？例如如何绑定 bb 始终高于 0？最小值和最大值类似
我认为 mle 不提供此功能。因此，您将需要使用不同的包或在您的 skewfun 中强制执行它。理想情况下，您应该这样做，以使对数似然在参数上是连续且可微的（我认为这种粗麻布错误是由于不可微性造成的）。
试试这个：y.pred

标签： r model-fitting log-likelihood

【解决方案1】：

问题来自 skewfun 中的数值错误，因为 x.number 可能大于 Tmax。您需要使 skewfun 对可能的参数值更加稳健。

对于这样的混乱数据，复杂的函数往往不会收敛到好的解决方案。在这种情况下，我更喜欢贝叶斯方法。

y.size <- c(2.69,4.1,8.04,3.1,5.27,5.033333333,3.2,7.25,6.29,4.55,6.1,2.65,3.145,3.775,3.46,5.73,5.31,4.425,3.725,4.32,5,3.09,5.25,5.65,3.48,6.1,10,9.666666667,6.06,5.9,2.665,4.32,3.816666667,3.69,5.8,5,3.72,3.045,4.485,3.642857143,5.5,6.333333333,4.75,6,7.466666667,5.03,5.23,4.85,5.59,5.96,5.33,4.92,4.255555556,6.346666667,4.13,6.33,4,7.35,6.35,4.63,5.13,7.4,4.28,4.233333333,4.3125,6.18,4.3,4.47,4.88,4.5,2.96,2.1,3.7,3.62,5.42,3.8,5.5,3.27,3.36,3.266666667,2.265,3.1,2.51,2.51,4.4,2.64,4.38,4.53,2.29,2.87,3.395,3.26,2.77,3.22,4.31,4.73,4.05,3.48,4.8,4.7,3.05,4.21,5.95,4.39,4.55,4.27,4.955,4.65,3.32,3.48,3.828571429,4.69,4.68,3.76,3.91,4,4.41,4.19,4.733333333,4.32,2.83,3.41,4.42,3.47,3.84,4.39)
x.number <- c(69,62,8,80,13,12,2,22,19,49,840,44,31,56,33,58,91,8,15,86,11,69,12,24,32,27,1,4,26,4,28,33,1516,41,20,58,44,29,58,14,3,3,6,3,26,52,26,29,92,30,18,11,27,19,38,78,57,52,17,45,56,7,37,7,14,13,164,76,82,14,273,122,662,434,126,374,1017,522,374,602,164,5,191,243,134,70,23,130,306,516,414,236,172,164,92,53,50,17,22,27,92,48,30,55,28,296,35,12,350,17,22,53,97,62,92,272,242,170,37,220,452,270,392,314,150,232)
range(x.number)
#> [1]    1 1516

require(bbmle)
#> Loading required package: bbmle
#> Loading required package: stats4

# linear
a = 5; b = -0.01; sigma = 1
linreg <- function(a, b, sigma){
  y.pred <- a + b * x.number
  ll <-   -sum(dnorm(y.size, mean = y.pred, sd = sigma, log=TRUE ))
  # cat(a, b, sigma, ll, "\n") # use this to track convergence
  ll
}

mle2.linreg.model <- mle(linreg, start = list(a = 5, b = -0.01 , sigma = 1))
#> Warning in dnorm(y.size, mean = y.pred, sd = sigma, log = TRUE): NaNs produced

#> Warning in dnorm(y.size, mean = y.pred, sd = sigma, log = TRUE): NaNs produced

summary(mle2.linreg.model)
#> Maximum likelihood estimation
#> 
#> Call:
#> mle(minuslogl = linreg, start = list(a = 5, b = -0.01, sigma = 1))
#> 
#> Coefficients:
#>           Estimate   Std. Error
#> a      4.702090332 0.1400488859
#> b     -0.001529292 0.0005658488
#> sigma  1.339407093 0.0843745782
#> 
#> -2 log L: 431.2136
-logLik(mle2.linreg.model)
#> 'log Lik.' 215.6068 (df=3)
AIC(mle2.linreg.model)
#> [1] 437.2136

a <- mle2.linreg.model@coef["a"]
b <- mle2.linreg.model@coef["b"]
sigma <- mle2.linreg.model@coef["sigma"]
y.pred <- a + b * x.number

plot(x.number, y.size)
lines(x.number, y.pred, col = "red")


# skew function
aa = 4/10^27; bb = 1/10; Tmin = 2; Tmax = 300; ymin = 2; sigma = 0.1
skewfun <- function(aa, bb, Tmin, Tmax, ymin, sigma){
  y.pred <- aa * x.number * pmax(0, x.number - Tmin) * pmax(0, Tmax - x.number) ^ (1/bb) + ymin
  ll <- -sum(dnorm(y.size, mean = y.pred, sd = sigma, log=TRUE ))
  # cat(aa, bb, Tmin, Tmax, sigma, ll, "\n") # use this to track convergence
  ll
}

mle2.skewfun.model <- mle(skewfun, start = list(aa = 4/10^27, bb = 1/10 , Tmin = 2 , Tmax = 300, ymin = 2, sigma = 0.1))
#> Warning in dnorm(y.size, mean = y.pred, sd = sigma, log = TRUE): NaNs produced

#> Warning in dnorm(y.size, mean = y.pred, sd = sigma, log = TRUE): NaNs produced

#> Warning in dnorm(y.size, mean = y.pred, sd = sigma, log = TRUE): NaNs produced

#> Warning in dnorm(y.size, mean = y.pred, sd = sigma, log = TRUE): NaNs produced

summary(mle2.skewfun.model)
#> Warning in sqrt(diag(object@vcov)): NaNs produced
#> Maximum likelihood estimation
#> 
#> Call:
#> mle(minuslogl = skewfun, start = list(aa = 4/10^27, bb = 1/10, 
#>     Tmin = 2, Tmax = 300, ymin = 2, sigma = 0.1))
#> 
#> Coefficients:
#>            Estimate   Std. Error
#> aa    -3.715992e-05 1.739826e-03
#> bb     2.645274e+00 1.322961e+02
#> Tmin   1.266428e+02          NaN
#> Tmax   1.396117e+02 3.574602e+02
#> ymin   4.504795e+00 1.235931e-01
#> sigma  1.377656e+00 8.678290e-02
#> 
#> -2 log L: 438.3107
-logLik(mle2.skewfun.model)
#> 'log Lik.' 219.1554 (df=6)
AIC(mle2.skewfun.model)
#> [1] 450.3107

aa <- mle2.skewfun.model@coef["aa"]
bb <- mle2.skewfun.model@coef["bb"]
Tmin <- mle2.skewfun.model@coef["Tmin"]
Tmax <- mle2.skewfun.model@coef["Tmax"]
ymin <- mle2.skewfun.model@coef["ymin"]
sigma <- mle2.skewfun.model@coef["sigma"]
y.pred <- abs(aa) * x.number * pmax(0, x.number - Tmin) * pmax(0, Tmax - x.number) ^ (1/bb) + ymin

plot(x.number, y.size)
lines(x.number, y.pred, col = "red")

^{由reprex package (v0.3.0) 于 2020 年 11 月 19 日创建}

【讨论】：

我尝试了您建议的具有更好估计值的方法，但使用新估计值的 MLE 函数仍然遇到类似错误。请参阅上面对我的代码的编辑
顺便说一句，我在函数外部设置参数值并在函数内部使用 cat() 的方式使调试变得更加容易。
而且，你的数据很乱，我不希望你会很合适。