[R] Issues with convolve

[Vincent Goulet]

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