2

我正在研究在许多夏天（尽管首先我只是想让一年的工作）随着时间的推移（以不规则的时间间隔拍摄）苍蝇的累积出现。累积出现遵循 sigmoid 模式，我想创建 3 参数 Weibull 累积分布函数的最大似然估计。我一直尝试在fitdistrplus包中使用的三参数模型不断给我一个错误。我认为这一定与我的数据结构有关，但我无法弄清楚。显然，我希望它将每个点读为x（学位日）和y(emergence) 值，但似乎无法读取两列。我得到的主要错误是“数学函数的非数字参数”或（代码略有不同）“数据必须是长度大于 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')
4

2 回答 2

3

这是我对你所追求的最好的猜测：

设置数据：

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)

我可能应该使用更系统的方法，例如矩方法（即，将 Weibull 分布的均值和方差与数据的均值和方差相匹配），但我只是稍微修改了一下以找到合理的起始值：

## 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) 增量法。
  • Usingbbmle::mle2可以很容易地做一些事情，比如获取配置文件置信区间：

 confint(m1)
 ##             2.5 %      97.5 %
 ## scale 1576.685652 1777.437283
 ## shape    4.223867    6.318481
于 2014-07-22T22:01:42.017 回答
0 
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))

在此处输入图像描述

于 2014-07-23T23:55:09.683 回答