# R fitting Poisson distribution with weighting

I am wanting to fit a Poisson distribution but I am wanting to weight some of the values as more important. Is there a way to do this in R with fitdistr() or some equivalent function? Here is what I have at the moment:

randoms <- rpois(15,10)
weighting <- seq(1, 100, by=1)
fit <- fitdistr(randoms, "poisson")


I would like to use the 'weighting' data as a way to emphasize the observation at the end of the data set.

• Other replies you received on Stack Overflow were merged with the ones here. Please, don't cross-post in the future. – chl Nov 2 '12 at 17:01
• Could you then please help me decide where to post this question? It is programming in R but it is also related to statistics which is what this website deals with. I know that cross-posting is frowned upon but I got 2 different answers which are both very valuable. I love these websites but I'm not always sure where to post questions. Is there a decision tree that can be referenced? Thanks for the help. – Wallhood Nov 3 '12 at 15:06
• Ah, I followed the link on cross-posting in your comment and cleared it up for me. Thanks for the help! – Wallhood Nov 3 '12 at 16:26
• No problem. One solution about those two valuable answers (you accepted both) is to mention that you like them in a comment (beneath each answer) or to give an equal amount of your reputation via bounties. – chl Nov 3 '12 at 21:05

Note that 'fitdistr does nothing but maximum likelihood estimation. That is to say, you can do it by yourself by writing down the likelihood. Below is an example for the poisson distribution in R. It can be adapted to upweight/downweight the contribution to the likelihood of each data.

density (as in R) $$f(x; \lambda) = \lambda^x \frac{\exp(-\lambda)}{x!}, \qquad \lambda > 0$$

likelihood $$L(\lambda; \mathbf{x}) = \prod_{i=1}^n \left\{ \lambda^{x_i} \frac{\exp(-\lambda)}{x_i!} \right\}$$

log-likelihood

$$\ell(\lambda; \mathbf{x}) = \sum_{i=1}^n \left\{ x_i \log(\lambda) - \lambda - \log(x_i!) \right\}$$

2nd derivatie of $\ell$:

$$\frac{d^2 \ell}{d\lambda^2}(\lambda; \mathbf{x}) = \sum_{i=1}^n - \frac{x_i}{\lambda^2} = -\frac{1}{\lambda^2} n\bar{x}$$

#------data------
set.seed(730)
sample <- rpois(1000, 10)
#----------------

################################################################################
# Using 'fitdistr'                                                             #
################################################################################

library(MASS)
fitdistr(x=sample, densfun="Poisson")
lambda
10.1240000
( 0.1006181)

################################################################################
# writing down the log-likelihood explicitly                                   #
################################################################################

#------minus log-likelihood------
mloglik <- function(lambda2, sample) #lambda2 = log(lambda) in (-\infty, \infty)
{
- sum(sample * lambda2 - exp(lambda2) - log(factorial(sample)))
}
#--------------------------------

#------optimisation------
res <- nlm(f=mloglik, p=1, sample=sample)
#------------------------

#------recover lambda------
lambda <- exp(res$estimate) round(lambda, 7) [1] 10.12399 #-------------------------- #------standard error------ #square root of negative inverse second derivative of the log-likelihood se <- lambda / sqrt(length(sample) * mean(sample)) round(se, 7) [1] 0.100618 #--------------------------  If you're looking to do this just for the poisson distribution (where you're just estimating the point estimate of the one parameter lambda), you could use the glm function: fit2 <- glm(randoms ~ 1, family = poisson(link = "log"), weights = weighting[1:15])  As a check, fit1 <- glm(randoms ~ 1, family = poisson(link = "log")) all.equal(unname(exp(coef(fit1))), unname(fit$estimate))
# [1] TRUE


And another check:

randoms2 <- rep(randoms, weighting[1:15])
fit3 <- glm(randoms2 ~ 1, family = poisson(link = "log"))

all.equal(exp(coef(fit3)), exp(coef(fit2)))
# [1] TRUE


If you look at the likelihood, the only bit of it that depends on both lambda and the x's is product lambda^x_i, which is lambda^{sum x_i}. Thus the sum of the x's, or equivalently the sample mean, is a sufficient statistic for lambda.

What that means is that you "fit a Poisson" by estimating the mean, using the sample mean. (And then you look at how well the Poisson with that mean fits your data.)

To go back to the OP's problem of weighting, all that means is that you calculate a weighted mean, and use that to estimate lambda with.