Convolution of a non-central chi-square with a central chi-square distribution

Fred Ahrens

5 July 2018

Abstract

I need an expression for the convolution of a non-central chi-square distribution with a central chi-square, where the two have different scale factors. The expression will be used to implement a Python function of the result's density function for use in maximum likelihood estimation of GPS movement data.

$\nu_{X},\nu_{Y}$ (nuX, nuY) Degrees of freedom of the two distributions.
$\sigma_{X}^{2},\sigma_{Y}^{2}$ (sigmaX2, sigmaY2) Scale factors of the two distributions.
$\delta_{X}$ (deltaX) Non-centrality parameter of the first distribution.

# These settings motivated by the GPS motion detection problem.
nuX = 1.
nuY = 1.
sigmaY2 = 10.*10 # m^2, variance of measuring device
sigmaX2 = 2.7*2.7 + sigmaY2 # combined variance of distance traveled and measurement error.
deltaX = 27.*27/sigmaX2

import numpy as np
from scipy.integrate import quad
from scipy.special import iv, gamma

Convolution Integral

For any z>0, the density of Z is

$$f_{Z}\left(z\right)=\frac{1}{2\sigma_{X}^{2}\sigma_{Y}^{2}2^{\frac{\nu_{Y}}{2}}\Gamma\left(\frac{\nu_{Y}}{2}\right)}\int_{0}^{z}\left(\frac{z-y}{\delta_{X}\sigma_{X}^{2}}\right)^{\frac{\nu_{x}}{4}-\frac{1}{2}}\left(\frac{y}{\sigma_{Y}^{2}}\right)^{\frac{\nu_{Y}}{2}-1}\exp\left(-\frac{1}{2}\left(\frac{z-y}{\sigma_{X}^{2}}+\delta_{X}+\frac{y}{\sigma_{Y}^{2}}\right)\right)I_{\frac{\nu_{x}}{2}-1}\left(\sqrt{\frac{\delta_{X}\left(z-y\right)}{\sigma_{X}^{2}}}\right)dy$$

The next cell is a function to calculate the integrand of this integral.

y is the variable of integration.

six parameters are passed in addition to y, which are needed to evaluate the integrand:
$\left(\nu_{x},\sigma_{X}^{2},\delta_{X},\nu_{Y},\sigma_{Y}^{2},z\right)$.

# this function definition is the integrand
def f_non_central_joint(y, nuX, sigmaX2, deltaX, nuY, sigmaY2, z):
    xnorm = (z - y)/sigmaX2
    ynorm = y/sigmaY2
    f = (xnorm/deltaX)**(nuX/4 - 0.5)*ynorm**(nuY/2 - 1)* \
            np.exp(-0.5*(xnorm + deltaX + ynorm))* \
            iv(nuX/2 - 1, np.sqrt(deltaX*xnorm))/ \
            2/sigmaX2/sigmaY2/2.**(nuY/2)/gamma(nuY/2)
    return f

# this function definition is the integral, the convolution of a non-central chi-square with a central chi-square
def f_non_central_convolved(z, nuX, sigmaX2, deltaX, nuY, sigmaY2):
    args = (nuX, sigmaX2, deltaX, nuY, sigmaY2, z)
    return quad(f_non_central_joint, 0., z, args)[0]

# Test 1: evaluate the convolved density function
# Means and variances of example parameters suggest a run of Z = X + Y from about 700 to 1400.
zs = np.arange(40., 3040., 40.)
fzs = np.zeros(zs.shape)
for i in range(0, zs.shape[0]):
    fzs[i] = f_non_central_convolved(zs[i], nuX, sigmaX2, deltaX, nuY, sigmaY2)

# plot
import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(8,8))
ax1 = fig.add_subplot(311)

ax1.plot(zs, fzs, color='dodgerblue')
ax1.set_xlabel('X + Y')
ax1.set_ylabel('pdf')

plt.tight_layout()
plt.show()

# Test 2: Does the density integrate to unity?
args = (nuX, sigmaX2, deltaX, nuY, sigmaY2)
quad(f_non_central_convolved, 0., np.inf, args)[0]

0.999998814116672

# Test 3: The mean of this distribution should be (nuX + deltaX)*sigmaX2 + nuY*sigmaY2 = 936.29
(nuX + deltaX)*sigmaX2 + nuY*sigmaY2
args = (nuX, sigmaX2, deltaX, nuY, sigmaY2)
quad(lambda z: z*f_non_central_convolved(z, nuX, sigmaX2, deltaX, nuY, sigmaY2), 0., np.inf)[0]

/home/fred/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)

936.2867839948735

!jupyter nbconvert noncentralX2convolution.ipynb --to html

[NbConvertApp] Converting notebook noncentralX2convolution.ipynb to html
[NbConvertApp] Writing 280852 bytes to noncentralX2convolution.html

(nuX + deltaX)*sigmaX2 + nuY*sigmaY2

936.2900000000001