????? This issue arose for me in simulations to estimate confidence, prediction, and tolerance intervals from glm(., family=poisson) fits embedded in a BMA::bic.glm fit using a simulate.bic.glm function I added to the development version of Ecfun, available at "https://github.com/sbgraves237/Ecfun".? This is part of a vignette I'm developing, available at "https://github.com/sbgraves237/Ecfun/blob/master/vignettes/time2nextNuclearWeaponState.Rmd". This includes a simulated mean of a mixture of Poissons that exceeds 2e22.? It doesn't seem unreasonable to me to have rpois output a numerics rather than integers when a number simulated exceeds .Machine$integer.max.? And it does seem to make less sense in such cases to return NAs. ?????? Alternatively, might it make sense to add another argument to rpois to give the user the choice?? E.g., an argument "bigOutput" with (I hope) default = "numeric" and "NA" as a second option.? Or NA is the default, so no code that relied that feature of the current code would be broken by the change.? If someone wanted to use arbitrary precision arithmetic, they could write their own version of this function with "arbitraryPrecision" as an optional value for the "bigOutput" argument. ????? Comments? ????? Thanks, ????? Spencer Graves On 2020-01-19 10:28, Avraham Adler wrote:> Technically, lambda can always be numeric. It is the observations > which must be integral. > > Would hitting everything larger than maxint or maxlonglong with floor > or round fundamentally change the distribution? Well, yes, but enough > that it would matter over process risk? > > Avi > > On Sun, Jan 19, 2020 at 11:20 AM Benjamin Tyner <btyner at gmail.com > <mailto:btyner at gmail.com>> wrote: > > So imagine rpois is changed, such that the storage mode of its return > value is sometimes integer and sometimes numeric. Then imagine the > case > where lambda is itself a realization of a random variable. Do we > really > want the storage mode to inherit that randomness? > > > On 1/19/20 10:47 AM, Avraham Adler wrote: > > Maybe there should be code for 64 bit R to use long long or the > like? > > > > On Sun, Jan 19, 2020 at 10:45 AM Spencer Graves > > <spencer.graves at prodsyse.com > <mailto:spencer.graves at prodsyse.com> > <mailto:spencer.graves at prodsyse.com > <mailto:spencer.graves at prodsyse.com>>> wrote: > > > > > > > >? ? ?On 2020-01-19 09:34, Benjamin Tyner wrote: > >? ? ?>> > > > ?------------------------------------------------------------------------ > >? ? ?>> Hello, All: > >? ? ?>> > >? ? ?>> > >? ? ?>> ? ????? Consider: > >? ? ?>> > >? ? ?>> > >? ? ?>> Browse[2]> set.seed(1) > >? ? ?>> Browse[2]> rpois(9, 1e10) > >? ? ?>> NAs produced[1] NA NA NA NA NA NA NA NA NA > >? ? ?>> > >? ? ?>> > >? ? ?>> ? ????? Should this happen? > >? ? ?>> > >? ? ?>> > >? ? ?>> ? ????? I think that for, say, lambda>1e6, rpois should > return > >? ? ?rnorm(., > >? ? ?>> lambda, sqrt(lambda)). > >? ? ?> But need to implement carefully; rpois should always return a > >? ? ?> non-negative integer, whereas rnorm always returns numeric... > >? ? ?> > > > >? ? ??????? Thanks for the reply. > > > > > >? ? ??????? However, I think it's not acceptable to get an NA from a > >? ? ?number > >? ? ?that cannot be expressed as an integer.? Whenever a randomly > >? ? ?generated > >? ? ?number would exceed .Machine$integer.max, the choice is between > >? ? ?returning NA or a non-integer numeric.? Consider: > > > > > >? ? ??> 2*.Machine$integer.max > >? ? ?[1] 4294967294 > >? ? ??> as.integer(2*.Machine$integer.max) > >? ? ?[1] NA > >? ? ?Warning message: > >? ? ?NAs introduced by coercion to integer range > > > > > >? ? ??????? I'd rather have the non-integer numeric. > > > > > >? ? ??????? Spencer > > > >? ? ?______________________________________________ > > R-devel at r-project.org <mailto:R-devel at r-project.org> > <mailto:R-devel at r-project.org <mailto:R-devel at r-project.org>> > mailing list > > https://stat.ethz.ch/mailman/listinfo/r-devel > > > > -- > > Sent from Gmail Mobile > > -- > Sent from Gmail Mobile[[alternative HTML version deleted]]
Crazy thought, but being that a sum of Poissons is Poisson in the sum, can you break your ?big? simulation into the sum of a few smaller ones? Or is the order of magnitude difference just too great? On Sun, Jan 19, 2020 at 1:58 PM Spencer Graves <spencer.graves at prodsyse.com> wrote:> This issue arose for me in simulations to estimate confidence, > prediction, and tolerance intervals from glm(., family=poisson) fits > embedded in a BMA::bic.glm fit using a simulate.bic.glm function I added to > the development version of Ecfun, available at > "https://github.com/sbgraves237/Ecfun" > <https://github.com/sbgraves237/Ecfun>. This is part of a vignette I'm > developing, available at > "https://github.com/sbgraves237/Ecfun/blob/master/vignettes/time2nextNuclearWeaponState.Rmd" > <https://github.com/sbgraves237/Ecfun/blob/master/vignettes/time2nextNuclearWeaponState.Rmd>. > This includes a simulated mean of a mixture of Poissons that exceeds 2e22. > It doesn't seem unreasonable to me to have rpois output a numerics rather > than integers when a number simulated exceeds .Machine$integer.max. And it > does seem to make less sense in such cases to return NAs. > > > Alternatively, might it make sense to add another argument to rpois > to give the user the choice? E.g., an argument "bigOutput" with (I hope) > default = "numeric" and "NA" as a second option. Or NA is the default, so > no code that relied that feature of the current code would be broken by the > change. If someone wanted to use arbitrary precision arithmetic, they > could write their own version of this function with "arbitraryPrecision" as > an optional value for the "bigOutput" argument. > > > Comments? > Thanks, > Spencer Graves > > > > On 2020-01-19 10:28, Avraham Adler wrote: > > Technically, lambda can always be numeric. It is the observations which > must be integral. > > Would hitting everything larger than maxint or maxlonglong with floor or > round fundamentally change the distribution? Well, yes, but enough that it > would matter over process risk? > > Avi > > On Sun, Jan 19, 2020 at 11:20 AM Benjamin Tyner <btyner at gmail.com> wrote: > >> So imagine rpois is changed, such that the storage mode of its return >> value is sometimes integer and sometimes numeric. Then imagine the case >> where lambda is itself a realization of a random variable. Do we really >> want the storage mode to inherit that randomness? >> >> >> On 1/19/20 10:47 AM, Avraham Adler wrote: >> > Maybe there should be code for 64 bit R to use long long or the like? >> > >> > On Sun, Jan 19, 2020 at 10:45 AM Spencer Graves >> > <spencer.graves at prodsyse.com <mailto:spencer.graves at prodsyse.com>> >> wrote: >> > >> > >> > >> > On 2020-01-19 09:34, Benjamin Tyner wrote: >> > >> >> > >> ------------------------------------------------------------------------ >> > >> Hello, All: >> > >> >> > >> >> > >> Consider: >> > >> >> > >> >> > >> Browse[2]> set.seed(1) >> > >> Browse[2]> rpois(9, 1e10) >> > >> NAs produced[1] NA NA NA NA NA NA NA NA NA >> > >> >> > >> >> > >> Should this happen? >> > >> >> > >> >> > >> I think that for, say, lambda>1e6, rpois should return >> > rnorm(., >> > >> lambda, sqrt(lambda)). >> > > But need to implement carefully; rpois should always return a >> > > non-negative integer, whereas rnorm always returns numeric... >> > > >> > >> > Thanks for the reply. >> > >> > >> > However, I think it's not acceptable to get an NA from a >> > number >> > that cannot be expressed as an integer. Whenever a randomly >> > generated >> > number would exceed .Machine$integer.max, the choice is between >> > returning NA or a non-integer numeric. Consider: >> > >> > >> > > 2*.Machine$integer.max >> > [1] 4294967294 >> > > as.integer(2*.Machine$integer.max) >> > [1] NA >> > Warning message: >> > NAs introduced by coercion to integer range >> > >> > >> > I'd rather have the non-integer numeric. >> > >> > >> > Spencer >> > >> > ______________________________________________ >> > R-devel at r-project.org <mailto:R-devel at r-project.org> mailing list >> > https://stat.ethz.ch/mailman/listinfo/r-devel >> > >> > -- >> > Sent from Gmail Mobile >> > -- > Sent from Gmail Mobile > > > --Sent from Gmail Mobile [[alternative HTML version deleted]]
On 2020-01-19 13:01, Avraham Adler wrote:> Crazy thought, but being that a sum of Poissons is Poisson in the sum, > can you break your ?big? simulation into the sum of a few smaller > ones? Or is the order of magnitude difference just too great?????? I don't perceive that as feasible.? Once I found what was generating NAs, it was easy to code a function to return pseudo-random numbers using the standard normal approximation to the Poisson for those extreme cases.? [For a Poisson with mean = 1e6, for example, the skewness (third standardized moment) is 0.001.? At least for my purposes, that should be adequate.][1] ????? What are the negative consequences of having rpois return numerics that are always nonnegative? ????? Spencer [1]? In the code I reported before, I just changed the threshold of 1e6 to 0.5*.Machine$integer.max.? On my Mac, .Machine$integer.max = 2147483647 = 2^31 > 1e9.? That still means that a Poisson distributed pseudo-random number just under that would have to be over 23000 standard deviations above the mean to exceed .Machine$integer.max.> > On Sun, Jan 19, 2020 at 1:58 PM Spencer Graves > <spencer.graves at prodsyse.com <mailto:spencer.graves at prodsyse.com>> wrote: > > ????? This issue arose for me in simulations to estimate > confidence, prediction, and tolerance intervals from glm(., > family=poisson) fits embedded in a BMA::bic.glm fit using a > simulate.bic.glm function I added to the development version of > Ecfun, available at "https://github.com/sbgraves237/Ecfun" > <https://github.com/sbgraves237/Ecfun>. This is part of a vignette > I'm developing, available at > "https://github.com/sbgraves237/Ecfun/blob/master/vignettes/time2nextNuclearWeaponState.Rmd" > <https://github.com/sbgraves237/Ecfun/blob/master/vignettes/time2nextNuclearWeaponState.Rmd>. > This includes a simulated mean of a mixture of Poissons that > exceeds 2e22.? It doesn't seem unreasonable to me to have rpois > output a numerics rather than integers when a number simulated > exceeds .Machine$integer.max.? And it does seem to make less sense > in such cases to return NAs. > > > ?????? Alternatively, might it make sense to add another argument > to rpois to give the user the choice?? E.g., an argument > "bigOutput" with (I hope) default = "numeric" and "NA" as a second > option.? Or NA is the default, so no code that relied that feature > of the current code would be broken by the change.? If someone > wanted to use arbitrary precision arithmetic, they could write > their own version of this function with "arbitraryPrecision" as an > optional value for the "bigOutput" argument. > > > ????? Comments? > ????? Thanks, > ????? Spencer Graves > > > > On 2020-01-19 10:28, Avraham Adler wrote: >> Technically, lambda can always be numeric. It is the observations >> which must be integral. >> >> Would hitting everything larger than maxint or maxlonglong with >> floor or round fundamentally change the distribution? Well, yes, >> but enough that it would matter over process risk? >> >> Avi >> >> On Sun, Jan 19, 2020 at 11:20 AM Benjamin Tyner <btyner at gmail.com >> <mailto:btyner at gmail.com>> wrote: >> >> So imagine rpois is changed, such that the storage mode of >> its return >> value is sometimes integer and sometimes numeric. Then >> imagine the case >> where lambda is itself a realization of a random variable. Do >> we really >> want the storage mode to inherit that randomness? >> >> >> On 1/19/20 10:47 AM, Avraham Adler wrote: >> > Maybe there should be code for 64 bit R to use long long or >> the like? >> > >> > On Sun, Jan 19, 2020 at 10:45 AM Spencer Graves >> > <spencer.graves at prodsyse.com >> <mailto:spencer.graves at prodsyse.com> >> <mailto:spencer.graves at prodsyse.com >> <mailto:spencer.graves at prodsyse.com>>> wrote: >> > >> > >> > >> >? ? ?On 2020-01-19 09:34, Benjamin Tyner wrote: >> >? ? ?>> >> > >> ?------------------------------------------------------------------------ >> >? ? ?>> Hello, All: >> >? ? ?>> >> >? ? ?>> >> >? ? ?>> ? ????? Consider: >> >? ? ?>> >> >? ? ?>> >> >? ? ?>> Browse[2]> set.seed(1) >> >? ? ?>> Browse[2]> rpois(9, 1e10) >> >? ? ?>> NAs produced[1] NA NA NA NA NA NA NA NA NA >> >? ? ?>> >> >? ? ?>> >> >? ? ?>> ? ????? Should this happen? >> >? ? ?>> >> >? ? ?>> >> >? ? ?>> ? ????? I think that for, say, lambda>1e6, rpois >> should return >> >? ? ?rnorm(., >> >? ? ?>> lambda, sqrt(lambda)). >> >? ? ?> But need to implement carefully; rpois should always >> return a >> >? ? ?> non-negative integer, whereas rnorm always returns >> numeric... >> >? ? ?> >> > >> >? ? ??????? Thanks for the reply. >> > >> > >> >? ? ??????? However, I think it's not acceptable to get an >> NA from a >> >? ? ?number >> >? ? ?that cannot be expressed as an integer. Whenever a randomly >> >? ? ?generated >> >? ? ?number would exceed .Machine$integer.max, the choice is >> between >> >? ? ?returning NA or a non-integer numeric. Consider: >> > >> > >> >? ? ??> 2*.Machine$integer.max >> >? ? ?[1] 4294967294 >> >? ? ??> as.integer(2*.Machine$integer.max) >> >? ? ?[1] NA >> >? ? ?Warning message: >> >? ? ?NAs introduced by coercion to integer range >> > >> > >> >? ? ??????? I'd rather have the non-integer numeric. >> > >> > >> >? ? ??????? Spencer >> > >> > ?______________________________________________ >> > R-devel at r-project.org <mailto:R-devel at r-project.org> >> <mailto:R-devel at r-project.org <mailto:R-devel at r-project.org>> >> mailing list >> > https://stat.ethz.ch/mailman/listinfo/r-devel >> > >> > -- >> > Sent from Gmail Mobile >> >> -- >> Sent from Gmail Mobile > > -- > Sent from Gmail Mobile[[alternative HTML version deleted]]