如何從R glm（）係數獲得泊松分佈「lambda」（）係數

-1

我的R腳本在下面產生glm（）coeffs。什麼是泊松拉姆達？這應該是〜3.0，因爲這是我用來創建分配的。如何從R glm（）係數獲得泊松分佈「lambda」（）係數

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

來源

2016-07-24 sermolin

你在這裏做什麼??? !!!您使用glm來配合分配？

嗯，這也不是不可能這樣做的，但它是通過這樣做：

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

即，我們有隻攔截迴歸模型擬合數據。

fit$fitted[1] 
# 3.005

這是一樣的：

mean(x) 
# 3.005

來源

2016-07-24 18:59:34

它看起來就像你正在試圖做的泊松適合彙總或分級數據;這不是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因爲離散數據柱狀圖是挑剔的（具有整數分割點往往使事情變得撲朔迷離，因爲你必須要小心右VS左封）

設立分級泊松數據（密度函數與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) 
}

飛度拉姆達（日誌鏈接有助於防止消極的λ值的數值問題/警告）：

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

來源

2016-07-24 19:29:08

如何從R glm（）係數獲得泊松分佈「lambda」（）係數

回答

相關問題