I'm trying to write a function to add two random varibles X1
and X2
. In my case, they are both uniform random variables from 0
to a1
and 0
to a2
. To compute the random variable Y = X1 + X2
, I need to take a convolution of the probability distributions of X1
and X2
.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
def convolution(f, g, x_range):
delta = (x_range[-1] - x_range[0])/len(x_range)
result = np.convolve(f(x_range), g(x_range), mode = 'full')*delta
return result
# Define uniform distribution for some a > 0. This part can be adapted to arbitrary distributions
def uniform_dist(x, a):
return np.where((x >= 0) & (x <= a), 1/a, 0)
# Set the range of x values, y values and constants
delta = 0.1
x_lim_low = -5
x_lim_upp = 5
a1 = 1
a2 = 1
x_range = np.arange(x_lim_low,x_lim_upp+delta,delta)
y_range = np.arange(2*x_lim_low,2*x_lim_upp+delta,delta)
# Perform convolution
convolution_pdf = convolution(lambda x: uniform_dist(x, a1), lambda x: uniform_dist(x, a2), x_range)
# Find mean of convolution
convolution_mean = np.sum(convolution_pdf*y_range)*delta
I've tried various combinations but have small errors in the mean. I think this is because the convolution is an array of dimension 2*len(x_range) - 1
and it's unclear how to deal with this off by one error.
What is the correct way to convolve to variables such that I can compute the mean of the convolution correctly?
in convolution
you calculate the delta
incorrect.
to get nicer sample points don't use np.arange
but np.linspace
now convolution_mean = 1.21
import numpy as np
def convolution(f, g, x_range):
delta = x_range[1]-x_range[0]
return np.convolve(f(x_range), g(x_range), mode = 'full') * delta
# Define uniform distribution for some a > 0. This part can be adapted to arbitrary distributions
def uniform_dist(x, a):
return np.where((x >= 0) & (x <= a), 1/a, 0)
# Set the range of x values, y values and constants
delta = 0.1
one_over_delta = 10
x_lim_low = -5
x_lim_upp = 5
a1 = 1
a2 = 1
x_range = np.linspace(x_lim_low*one_over_delta, x_lim_upp*one_over_delta, (x_lim_upp - x_lim_low)*one_over_delta + 1) / one_over_delta
y_range = np.linspace(2*x_lim_low*one_over_delta, 2*x_lim_upp*one_over_delta, 2*(x_lim_upp - x_lim_low)*one_over_delta + 1) / one_over_delta
# Perform convolution
convolution_pdf = convolution(lambda x: uniform_dist(x, a1), lambda x: uniform_dist(x, a2), x_range)
# Find mean of convolution
convolution_mean = np.sum(convolution_pdf * y_range) * delta
print(convolution_mean)