How to Calculate the Characteristic Function

Key Insights

Characteristic functions uniquely identify probability distributions and transform convolutions into simple multiplications, making them invaluable for analyzing sums of random variables
Analytical solutions exist for common distributions (normal, exponential, Poisson), but numerical methods using scipy.integrate or FFT handle arbitrary distributions effectively
Numerical stability requires careful selection of the evaluation range and handling of complex exponentials, particularly for heavy-tailed distributions

Introduction to Characteristic Functions

The characteristic function is the Fourier transform of a probability distribution. While moment generating functions get more attention in introductory courses, characteristic functions are more fundamental—they exist for all probability distributions, even those with undefined moments.

You’ll need characteristic functions when analyzing sums of independent random variables, identifying distributions from their properties, or working with heavy-tailed distributions where moment generating functions don’t exist. In signal processing and time series analysis, they’re essential for understanding the frequency domain representation of random processes.

The key advantage: if X and Y are independent random variables, the characteristic function of X + Y is simply the product of their individual characteristic functions. This transforms a difficult convolution integral into straightforward multiplication.

Mathematical Foundation

The characteristic function of a random variable X is defined as:

φ_X(t) = E[e^(itX)] = E[cos(tX) + i·sin(tX)]

For continuous distributions with PDF f(x):

φ_X(t) = ∫_{-∞}^{∞} e^(itx) f(x) dx

For discrete distributions with PMF p(x):

φ_X(t) = Σ e^(itx) p(x)

Key properties every characteristic function satisfies:

φ(0) = 1 (always)
|φ(t)| ≤ 1 for all t
φ(-t) = φ̄(t) (conjugate symmetry for real-valued random variables)
Continuity at t = 0

Let’s visualize this relationship for a standard normal distribution:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Standard normal PDF
x = np.linspace(-4, 4, 1000)
pdf = norm.pdf(x)

# Characteristic function: φ(t) = exp(-t²/2) for standard normal
t = np.linspace(-3, 3, 1000)
char_func = np.exp(-t**2 / 2)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

ax1.plot(x, pdf, 'b-', linewidth=2)
ax1.set_title('Probability Density Function')
ax1.set_xlabel('x')
ax1.set_ylabel('f(x)')
ax1.grid(True, alpha=0.3)

ax2.plot(t, char_func, 'r-', linewidth=2)
ax2.set_title('Characteristic Function (Real Part)')
ax2.set_xlabel('t')
ax2.set_ylabel('φ(t)')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Analytical Calculation Methods

For many standard distributions, you can derive the characteristic function analytically through direct integration or summation.

Exponential Distribution: For X ~ Exp(λ) with PDF f(x) = λe^(-λx) for x ≥ 0:

φ(t) = ∫₀^∞ e^(itx) λe^(-λx) dx = λ ∫₀^∞ e^((it-λ)x) dx = λ/(λ - it)

Normal Distribution: For X ~ N(μ, σ²):

φ(t) = exp(iμt - σ²t²/2)

Poisson Distribution: For X ~ Poisson(λ):

φ(t) = Σ_{k=0}^∞ e^(itk) · (e^(-λ)λ^k)/k! = exp(λ(e^(it) - 1))

Here’s a Python implementation for the exponential distribution with verification:

import numpy as np
from scipy.integrate import quad

def exponential_char_func_analytical(t, lambda_param):
    """Analytical characteristic function for Exp(λ)"""
    return lambda_param / (lambda_param - 1j * t)

def exponential_char_func_numerical(t, lambda_param):
    """Numerical verification using integration"""
    def integrand_real(x):
        return np.cos(t * x) * lambda_param * np.exp(-lambda_param * x)
    
    def integrand_imag(x):
        return np.sin(t * x) * lambda_param * np.exp(-lambda_param * x)
    
    real_part, _ = quad(integrand_real, 0, np.inf)
    imag_part, _ = quad(integrand_imag, 0, np.inf)
    
    return real_part + 1j * imag_part

# Verification
lambda_param = 2.0
t_values = np.linspace(-5, 5, 11)

print("t\t\tAnalytical\t\tNumerical\t\tDifference")
print("-" * 70)
for t in t_values:
    analytical = exponential_char_func_analytical(t, lambda_param)
    numerical = exponential_char_func_numerical(t, lambda_param)
    diff = np.abs(analytical - numerical)
    print(f"{t:6.2f}\t\t{analytical:.6f}\t{numerical:.6f}\t{diff:.2e}")

Numerical Computation Approaches

When analytical solutions aren’t available or you’re working with empirical data, numerical methods become necessary.

For distributions with known PDFs, use numerical integration:

from scipy.integrate import quad
import numpy as np

def characteristic_function_continuous(t, pdf, x_range=(-10, 10)):
    """
    Compute characteristic function numerically for continuous distribution
    
    Args:
        t: Point at which to evaluate
        pdf: Probability density function
        x_range: Integration range tuple
    """
    def integrand_real(x):
        return np.cos(t * x) * pdf(x)
    
    def integrand_imag(x):
        return np.sin(t * x) * pdf(x)
    
    real_part, real_err = quad(integrand_real, x_range[0], x_range[1])
    imag_part, imag_err = quad(integrand_imag, x_range[0], x_range[1])
    
    return complex(real_part, imag_part), (real_err, imag_err)

# Example: Gamma distribution
from scipy.stats import gamma

shape, scale = 2.0, 2.0
gamma_pdf = lambda x: gamma.pdf(x, shape, scale=scale)

t_vals = np.linspace(-2, 2, 50)
char_vals = [characteristic_function_continuous(t, gamma_pdf, (0, 50))[0] 
             for t in t_vals]

For empirical distributions from sample data, FFT provides an efficient approach:

def characteristic_function_empirical(data, t_values):
    """
    Compute characteristic function from sample data using FFT
    
    Args:
        data: Sample data array
        t_values: Array of t values at which to evaluate
    """
    n = len(data)
    char_func = np.zeros(len(t_values), dtype=complex)
    
    for i, t in enumerate(t_values):
        # Direct computation: average of e^(itX)
        char_func[i] = np.mean(np.exp(1j * t * data))
    
    return char_func

# Example with normal samples
np.random.seed(42)
samples = np.random.normal(0, 1, 10000)

t_range = np.linspace(-3, 3, 100)
empirical_cf = characteristic_function_empirical(samples, t_range)

# Compare with theoretical
theoretical_cf = np.exp(-t_range**2 / 2)

plt.figure(figsize=(10, 5))
plt.plot(t_range, np.real(empirical_cf), 'b-', label='Empirical (Real)', alpha=0.7)
plt.plot(t_range, theoretical_cf, 'r--', label='Theoretical', linewidth=2)
plt.xlabel('t')
plt.ylabel('φ(t)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title('Empirical vs Theoretical Characteristic Function')
plt.show()

Practical Applications

The most powerful application is computing the distribution of sums. If X and Y are independent with characteristic functions φ_X and φ_Y, then:

φ_{X+Y}(t) = φ_X(t) · φ_Y(t)

def sum_of_distributions_example():
    """
    Demonstrate using characteristic functions to find 
    distribution of sum of independent RVs
    """
    # X ~ Exp(1), Y ~ Exp(2)
    lambda1, lambda2 = 1.0, 2.0
    
    def char_func_sum(t):
        cf_x = lambda1 / (lambda1 - 1j * t)
        cf_y = lambda2 / (lambda2 - 1j * t)
        return cf_x * cf_y
    
    # Verify with simulation
    np.random.seed(42)
    n_samples = 100000
    x_samples = np.random.exponential(1/lambda1, n_samples)
    y_samples = np.random.exponential(1/lambda2, n_samples)
    sum_samples = x_samples + y_samples
    
    t_vals = np.linspace(-2, 2, 50)
    
    # Characteristic function from theory
    cf_theoretical = np.array([char_func_sum(t) for t in t_vals])
    
    # Characteristic function from samples
    cf_empirical = characteristic_function_empirical(sum_samples, t_vals)
    
    # Plot comparison
    plt.figure(figsize=(10, 5))
    plt.plot(t_vals, np.real(cf_theoretical), 'r-', 
             label='Theoretical (Product)', linewidth=2)
    plt.plot(t_vals, np.real(cf_empirical), 'b--', 
             label='Empirical (Sum of Samples)', alpha=0.7)
    plt.xlabel('t')
    plt.ylabel('Re(φ(t))')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.title('Sum of Independent Exponentials')
    plt.show()

sum_of_distributions_example()

Implementation Best Practices

A production-ready implementation requires careful attention to numerical stability and edge cases:

class CharacteristicFunction:
    """Robust characteristic function calculator"""
    
    def __init__(self, pdf=None, data=None, support=None):
        """
        Initialize with either PDF or empirical data
        
        Args:
            pdf: Callable PDF function
            data: Empirical data array
            support: Tuple (min, max) for PDF support
        """
        if pdf is None and data is None:
            raise ValueError("Must provide either pdf or data")
        
        self.pdf = pdf
        self.data = data
        self.support = support if support else (-np.inf, np.inf)
    
    def evaluate(self, t, method='auto'):
        """
        Evaluate characteristic function at t
        
        Args:
            t: Point(s) at which to evaluate
            method: 'analytical', 'numerical', or 'auto'
        """
        t = np.atleast_1d(t)
        
        if self.data is not None:
            return self._evaluate_empirical(t)
        else:
            return self._evaluate_from_pdf(t)
    
    def _evaluate_empirical(self, t):
        """Compute from empirical data"""
        result = np.zeros(len(t), dtype=complex)
        for i, t_val in enumerate(t):
            if np.abs(t_val) < 1e-10:
                result[i] = 1.0 + 0j
            else:
                result[i] = np.mean(np.exp(1j * t_val * self.data))
        return result
    
    def _evaluate_from_pdf(self, t):
        """Compute from PDF using numerical integration"""
        result = np.zeros(len(t), dtype=complex)
        
        # Determine appropriate integration range
        if np.isfinite(self.support[0]) and np.isfinite(self.support[1]):
            x_min, x_max = self.support
        else:
            x_min, x_max = -100, 100  # Default range
        
        for i, t_val in enumerate(t):
            if np.abs(t_val) < 1e-10:
                result[i] = 1.0 + 0j
            else:
                real_part, _ = quad(
                    lambda x: np.cos(t_val * x) * self.pdf(x),
                    x_min, x_max
                )
                imag_part, _ = quad(
                    lambda x: np.sin(t_val * x) * self.pdf(x),
                    x_min, x_max
                )
                result[i] = complex(real_part, imag_part)
        
        return result
    
    def validate(self, t_range=(-5, 5), n_points=50):
        """Validate characteristic function properties"""
        t_vals = np.linspace(t_range[0], t_range[1], n_points)
        cf_vals = self.evaluate(t_vals)
        
        # Check φ(0) = 1
        phi_0 = self.evaluate(np.array([0.0]))[0]
        assert np.abs(phi_0 - 1.0) < 1e-6, f"φ(0) = {phi_0}, expected 1.0"
        
        # Check |φ(t)| ≤ 1
        magnitudes = np.abs(cf_vals)
        assert np.all(magnitudes <= 1.01), f"Max |φ(t)| = {magnitudes.max()}"
        
        print("Validation passed!")
        return True

# Example usage
from scipy.stats import norm

cf_normal = CharacteristicFunction(
    pdf=lambda x: norm.pdf(x, 0, 1),
    support=(-10, 10)
)

cf_normal.validate()
t_test = np.linspace(-3, 3, 100)
values = cf_normal.evaluate(t_test)

Conclusion

Characteristic functions provide a powerful tool for probability theory and statistical analysis. Use analytical formulas when available—they’re exact and computationally cheap. For standard distributions (normal, exponential, gamma, Poisson), memorize or reference the closed-form expressions.

When working with custom distributions or empirical data, numerical integration via scipy.integrate.quad offers good accuracy for moderate dimensions. For large sample sizes, the empirical approach (direct averaging of e^(itX)) is both simple and effective.

The key decision points: If you have a known distribution with a standard form, use the analytical formula. If you have a custom PDF, use numerical integration. If you only have sample data, use the empirical estimator. Always validate that φ(0) = 1 and |φ(t)| ≤ 1 to catch implementation errors.

For production systems, wrap your implementation in a class with proper error handling, validate inputs, and choose integration ranges carefully based on the distribution’s support. The code examples here provide a solid foundation for most practical applications.