Thierry Onkelinx
2016-Mar-22 12:58 UTC
[R] Is there dpois equivalent for zero-inflated Poisson?
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) 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]]
Martin Maechler
2016-Mar-23 15:27 UTC
[R] Is there dpois equivalent for zero-inflated Poisson?
>>>>> 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.
Michael Friendly
2016-Mar-24 13:06 UTC
[R] Is there dpois equivalent for zero-inflated Poisson?
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. >