How to get a convolution of 3 or more continuous random variables to obtain the average of the rando-CodePudding

Say I have three random variables. I would like to do a convolution to obtain the average. How do I do this in Python and or R?

Edit 1

Also. It seems the default behavior is to have the convolution size larger than any of the inputs. I will assume that all of the inputs are the same size. Is it possible to have the resulting convolution the same size as the vectors which are being used as inputs to the convolution?

For example, if x1 is n=100 then I would like the resulting convolution to be n=100

Edit 2 - Added Example

I theory the convolution should be close to what I can calculate analytically.

import numpy as np
rng = np.random.default_rng(42)
n, u1, u2, u3, sd = 100, 10, 20, 6, 5
u_avg = np.mean([u1,u2,u3])
a = rng.normal(u1, sd, size=n)
b = rng.normal(u2, sd, size=n)
c = rng.normal(u3, sd, size=n)
z = rng.normal(u_avg, sd/np.sqrt(2), size=n)
convolution = rng.choice(reduce(np.convolve, [a, b, c]), size=n)

print("true distribution")
print(np.round(np.quantile(z, [0.01, 0.25, 0.5, 0.75, 0.99]), 2))
print("convolution")
print(np.round(np.quantile(convolution, [0.01, 0.25, 0.5, 0.75, 0.99]),2))

If the convolution is working then the convolution should be close to the true distribution.

true distribution
[ 3.9   9.84 12.83 14.89 18.45]
convolution
[5.73630000e 03 5.47855750e 05 2.15576037e 06 6.67763665e 06
 8.43843281e 06]

It looks like the convolution is not even close.

CodePudding user response：

In Python:

import numpy as np
from functools import reduce

x1 = np.random.normal(0, 1, 100)
x2 = np.random.normal(0, 1, 100)
x3 = np.random.normal(0, 1, 100)

result = reduce(np.convolve, [x1, x2, x3])

In R:


convolve3 <- function(x1, x2, x3) {
  result <- convolve(x1, x2)
  result <- convolve(result, x3)
  return(result)
}

x1 <- rnorm(100)
x2 <- rnorm(100)
x3 <- rnorm(100)

result <- Reduce(convolve3, list(x1, x2, x3))

CodePudding user response：

You can try conv from package pracma

> library(pracma)

> x <- c(1, 2)

> y <- c(2, 3, 4)

> z <- c(-1, 0, 1, 5)

> Reduce(conv, list(x, y, z))
[1] -2 -7 -8  9 45 58 40