如何从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”（）系数

回答

相关问题