# [S] SUMMARY: Numerical convolution

joyet pierre (pierre.joyet@basler.ch)
Fri, 31 Jul 1998 12:46:43 +0200

Dies ist eine mehrteilige Nachricht im MIME-Format.
--------------735FDDDB08298620660E04B5
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Yesterday I asked the following question:

>
> Using S-Plus 4.5, what is the most efficient numerical method for
> computing (in 1000 points, say) the density or probability function of
> X
> + Y, where X and Y are random variables whose joint density or
> probability function is given?
>
> Is there a special method when X and Y are independent?
>
> Please tell me if You know of a good reference about such questions.
>

Many thanks to Andreas Krause, Bob Henery, Andrey Feuerverger and Peter
Perkins for their helpful suggestions which I reproduce below.
The fast Fourier transform method (cf. Bob's answer) is what I was
looking for.

Regards

Pierre Joyet
Baloise Insurance Company
Switzerland

-------------------------------------------------------------------
>From Andreas Krause:

If you want an exact result, you need to integrate the other variable
out. f(x)
= integral f(x,y) dy. You can do that to some extent in Splus. For more
complicate functions, I'd rather use a symbolic computation package like
mathematica.
If an approximation is enough, just draw a big sample of (x,y) jointly
and plot
a density estimate of the x component (like a hist or a nonparametric
density
estimate like ksmooth).

A good reference for you might be Wolfgang Haerdle, Smoothing Techniques
with
Implementation in S, Springer, 1991, although it is somewhat outdated
with
respect to S/Splus.

-------------------------------------------------------------------
>From Bob Henery:

There are various interpretations of "easy".

Here are two "solutions", illustrated for the convolution of two
discrete
binomials X1 ~ binom(10,0.3), X2 ~ binom(15,0.3) with X1, X2
independent.

d1 <- dbinom(c(0:10),10,0.3) # density X1
d2 <- dbinom(c(0:15),15,0.3) # density X2
d3 <- dbinom(c(0:25),25,0.3) # the correct answer for density
X1+X2

1. small number of Splus commands but also works for non-independent
RV's

dd <- kronecker(d1,d2)
ss <- kronecker(c(0:10),c(0:15),fun="+")
ddd3 <- tapply(dd,ss,sum)

2. computationally fast for large vectors:

# pad out the densities with zeros to give
# length 128 (more than 2*length(d1)+length(d2))
ed1 <- c(d1,rep(0,117))
ed2 <- c(d2,rep(0,112))
ed3 <- c(d3,rep(0,102))

Td1 <- fft(ed1)
Td2 <- fft(ed2)

FFTd3 <- fft(Td1*Td2,inv=T)/128

fffd3 <- Re(FFTd3)[1:26]

-------------------------------------------------------------------
>From Andrey Feuerverger:

a good reference in general for random variable generation
is the book by luc devroye called something like `nonuniform
random variable generation...'
usually things will be best if you can work out the distribution
of X+Y analytically.

-------------------------------------------------------------------
>From Peter Perkins:

not exactly sure what you're trying to do, but the following works as an
approximation to a convolution if x and y are independent:

# some dummy data
x <- seq(-5, 5, length = 101)
px <- dnorm(x)
y <- seq(0, 10, length = 101)
py <- dnorm(y, mean = 5)
#
# the convolution itself
pz <- tapply(outer(px, py, "*"), round(outer(x, y, "+"), 6),
sum)
z <- as.numeric(names(pz))

you're going to run into memory problems if your bin size gets too
small.
adjust the digits arg to round() as needed.

if x and y are not independent, it may be easiest to just use
monte-carlo.
--------------735FDDDB08298620660E04B5
Content-Type: message/rfc822
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

Message-ID: <35C02F75.1A4C58F5@basler.ch>
Date: Thu, 30 Jul 1998 10:31:49 +0200
From: Pierre Joyet <pierre.joyet@basler.ch>
X-Mailer: Mozilla 4.01 [de] (Win95; I)
MIME-Version: 1.0
To: s-news@wubios.wustl.edu
Subject: Numerical convolution
X-Priority: 3 (Normal)
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Using S-Plus 4.5, what is the most efficient numerical method for
computing (in 1000 points, say) the density or probability function of X
+ Y, where X and Y are random variables whose joint density or
probability function is given?

Is there a special method when X and Y are independent?

Please tell me if You know of a good reference about such questions.