[Rd] rpois gives a large number repeatedly (PR#439)

[kjetikj@astro.uio.no]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._