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.onkel...@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 <mvilja...@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 <mvilja...@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.onkel...@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 <mvilja...@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@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@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@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@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.
