查找具有最低组件值的变量名称的快速方法

Question

我有一个拟合分布的函数和 returning 一个由分布名称、均值、sd 等组成的向量。我正在测试几个发行版，但我不能依赖 gofstat() 因为当要考虑的零太多时它会发疯。

因此我必须手动比较几个变量的 AIC，决定哪些实际上是 "fitdist" class 和 return 具有最低 AIC 的变量名称。一旦有了它，我就会计算平均值、sd 等和 return.

代码目前如下所示：

library(fitdistrplus)

fit_distr <- function(data){

fe <- tryCatch(fitdist(data, "exp"), error = function(e) FALSE )
flogis <- tryCatch(fitdist(data, "logis"), error = function(e) FALSE )
fn <- tryCatch(fitdist(data, "norm"), error = function(e) FALSE)
fp <- tryCatch(fitdist(data, "pois"), error = function(e) FALSE)
fg <- tryCatch(fitdist(data, "gamma"), error = function(e) FALSE)

classFitDist <- c(class(fe), class(flogis), class(fn), class(fp),class(fg))

distributions <- classFitDist == "fitdist"


AIC <- data.frame()
if(class(fe)=="fitdist") {AIC[1,ncol(AIC)+1] <- fe$aic}
if(class(flogis)=="fitdist") {AIC[1,ncol(AIC)+1] <- flogis$aic}
if(class(fn)=="fitdist") {AIC[1,ncol(AIC)+1] <- fn$aic}
if(class(fp)=="fitdist") {AIC[1,ncol(AIC)+1] <- fp$aic}
if(class(fg)=="fitdist") {AIC[1,ncol(AIC)+1] <- fg$aic}

names(AIC) <- c("exp", "logis", "norm", "pois", "gamma")[distributions]

fit <- names(AIC[which.min(AIC)])


mean <- switch (fit,
                          exp = 1/fe$estimate[[1]],
                          logis = flogis$estimate[[1]],
                          norm = fn$estimate[[1]],
                          pois = fp$estimate[[1]],
                          gamma = fg$estimate[[1]]/fg$estimate[[2]]
          )

sd <- switch (fit,
                        exp = mean,
                        logis = (flogis$estimate[[2]]*pi)/sqrt(3),
                        norm = fn$estimate[[2]],
                        pois = sqrt(mean),
                        gamma = sqrt(fg$estimate[[1]]/(fg$estimate[[2]]^2))
          )

return(c(fit,mean,sd))
}

它可以工作，但是考虑数千个样本非常慢。我欢迎任何关于如何优化它并使其 'cleaner' 和更快的建议。

顺便说一句，这是我以前的，但是就像我提到的那样 - 样本包含太多零，它很合适（双关语无意！）

goodnessoffit <- gofstat(list(fe, flogis, fn, fp, fg)[distributions],  fitnames = c("exp", "logis", "norm", "pois","gamma")[distributions])
fit <- names(which(goodnessoffit$aic == min(goodnessoffit$aic)))

Error in ans[!test & ok] <- rep(no, length.out = length(ans))[!test & : replacement has length zero

Answer 1

这种方法的问题是 fitdist 效率低下。你需要通过编写更好的算法来想出更快的方法来找到 AIC。一种方法是拟合 glm。

AIC.fitdist <- function(x, ...) x$aic

x <- rnorm(100, mean=20)

AIC(fitdist(x, 'norm'))
AIC(glm(x ~ 1 , family=gaussian)) ## same

AIC(fitdist(x, 'gamma'))
AIC(glm(x ~ 1 , family=Gamma)) ## same

一些分析显示 fitdist 与 glm 具有相同的计算时间。这对 fitdist 来说是个坏消息，因为 glm 只是 glm.fit 的臃肿包装。使用 glm.fit 可以为您节省大量时间。最后，如果您真的不得不减少模型的时间（数百万，而不是数千），您可以使用

的一步估计器

> benchmark(
+     fitdist(x, 'gamma'),
+     glm(x ~ 1, family=Gamma),
+     glm.fit(rep(1, length(x)), x, family=Gamma()),
+     glm.fit(rep(1, length(x)), x, family=Gamma(), control = glm.control(maxit=1))
+ )
                                                                               test replications elapsed relative user.self sys.self user.child
1                                                               fitdist(x, "gamma")          100    0.42    7.000      0.42        0         NA
2                                                        glm(x ~ 1, family = Gamma)          100    0.17    2.833      0.17        0         NA
3                                   glm.fit(rep(1, length(x)), x, family = Gamma())          100    0.06    1.000      0.07        0         NA
4 glm.fit(rep(1, length(x)), x, family = Gamma(), control = glm.control(maxit = 1))          100    0.06    1.000      0.06        0         NA
  sys.child
1        NA
2        NA
3        NA
4        NA

aic 是 glm.fit 输出中的存储对象。

生存包中的survreg可以进行指数分布拟合：survreg(rep(1,100), x, dist='exponential).

最后，由于这些都是正则指数族，您可以使用充分的统计数据来得出概率分布。例如：

normaic <- function(x) {
  4 - 2*sum(dnorm(x, mean(x), sd(x), log=T))
}

> benchmark(normaic(x), glm.fit(rep(1, 100), x)$aic)
                         test replications elapsed relative user.self sys.self user.child sys.child
2 glm.fit(rep(1, 100), x)$aic          100    0.04       NA      0.05        0         NA        NA
1                  normaic(x)          100    0.00       NA      0.00        0         NA        NA

查找具有最低组件值的变量名称的快速方法

Quick way of finding a name of a variable with lowest component value

r

optimization

goodness-of-fit

fitdistrplus