Как получить распределение Пуассона «лямбда» из коэффициентов R glm()

Question

Мой R-скрипт производит коэффициенты glm() ниже. Что такое лямбда Пуассона? Это должно быть ~ 3.0, так как это то, что я использовал для создания дистрибутива.

Call: 
glm(formula = h_counts ~ ., family = poisson(link = log), data = pois_ideal_data) 

Deviance Residuals: 
    Min  1Q Median  3Q  Max 
-22.726 -12.726 -8.624 6.405 18.515 

Coefficients: 
      Estimate Std. Error z value Pr(>|z|)  
(Intercept) 8.222532 0.015100 544.53 <2e-16 *** 
h_mids  -0.363560 0.004393 -82.75 <2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1) 

Null deviance: 11451.0 on 10 degrees of freedom 
Residual deviance: 1975.5 on 9 degrees of freedom 
AIC: 2059 

Number of Fisher Scoring iterations: 5 


random_pois = rpois(10000,3) 
h=hist(random_pois, breaks = 10) 
mean(random_pois) #verifying that the mean is close to 3. 
h_mids = h$mids 
h_counts = h$counts 
pois_ideal_data <- data.frame(h_mids, h_counts) 
pois_ideal_model <- glm(h_counts ~ ., pois_ideal_data, family=poisson(link=log)) 
summary_ideal=summary(pois_ideal_model) 
summary_ideal 

Answer 1

Что вы здесь делаете ??? !!! Вы использовали glm, чтобы соответствовать дистрибутиву ???

Ну, это не невозможно сделать, но это делается с помощью этого:

set.seed(0) 
x <- rpois(10000,3) 
fit <- glm(x ~ 1, family = poisson()) 

т.е. мы вписываемся данные с регрессионной модели перехватывают только.

fit$fitted[1] 
# 3.005 

Это то же самое, как:

mean(x) 
# 3.005 

Answer 2

Похоже, что вы пытаетесь сделать пуассоново впору агрегированном или Binned данных; это не то, что glm. Я быстро взглянул на консервированные способы подгонки дистрибутивов к консервированным данным, но не смог их найти; похоже, более ранние версии пакета bda, возможно, предложили это, но не сейчас.

У корня, что вам нужно сделать, это настроить отрицательную функцию логарифмического правдоподобия, которая вычисляет (# counts)*prob(count|lambda) и сводит к минимуму ее с помощью optim(); решение дано ниже с использованием bbmle пакета является немного более сложным, авансовые, но дает дополнительные преимущества, как легко вычислительными доверительные интервалы и т.д ..

Настройка данных:

set.seed(101) 
random_pois <- rpois(10000,3) 
tt <- table(random_pois) 
dd <- data.frame(counts=unname(c(tt)), 
       val=as.numeric(names(tt))) 

Здесь я использую table а не потому, что hist гистограммы на дискретных данных суетливы (с целым часто делает сочленение вещи запутанной, потому что вы должны быть осторожны право- против левого закрытия)

Настройки функции плотности для Binned данных Пуассона (для работы с bbmle " формы ula, первый аргумент должен быть вызван x, и он должен иметь аргумент log).

dpoisbin <- function(x,val,lambda,log=FALSE) { 
    probs <- dpois(val,lambda,log=TRUE) 
    r <- sum(x*probs) 
    if (log) r else exp(r) 
} 

Fit лямбда (срубы ссылка помогает предотвратить численные проблемы/предупреждение от отрицательных значений лямбды):

library(bbmle) 
m1 <- mle2(counts~dpoisbin(val,exp(loglambda)), 
     data=dd, 
     start=list(loglambda=0)) 
all.equal(unname(exp(coef(m1))),mean(random_pois),tol=1e-6) ## TRUE 
exp(confint(m1)) 
## 2.5 % 97.5 % 
## 2.972047 3.040009