[R] Issues with convolve
Gabor Grothendieck
ggrothendieck at gmail.com
Thu Jul 21 01:29:24 CEST 2005
Check out sum.exact in the caTools package.
On 7/20/05, Vincent Goulet <vincent.goulet at 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 at 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
>
More information about the R-help
mailing list