kjetikj@astro.uio.no
2000-Feb-15 11:18 UTC
[Rd] rpois gives a large number repeatedly (PR#439)
Full_Name: Kjetil Kjernsmo
Version: 0.65.1
OS: Digital UNIX 4.0
Submission from: (NULL) (129.240.28.172)
I'm experiencing problems with rpois. Occasionally, it draws a very high
number.
Yeah, I know, this is statistics, things like that does happen, but this really
strange because a poisson distribution with a parameter of 3 shouldn't see
the
number 1932 very often, but the same, incredibly high number is drawn
repeatedly.
These are the two relevant functions, the latter somewhat expanded for debugging
purposes.
ncloudsbin <- function(binno, ntotalclouds, numberofbins = 100, linewidth
numberofbins / 6)
return(ntotalclouds * dnorm(binno, sd = linewidth))
photonsincidenttodetectorbin <- function(binno, ntotalclouds,
avgphotonfromcloud, ampmode = 1, ampmean = 2, numberofbins = 100, linewidth
numberofbins / 6)
{
sc <- ampmean - ampmode
if(sc <= 0) stop("Mean Amplification must be greater than Mode
Amplification")
ncb <- ncloudsbin(binno, ntotalclouds, numberofbins, linewidth)
rg <- rgamma(ncb, ampmean / sc, sc)
rp <- rpois(ncb, avgphotonfromcloud)
dr <- rg * rp
if(sum(dr) > 25000) print(cbind(rg, rp, dr, avgphotonfromcloud, ncb))
return(sum(dr))
}
The if-line is included to trigger unusual events, and I just had output like
this:
[683,] 2.09611150 1932 4.049687e+03 3 1682.063
[684,] 0.40671731 2 8.134346e-01 3 1682.063
[685,] 2.92849186 3 8.785476e+00 3 1682.063
[686,] 1.18576903 2 2.371538e+00 3 1682.063
[687,] 1.57518417 2 3.150368e+00 3 1682.063
[688,] 3.72793434 2 7.455869e+00 3 1682.063
[689,] 0.49290765 3 1.478723e+00 3 1682.063
[690,] 2.25682070 2 4.513641e+00 3 1682.063
[691,] 4.37234292 0 0.000000e+00 3 1682.063
[692,] 2.86349873 6 1.718099e+01 3 1682.063
[693,] 2.04314101 5 1.021571e+01 3 1682.063
[694,] 1.74437433 1932 3.370131e+03 3 1682.063
[695,] 0.71229856 1 7.122986e-01 3 1682.063
[696,] 4.08736383 1 4.087364e+00 3 1682.063
[697,] 3.28338545 2 6.566771e+00 3 1682.063
[698,] 2.74590789 7 1.922136e+01 3 1682.063
[699,] 0.81551877 1932 1.575582e+03 3 1682.063
Actually, it seems like this large number may dependent of the session one is
in.
I get 1932 on an XP1000 with an 500MHz alphaev6, but on my own computer, which
is an PWS500au with
an 500MHz alphaev56, I have seen the number 11344 repeatedly. It could of course
be the
seed or something, I guess.
Kjetil
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Peter Dalgaard BSA
2000-Feb-15 11:34 UTC
[Rd] rpois gives a large number repeatedly (PR#439)
kjetikj@astro.uio.no writes:> Full_Name: Kjetil Kjernsmo > Version: 0.65.1!!!!!!> OS: Digital UNIX 4.0 > Submission from: (NULL) (129.240.28.172) >> ncloudsbin <- function(binno, ntotalclouds, numberofbins = 100, linewidth > numberofbins / 6) > return(ntotalclouds * dnorm(binno, sd = linewidth)) > > photonsincidenttodetectorbin <- function(binno, ntotalclouds, > avgphotonfromcloud, ampmode = 1, ampmean = 2, numberofbins = 100, linewidth > numberofbins / 6) > { > sc <- ampmean - ampmode > if(sc <= 0) stop("Mean Amplification must be greater than Mode Amplification") > > ncb <- ncloudsbin(binno, ntotalclouds, numberofbins, linewidth) > rg <- rgamma(ncb, ampmean / sc, sc) > rp <- rpois(ncb, avgphotonfromcloud) > dr <- rg * rp > if(sum(dr) > 25000) print(cbind(rg, rp, dr, avgphotonfromcloud, ncb)) > return(sum(dr)) > } >I can't seem to generate anything unusual with the latest R snapshot. Perhaps it would help if you told us what parameters to call photonsincidenttodetectorbin with? -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._