r - 为什么 nls 和 nlsLM 可以正确拟合泊松分布，但无法拟合负二项分布？

Question

我有两个人工生成的概率质量分布，彼此非常相似，只是一个是泊松分布，而另一个是负二项分布，方差比泊松分布稍大。我使用下面的 R 代码示例生成它们，然后尝试使用 nls 或 nlsLM 函数重新估计初始输入参数值:

library(minpack.lm)
library(ggplot2)

# Number of samples (identical for both distributions)
n <- 10000
# Distribution mean (identical for both distributions)
mn <- 10
# Size variable: relevant to negative binomial only; sets the level of
# over-dispersion relative to the Poisson distribution.  Reproduces
# Poisson in the limit that size --> Inf
sz <- 5

# Generate n random samples
psx <- rpois(n, lambda=mn)                    # Poisson
nbx <- rnbinom(n, size=sz, mu=mn)             # negative binomial
# Sort into sample quantiles
psqnt <- unique(sort(psx))                    # Poisson
nbqnt <- unique(sort(nbx))                    # negative binomial
# Generate empirical cdf functions
pscdf <- ecdf(psx)                            # Poisson
pscumdist <- pscdf(psqnt)
nbcdf <- ecdf(nbx)                            # negative binomial
nbcumdist <- nbcdf(nbqnt)
# Place quantiles and cdf into data frame
psdata <- data.frame(q=psqnt, cdf=pscumdist)  # Poisson
nbdata <- data.frame(q=nbqnt, cdf=nbcumdist)  # negative binomial
# Generate estimated starting values that are modified from true values by
# modest amounts
psstart <- list(lambda=0.8*mn)                # Poisson
nbstart <- list(size=0.8*sz, mu=0.8*mn)       # negative binomial

# Plot the sample density functions
pldata <- rbind(data.frame(x=psx, type="Poisson"),
                data.frame(x=nbx, type="Negative Binomial"))
pldata$type <- factor(pldata$type, levels=c("Poisson", "Negative Binomial"))
hst <- ggplot(pldata, aes(x)) +
       geom_histogram(binwidth=1) +
       facet_grid(type ~ .) +
       theme_gray(base_size=18)
print(hst)

# Re-estimate the Poisson distribution parameter, lambda, using either
# nls or nlsLM
print("Starting Poisson fit now...")
#psfit <- nls(cdf ~ ppois(q, lambda), data=psdata, start=psstart, trace=TRUE)
psfit <- nlsLM(cdf ~ ppois(q, lambda), data=psdata, start=psstart, trace=TRUE)
print(coef(psfit))

# Re-estimate the two negative binomial distribution parameters, size and mu,
# using the same technique
print("Starting negative binomial fit now...")
#nbfit <- nls(cdf ~ pnbinom(q, size, mu), data=nbdata, start=nbstart, trace=TRUE)
nbfit <- nlsLM(cdf ~ pnbinom(q, size, mu), data=nbdata, start=nbstart, trace=TRUE)
print(coef(nbfit))

对 ggplot 的调用生成一对直方图，显示两个显然非常相似的离散概率质量分布:

以下是运行 nlsLM 的结果(nls 也给出了非常相似的结果，只是跟踪提供的信息稍少一些):

> source('~/Desktop/nls_error.R')
[1] "Starting Poisson fit now..."
It.    0, RSS =   0.369437, Par. =          8
It.    1, RSS = 0.00130718, Par. =     9.8698
It.    2, RSS = 9.26239e-05, Par. =    9.98602
It.    3, RSS = 9.26083e-05, Par. =     9.9856
It.    4, RSS = 9.26083e-05, Par. =     9.9856
  lambda 
9.985601 
[1] "Starting negative binomial fit now..."
It.    0, RSS =        nan, Par. =          4          8
It.    1, RSS = 2.122e-314, Par. =          4          8
Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model
In addition: Warning messages:
1: In pnbinom(q, size, mu) : NaNs produced
2: In pnbinom(q, size, mu) : NaNs produced
3: In pnbinom(q, size, mu) : NaNs produced
4: In pnbinom(q, size, mu) : NaNs produced
5: In pnbinom(q, size, mu) : NaNs produced
6: In pnbinom(q, size, mu) : NaNs produced

我的问题:我故意将两个示例构建得尽可能相似，那么为什么一个成功而另一个失败呢？

Answer 1

那是因为你被 R 中的默认参数顺序欺骗了。从 pnbinom() 的帮助页面我们看到 pnbinom() 具有以下语法:

pnbinom(q, size, prob, mu, lower.tail = TRUE, log.p = FALSE)

您在调用中向 pnbinom() 提供三个参数

nls(cdf ~ pnbinom(q, size, mu), data=nbdata, start=nbstart, trace=TRUE)

但即使您的参数称为 mu，它也是第三个参数，因此对应于 pnbinom() 中的 prob。由于您使用的是替代公式，因此您需要命名参数以确保其被解释为mu。以下行按您的预期工作

> nbfit <- nls(cdf ~ pnbinom(q, size, mu=mu), data=nbdata, start=nbstart, trace=TRUE)
0.2185854 :  4 8
0.004568844 :  4.069641 9.972202
0.0001207377 :  4.921435 9.961606
3.952388e-05 :  5.068563 9.966108
3.948957e-05 :  5.071698 9.966222
3.948957e-05 :  5.071696 9.966224

如果 nls() 尝试强制将大小设置为负数，则可能会遇到问题。通过对输入参数求幂可以使其更加稳定

> nbfit <- nls(cdf ~ pnbinom(q, exp(size), 
               mu=exp(mu)), data=nbdata, start=nbstart2, trace=TRUE)
0.2971457 :  3.688879 2.079442
0.2622977 :  0.4969337 2.1664490
0.00517649 :  1.408688 2.316948
6.196776e-05 :  1.610170 2.298254
3.948972e-05 :  1.623637 2.299200
3.948957e-05 :  1.623675 2.299202  

其中 nbstart2 与 nbstart 相同，只是记录了启动参数。