pythonprobabilityconvolutionoff-by-one

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


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?


Solution

  • 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)