Почему nls и nlsLM работают правильно для подгонки распределения Пуассона, но не работают для отрицательного бинома? - PullRequest
Новая
       21

Почему nls и nlsLM работают правильно для подгонки распределения Пуассона, но не работают для отрицательного бинома?

0 голосов
/ 11 января 2019

У меня есть два искусственно сгенерированных распределения вероятности по массе, очень похожих друг на друга, за исключением того, что одно является распределением Пуассона, а другое - отрицательным биномом, с немного большей дисперсией, чем Пуассон. Я генерирую их, используя приведенный ниже пример кода 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 генерируется пара гистограмм, показывающих два распределения массы дискретной вероятности, которые, очевидно, очень похожи: Poisson distribution and negative binomial distribution

Вот результаты запуска 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

Мой вопрос: Я сознательно построил два примера, чтобы они были как можно более похожими, так почему же один из них успешен, а другой - нет?

1 Ответ

0 голосов
/ 11 января 2019

Это потому, что вас обманывают в порядке аргументов по умолчанию в 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, это третий аргумент и, следовательно, соответствует prob в pnbinom(). Поскольку вы используете альтернативную формулировку, вам нужно назвать аргумент, чтобы убедиться, что он интерпретируется как 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, за исключением того, что регистрируются начальные параметры.

...