[R] Is there dpois equivalent for zero-inflated Poisson?
Michael Friendly
friendly at yorku.ca
Thu Mar 24 14:06:06 CET 2016
Neat!
It would be nice to complete dzipois() with the corresponding rzipois()
and pzipois() functions. I would have found these useful in my new book,
http://ddar.datavis.ca
-Michael
On 3/23/2016 11:27 AM, Martin Maechler wrote:
>>>>>> Thierry Onkelinx <thierry.onkelinx at inbo.be>
>>>>>> on Tue, 22 Mar 2016 13:58:09 +0100 writes:
>
> > dpois(0, lambda) == e^(-lambda)
> > The wikipedia formula is
> > ifelse(x == 0, zero + dpois(0, lambda) * (1-zero), dpois(x, lambda) *
> > (1-zero))
>
> > or
>
> > ifelse(x == 0, zero + dpois(x, lambda) * (1-zero), dpois(x, lambda) *
> > (1-zero))
>
> > so we can move the dpois() out of the ifelse()
>
> > ifelse(x == 0, zero, 0) + dpois(x, lambda) * (1-zero)
>
> Nice! Thank you, Thierry.
>
> Even nicer for symmetry reasons (and much faster) is *not* to use ifelse(),
> so we'd get
>
> dzipois <- function(x, lambda, zero)
> (x == 0) * zero + dpois(x, lambda) * (1-zero)
>
> However, numerically correctly adding the 'log = FALSE'
> argument which all good "density" functions have in R --
> ((and which you *should* use as log = TRUE for the
> log-likelihood if you want to become more professional))
> is a bit tricky.
>
> The x == 0 case at least needs care for large lambda and/or
> small 'zero' (= pi).
> All this is related on how to accurately compute f(L) = log(1 - exp(- L))
> which I call log1mexp(L) in my still-not-submitted paper
> available as one of the Rmpfr vignettes at
> https://cloud.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
>
> BTW: There are at least 3-4 R packages which deal with zero-inflated
> poisson in some ways, notably the recommended 'mgcv'
> package which comes bundled with every R.
>
>
> Martin Maechler,
> ETH ZUrich
>
> > ir. Thierry Onkelinx
> > Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> > Forest
> > team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> > Kliniekstraat 25
> > 1070 Anderlecht
> > Belgium
>
> > To call in the statistician after the experiment is done may be no more
> > than asking him to perform a post-mortem examination: he may be able to say
> > what the experiment died of. ~ Sir Ronald Aylmer Fisher
> > The plural of anecdote is not data. ~ Roger Brinner
> > The combination of some data and an aching desire for an answer does not
> > ensure that a reasonable answer can be extracted from a given body of data.
> > ~ John Tukey
>
> > 2016-03-22 13:50 GMT+01:00 Matti Viljamaa <mviljamaa at kapsi.fi>:
>
> >> And why is the first term of ifelse(x == 0, zero, 0) + dpois(x, lambda) /
> >> (1 - zero)
> >>
> >> ifelse(x == 0, zero, 0)
> >>
> >> rather than something corresponding to
> >>
> >> zero+(1-zero)e^{-lambda}
> >>
> >> https://en.wikipedia.org/wiki/Zero-inflated_model#Zero-inflated_Poisson
> >>
> >> On 22 Mar 2016, at 14:25, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:
> >>
> >> Could you clarify what are the parameters and why it’s formulated that way?
> >>
> >> -Matti
> >>
> >> On 22 Mar 2016, at 14:17, Thierry Onkelinx <thierry.onkelinx at inbo.be>
> >> wrote:
> >>
> >> Dear Matti,
> >>
> >> What about this?
> >>
> >> dzeroinflpois <- function(x, lambda, zero){
> >> ifelse(x == 0, zero, 0) + dpois(x, lambda) / (1 - zero)
> >> }
> >> plot(x, dzeroinflpois(x, lambda = 10, zero = 0.2), type = "l")
> >>
> >>
> >>
> >> ir. Thierry Onkelinx
> >> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> >> Forest
> >> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> >> Kliniekstraat 25
> >> 1070 Anderlecht
> >> Belgium
> >>
> >> To call in the statistician after the experiment is done may be no more
> >> than asking him to perform a post-mortem examination: he may be able to say
> >> what the experiment died of. ~ Sir Ronald Aylmer Fisher
> >> The plural of anecdote is not data. ~ Roger Brinner
> >> The combination of some data and an aching desire for an answer does not
> >> ensure that a reasonable answer can be extracted from a given body of data.
> >> ~ John Tukey
> >>
> >> 2016-03-22 13:04 GMT+01:00 Matti Viljamaa <mviljamaa at kapsi.fi>:
> >>
> >>> I’m doing some optimisation that I first did with normal Poisson (only
> >>> parameter theta was estimated), but now I’m doing the same with a
> >>> zero-inflated Poisson model which
> >>> gives me two estimated parameters theta and p (p is also pi in some
> >>> notation).
> >>>
> >>> My question is, is there something equivalent to dpois that would use
> >>> both of the parameters (or is the p parameter possibly unnecessary)?
> >>>
> >>> I’m calculating the “fit” of the Poisson model
> >>>
> >>> i.e. like
> >>>
> >>> x = c(0,1,2,3,4,5,6)
> >>> y = c(3062,587,284,103,33,4,2)
> >>> fit1 <- sum(y)*dpois(x, est_theta)
> >>>
> >>> and then comparing fit1 to the real observations.
> >>> [[alternative HTML version deleted]]
> >>>
> >>> ______________________________________________
> >>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide
> >>> http://www.R-project.org/posting-guide.html
> >>> <http://www.r-project.org/posting-guide.html>
> >>> and provide commented, minimal, self-contained, reproducible code.
> >>
> >>
> >>
> >>
>
> > [[alternative HTML version deleted]]
>
> > ______________________________________________
> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
More information about the R-help
mailing list