拟合三角分布

Question

我有以下数据

dat<-c(16.254884, 14.077510, 12.851675, 19.152597, 11.511230, 
   16.122911, 16.099962,  9.670949, 12.523661, 15.257432, 13.603848,
   14.118873, 12.632340, 15.413753,  5.426383, 11.369880, 12.895920, 
   13.635134, 15.118388,13.154107, 8.913164, 17.302810, 14.968054, 
   16.200151, 16.068944, 18.571952, 15.247535, 15.018281)

我正在使用此代码查找模式：

Mode_fc <- function(x) {
     ux <- unique(x)
     ux[which.max(tabulate(match(x, ux)))]
 }

使用MyParam，我可以获得最小值、最大值和众数

MyParam <- c(min= min(dat), max= max(dat), mode= Mode_fc(dat))

当我在下面的代码中输入这些值时 fitdist 按预期工作

fitdist(dat, "triang", start = list(min=5.4, max=19.2, mode=16.3))

但是，当我尝试读取 MyParam 时，我遇到了各种错误

fitdist(dat, "triang", 
 start = list(min=MyParam[[1]], max=MyParam[[2]], mode=MyParam[[3]]))

我知道问题出在 optim()，但我还没弄清楚如何解决这个问题。任何建议表示赞赏！

Answer 1

你的问题（诚然相当微妙）是如果任何数据在三角分布下数据的可能性为零（因此对数似然为负无限），或在的边界上分布。插图：

library(fitdistrplus)
library(mc2d)  ## needed for dtriang

试穿（如上例所示）：

L1 <- list(min=5.4, max=19.2, mode=16.3)
fitdist(dat, "triang", start = L1)  ## works
L2 <- list(min=MyParam[[1]], max=MyParam[[2]], mode=MyParam[[3]])
fitdist(dat, "triang",  start = L2)  ## fails

让我们稍微分解一下，看看每组参数的实际对数似然是什么：

do.call(dtriang,c(list(x=dat,log=TRUE),L1))
##  [1] -1.935669 -2.159550 -2.311845 -6.045302 -2.510156 -1.947902 -1.950044
##  [8] -2.868448 -2.356862 -2.032059 -2.215681 -2.154794 -2.341722 -2.016325
## [15] -7.955320 -2.533557 -2.305925 -2.211875 -2.046264 -2.272062 -3.063767
## [22] -2.355858 -2.061854 -1.940724 -1.952947 -3.461371 -2.033063 -2.056619

所有有限值。

(test2 <- do.call(dtriang,c(list(x=dat,log=TRUE),L2)))
##  [1] -1.926160 -2.150652 -2.303450      -Inf -2.502540 -1.938423 -1.940570
##  [8] -2.862702 -2.348631 -2.022796 -2.206960 -2.145882 -2.333434 -2.007021
## [15]      -Inf -2.526044 -2.297509 -2.203141 -2.037041 -2.263528 -3.059363
## [22] -2.375012 -2.052673 -1.931228 -1.943481 -3.533698 -2.023803 -2.047423

两个无穷大值，分别对应最小值和最大值。

which(!is.finite(test2))  ## 4 15
which.min(test2)  ## 4
which.max(test2)  ## 5

我们可以通过将观察值的最小值向下调整一点并将最大值向上调整一点来轻松解决此问题：

eps <- 0.100
L3 <- list(min=MyParam[[1]]-eps, max=MyParam[[2]]+eps, mode=MyParam[[3]])
fitdist(dat, "triang",  start = L3)

这很好用。

拟合三角分布

Fit triangular distribution

r

distribution

fitdistrplus