From: Gabor Grothendieck <ggrothendieck_at_gmail.com>

Date: Thu 21 Jul 2005 - 09:29:24 EST

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Jul 21 09:39:52 2005

Date: Thu 21 Jul 2005 - 09:29:24 EST

Check out sum.exact in the caTools package.

On 7/20/05, Vincent Goulet <vincent.goulet@act.ulaval.ca> wrote:

> We obtained some disturbing results from convolve() (inaccuracies and negative

*> probabilities). We'll try to make the context clear in as few lines as
**> possible...
**>
**> Our function panjer() (code below) basically computes recursively the
**> probability mass function of a compound Poisson distribution. When the
**> Poisson parameter lambda is very large, the starting value of the recursive
**> scheme --- the mass at 0 --- is 0 and the recursion fails. One way to solve
**> this problem is to divide lambda by 2^n, apply the panjer() function and then
**> convolve the result with itself n times.
**>
**> We applied the panjer() function with a lambda such that the mass at 0 is just
**> larger than .Machine$double.xmin. We thus know that once this pmf is
**> convoluted with itself, the first probabilities will be 0 (for the computer).
**>
**> Here are the two issues we have with convolve():
**>
**> 1. The probabilities we know should be 0 are rather in the vicinity of 1E-19,
**> as if convolve() could not "go lower". Using a hand made convolution function
**> (not given here), we obtained the correct values. When probabilities get
**> around 1E-12, results from convolve() and our home made function are
**> essentially identical.
**>
**> 2. We obtained negative probabilities. More accurately, the same example
**> returns negative probabilities under Windows, but not under Linux. We also
**> obtained negative probabilities for another example under Linux, though.
**>
**> We understand that convolve() computes the convolutions using fft(), but we
**> are not familiar enough with the latter to assess if the above issues are
**> some sort of bugs or normal behavior. In the latter case, is there is any
**> workaround?
**>
**> Any help/comments/ideas appreciated.
**>
**> === EXAMPLES ===
**> panjer <- function(fx, lambda)
**> {
**> fx0 <- fx[1]
**> fx <- fx[-1]
**> r <- length(fx)
**> cumul <- fs <- exp(-lambda * (1 - fx0))
**> repeat
**> {
**> x <- length(fs)
**> m <- min(x, r)
**> fs <- c(fs, sum((lambda * 1:m / x) * head(fx, m) * rev(tail(fs, m))))
**> if ((1-1E-8) < (cumul <- cumul + tail(fs, 1)))
**> break
**> }
**> fs
**> }
**>
**> ### Linux example ###
**> > R.version
**> _
**> platform i386-pc-linux-gnu
**> arch i386
**> os linux-gnu
**> system i386, linux-gnu
**> status
**> major 2
**> minor 1.0
**> year 2005
**> month 04
**> day 18
**> language R
**> > fs <- panjer(rep(0.2, 5), 600) # compute pmf
**> > fs[1:10] # small values
**> [1] 3.456597e-209 4.147916e-207 2.530229e-205 1.045690e-203
**> [5] 3.292657e-202 8.422924e-201 1.822768e-199 3.431181e-198
**> [9] 5.733481e-197 8.637025e-196
**> > fsc <- convolve(fs, rev(fs), type="o") # convolution
**> > sum(fsc < 0) # no negatives under Linux
**> [1] 0
**> > fsc[1:10] # these should be 0
**> [1] 4.819870e-19 4.135471e-19 4.511455e-19 3.994485e-19 3.820177e-19
**> [6] 4.738272e-19 3.885447e-19 3.167399e-19 3.695636e-19 3.701792e-19
**> > sum(fsc) # no impact on the sum
**> [1] 1
**>
**> ### Windows example ###
**> > R.version
**> _
**> platform i386-pc-mingw32
**> arch i386
**> os mingw32
**> system i386, mingw32
**> status
**> major 2
**> minor 1.1
**> year 2005
**> month 06
**> day 20
**> language R
**> > fs <- panjer(rep(0.2, 5), 600) # compute pmf
**> > fs[1:10] # small values
**> [1] 3.456597e-209 4.147916e-207 2.530229e-205 1.045690e-203
**> [5] 3.292657e-202 8.422924e-201 1.822768e-199 3.431181e-198
**> [9] 5.733481e-197 8.637025e-196
**> > fsc <- convolve(fs, rev(fs), type="o") # convolution
**> > sum(fsc < 0) # negatives under Windows
**> [1] 39
**> > fsc[1:10] # these should be 0
**> [1] 4.049552e-19 3.345584e-19 4.126913e-19 3.351211e-19 2.758626e-19
**> [6] 3.530111e-19 2.735041e-19 2.376711e-19 2.591287e-19 3.196405e-19
**> > sum(fsc) # no impact on the sum
**> [1] 1
**>
**> --
**> Vincent Goulet, Associate Professor
**> École d'actuariat
**> Université Laval, Québec
**> Vincent.Goulet_at_act.ulaval.ca http://vgoulet.act.ulaval.ca
**>
**> ______________________________________________
**> R-help@stat.math.ethz.ch mailing list
**> https://stat.ethz.ch/mailman/listinfo/r-help
**> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
*

>

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Jul 21 09:39:52 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:33:53 EST
*