錯誤運行最大似然估計三參數威布爾cdf

Question

我正在與蒼蠅隨着時間的推移（不規則間隔）在許多夏季（儘管首先我只是想讓一年的工作）累積出現。累積出現遵循S形模式，我想創建一個3參數威布爾累積分佈函數的最大似然估計。我一直試圖在fitdistrplus包中使用的三參數模型不斷給我一個錯誤。我認爲這一定與我的數據結構有關，但我無法弄清楚。顯然，我希望它讀取每個點作爲x（度數天）和y（涌現）值，但它似乎無法讀取兩列。我得到的主要錯誤是「數學函數的非數字參數」或（具有稍微不同的代碼）「數據必須是長度大於1的數字向量」。以下是我的代碼，其中包括在df_dd_em數據框中添加的列，用於累計出現和出現百分比，以防有用。

degree_days <- c(998.08,1039.66,1111.29,1165.89,1236.53,1293.71, 
         1347.66,1387.76,1445.47,1493.44,1553.23,1601.97, 
         1670.28,1737.29,1791.94,1849.20,1920.91,1967.25, 
         2036.64,2091.85,2152.89,2199.13,2199.13,2263.09, 
         2297.94,2352.39,2384.03,2442.44,2541.28,2663.90, 
         2707.36,2773.82,2816.39,2863.94) 
    emergence <- c(0,0,0,1,1,0,2,3,17,10,0,0,0,2,0,3,0,0,1,5,0,0,0,0, 
        0,0,0,0,1,0,0,0,0,0) 
    cum_em <- cumsum(emergence) 
    df_dd_em <- data.frame (degree_days, emergence, cum_em) 
    df_dd_em$percent <- ave(df_dd_em$emergence, FUN = function(df_dd_em) 100*(df_dd_em)/46) 
    df_dd_em$cum_per <- ave(df_dd_em$cum_em, FUN = function(df_dd_em) 100*(df_dd_em)/46) 
    x <- pweibull(df_dd_em[c(1,3)],shape=5) 
    dframe2.mle <- fitdist(x, "weibull",method='mle') 

Answer 1

這是我最好的猜測在你以後：

設置數據：

dd <- data.frame(degree_days=c(998.08,1039.66,1111.29,1165.89,1236.53,1293.71, 
         1347.66,1387.76,1445.47,1493.44,1553.23,1601.97, 
         1670.28,1737.29,1791.94,1849.20,1920.91,1967.25, 
         2036.64,2091.85,2152.89,2199.13,2199.13,2263.09, 
         2297.94,2352.39,2384.03,2442.44,2541.28,2663.90, 
         2707.36,2773.82,2816.39,2863.94), 
       emergence=c(0,0,0,1,1,0,2,3,17,10,0,0,0,2,0,3,0,0,1,5,0,0,0,0, 
       0,0,0,0,1,0,0,0,0,0)) 
dd <- transform(dd,cum_em=cumsum(emergence)) 

實際上，我們要去適應一個「了區間」分佈（即連續度日觀察之間出現的概率：該版本假設第一次觀測是指觀測之前的第一次度日觀測，您可以將其更改爲參考之後的觀測最後一次觀察）。

library(bbmle) 
## y*log(p) allowing for 0/0 occurrences: 
y_log_p <- function(y,p) ifelse(y==0 & p==0,0,y*log(p)) 
NLLfun <- function(scale,shape,x=dd$degree_days,y=dd$emergence) { 
    prob <- pmax(diff(pweibull(c(-Inf,x),  ## or (c(x,Inf)) 
      shape=shape,scale=scale)),1e-6) 
    ## multinomial probability 
    -sum(y_log_p(y,prob)) 
}  
library(bbmle) 

我也許應該使用了更多的東西系統化像矩量法（即匹配的均值與均值和數據的方差威布爾分佈的方差），但我只是砍死了一下週圍找到合理的初始值：

## preliminary look (method of moments would be better) 
scvec <- 10^(seq(0,4,length=101)) 
plot(scvec,sapply(scvec,NLLfun,shape=1)) 

使用parscale到令R知道的參數是非常不同的尺度是非常重要的：

startvals <- list(scale=1000,shape=1) 
m1 <- mle2(NLLfun,start=startvals, 
    control=list(parscale=unlist(startvals))) 

現在與三參數Weibull嘗試（如原先要求） - 只需要的是什麼，我們已經有一個輕微的修改：

library(FAdist) 
NLLfun2 <- function(scale,shape,thres, 
        x=dd$degree_days,y=dd$emergence) { 
    prob <- pmax(diff(pweibull3(c(-Inf,x),shape=shape,scale=scale,thres)), 
       1e-6) 
    ## multinomial probability 
    -sum(y_log_p(y,prob)) 
}  
startvals2 <- list(scale=1000,shape=1,thres=100) 
m2 <- mle2(NLLfun2,start=startvals2, 
    control=list(parscale=unlist(startvals2))) 

看起來像三個參數擬合好得多：

library(emdbook) 
AICtab(m1,m2) 
## dAIC df 
## m2 0.0 3 
## m1 21.7 2 

而這裏的圖形概括：

with(dd,plot(cum_em~degree_days,cex=3)) 
with(as.list(coef(m1)),curve(sum(dd$emergence)* 
          pweibull(x,shape=shape,scale=scale),col=2, 
          add=TRUE)) 
with(as.list(coef(m2)),curve(sum(dd$emergence)* 
          pweibull3(x,shape=shape, 
             scale=scale,thres=thres),col=4, 
          add=TRUE)) 

（也可以做到這一點更優雅與ggplot2 ...）

• 這些似乎並不像奇佳配合，但他們是理智的。 （原則上可以基於每個區間的預期出現次數進行卡方擬合優度檢驗，並說明您已經擬合了三參數模型，儘管這些值可能有點低...）
• 擬合的置信區間有點麻煩;你的選擇是（1）自舉; （2）參數自舉（假設數據的多變量正態分佈的重採樣參數）; （3）三角洲法。
• 使用bbmle::mle2可以很容易做的事情一樣獲得簡介置信區間：

confint(m1) 
##    2.5 %  97.5 % 
## scale 1576.685652 1777.437283 
## shape 4.223867 6.318481 

Answer 2

dd <- data.frame(degree_days=c(998.08,1039.66,1111.29,1165.89,1236.53,1293.71, 
          1347.66,1387.76,1445.47,1493.44,1553.23,1601.97, 
          1670.28,1737.29,1791.94,1849.20,1920.91,1967.25, 
          2036.64,2091.85,2152.89,2199.13,2199.13,2263.09, 
          2297.94,2352.39,2384.03,2442.44,2541.28,2663.90, 
          2707.36,2773.82,2816.39,2863.94), 
      emergence=c(0,0,0,1,1,0,2,3,17,10,0,0,0,2,0,3,0,0,1,5,0,0,0,0, 
         0,0,0,0,1,0,0,0,0,0)) 

dd$cum_em <- cumsum(dd$emergence) 

dd$percent <- ave(dd$emergence, FUN = function(dd) 100*(dd)/46) 

dd$cum_per <- ave(dd$cum_em, FUN = function(dd) 100*(dd)/46) 

dd <- transform(dd) 


#start 3 parameter model 

library(FAdist) 

## y*log(p) allowing for 0/0 occurrences: 

y_log_p <- function(y,p) ifelse(y==0 & p==0,0,y*log(p)) 

NLLfun2 <- function(scale,shape,thres, 
       x=dd$degree_days,y=dd$percent) { 
    prob <- pmax(diff(pweibull3(c(-Inf,x),shape=shape,scale=scale,thres)), 
      1e-6) 
    ## multinomial probability 
    -sum(y_log_p(y,prob)) 
} 

startvals2 <- list(scale=1000,shape=1,thres=100) 

m2 <- mle2(NLLfun2,start=startvals2, 
     control=list(parscale=unlist(startvals2))) 

summary(m2) 

#graphical summary 

windows(5,5) 

with(dd,plot(cum_per~degree_days,cex=3)) 

with(as.list(coef(m2)),curve(sum(dd$percent)* 
          pweibull3(x,shape=shape, 
            scale=scale,thres=thres),col=4, 
         add=TRUE))