2016-04-06 9 views
2

После генерации отрицательных биномиальных данных с набором задач, равным 0,001, я получаю это число обратно из glm.nb(), но только путем обмана.Как я могу получить параметр «prob» из glm.nb()?

library(MASS) 
counts<-data.frame(as.matrix(rnbinom(10000, prob = .007, size = 247))) 
names(counts)<-"y" 

head(counts) 

fitted_model<-glm.nb(y ~ 1, data = counts, link="identity") 

#Theta is the shape parameter of the negative binomial distribution. So this is "r". 
r<-theta.ml(fitted_model$y, fitted(fitted_model))[1]  
# the parameter r is referred to as the “dispersion parameter” or “shape parameter” 

mu<-coef(fitted_model) #This is the mean 

# mu=prob*r/(1-prob) according to https://en.wikipedia.org/wiki/Negative_binomial_distribution 
# so prob = 1/(r + mu) ? 
1/(r + mu) # Wrong! This isn't the prob I used to generate th data! 
r/(r + mu) # Right! But why does this get me the correct value of prob? 

#This has hints: http://www.wright.edu/~thaddeus.tarpey/ES714glm.pdf 

Я не хочу обманывать, чтобы получить значение «prob» из установленной модели. Может ли кто-нибудь объяснить, почему r/(r + mu) = prob?

ответ

2

Если сравнить определение Википедии

C(k+r-1,k) (1-p)^r p^k 

с определением, данным в ?NegBinomial

Gamma(x+n)/(Gamma(n) x!) p^n (1-p)^x 

вы увидите, что роли p и 1-p включаются; если мы определим NB как «вероятность n успехов, возникающих до одного отказа», то Википедия определяет p как вероятность «отказа», в то время как R определяет p как вероятность «успеха». Я получаю правильный результат от r/(r+mu), а не от mu/(r+mu) ...

+0

Теперь я узнал две вещи: решение этой проблемы и правило, что всегда следует сравнивать формулы pdf в обеих системах. – rwinkel2000