How to deal with off-by-one issues in convolution (Python)?

Question

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?

Answer 1

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)