Characteristic and distribution functions

Given a continuous rv X, its characteristic function (CF) is defined as \(\phi(t)=E(e^{itX})\).
- MFG: \(M(t)=\phi(-it)\)
- Exp dist \(\phi(t)=(1-it/\lambda)^{-1}\)
- \(\chi_k^2\) dist: \(\phi(t)=(1-2it)^{-k/2}\)
- Gamma(k,\(\theta\)) dist: \(\phi(t)=(1-it\theta)^{-k}\)
Note the CF of sum of independent rvs is just the product of individual CFs.
Given \(\phi(t)\), we can compute dist function
- PDF: \(f(x) = \frac{1}{2\pi}\int_{R} e^{-itx}\phi(t)dt\)
- CDF: \(F(x) = \frac{1}{2}-\frac{1}{\pi}\int_0^{\infty} \frac{Im[e^{-itx}\phi(t)]}{t}dt\).

A numerical example for sum of independent Exp rvs

We consider \(\sum_{i=1}^n Exp(\lambda_i)\)
- if all \(\lambda_i\) are equal, we will have a Gamma dist.
- if any two \(\lambda_i\) are different, we have the hypoexponential dist.
- For both scenarios, we have analytical forms of PDF.

## Compute PDF from characteristic function based on the FFT
## phi: characteristic function
## n: Number of points, ideally a power of 2
## ab: evaluate density on interval ab
CF2PDF <- function(phi, n=2^10, ab){
  a = ab[1]; b = ab[2]
  i = 0:(n-1);  dx = (b-a)/n
  x = a + i * dx;  dt = 2*pi / ( n * dx )
  c1 = -n/2 * dt;  d1 =  n/2 * dt
  t = c1 + i * dt; phi.t = phi(t)
  X = exp(-(0+1i)*i*dt*a )*phi.t
  Y = fft(X)  # FFT
  density = dt/(2*pi)*exp(-(0+1i)*c1*x)*Y
  res = data.frame(CF= phi.t, x=x, density=Re(density) )
  return(res)
}
## Compute the Characteristic function for sum of independent Exp rvs
## x: input values
## lambda: coefficients for the weighted sums
CFExpSums <- function(x, lambda){
  sapply(x, function(t){
    prod( (1-(0+1i)*t/lambda)^(-1) )
  })
}
## PDF for hypoexponential dist
dExpSums <- function(x, lambda){
  K = length(lambda)
  Cik = rep(0, K)
  for(i in 1:K){
    Cik[i] = prod( lambda[-i]/(lambda[-i]-lambda[i]) )
  }
  res = Cik[1]*dexp(x,lambda[1])
  for(i in 2:K) res = res + Cik[i]*dexp(x,lambda[i])
  return(res)
}
library(distr)
## Direct calculation of convolution of probability distributions
## Efficient implementation based on the FFT.
DistExpSums <- function(lambda){
  K = length(lambda)
  Ex = Exp(rate=lambda[1])
  for(i in 2:K) Ex = Ex + Exp(rate=lambda[i])
  return(Ex)
}
## PDF comparison
distroptions(DefaultNrGridPoints=2^15)
ab = c(8,25)
lambda1 = rep(1,10)
d1 = CF2PDF(function(x) CFExpSums(x,lambda1), n=2^15, ab)
Ex1 = DistExpSums(lambda1)
y1 = d(Ex1)(d1$x)
lambda2 = 1:8
d2 = CF2PDF(function(x) CFExpSums(x,lambda2), n=2^15, ab)
Ex2 = DistExpSums(lambda2)
y2 = d(Ex2)(d2$x)
## Gamma
par(mfrow=c(1,2))
plot(d1$x, d1$density, las=1, log='y', ylim=c(1e-3,0.15))
points(d1$x, y1, col=3, pch=3, cex=0.5)
curve(dgamma(x,length(lambda1),1), add=TRUE, col=2, lwd=2)
## Hypoexp
plot(d2$x, d2$density, las=1, log='y', ylim=c(1e-14,0.5))
points(d2$x, y2, col=3, pch=3, cex=0.5)
curve(dExpSums(x,lambda2), add=TRUE, col=2, lwd=2)

### CDF comparison
ab = c(0,25); n = 2^15
lambda1 = rep(1,10)
d1 = CF2PDF(function(x) CFExpSums(x,lambda1), n=n, ab)
pr1 = (cumsum(d1$density)-d1$density[1]/2-d1$density/2)*(d1$x[2]-d1$x[1])
Ex1 = DistExpSums(lambda1)
y1 = p(Ex1)(d1$x)
lambda2 = 1:8
d2 = CF2PDF(function(x) CFExpSums(x,lambda2), n=2^15, ab)
pr2 = (cumsum(d2$density)-d2$density[1]/2-d2$density/2)*(d2$x[2]-d2$x[1])
Ex2 = DistExpSums(lambda2)
y2 = p(Ex2)(d2$x)
dy2 = dExpSums(d2$x,lambda2)
ip2 = (cumsum(dy2)-dy2[1]/2-dy2/2)*(d2$x[2]-d2$x[1])

## Gamma
par(mfrow=c(2,2))
plot(d1$x, pr1, las=1, xlab='x', ylab='CDF')
points(d1$x, y1, col=3, pch=3, cex=0.5)
curve(pgamma(x,length(lambda1),1), add=TRUE, col=2, lwd=2)
## Hypoexp
plot(d2$x, pr2, las=1, xlab='x', ylab='CDF')
points(d2$x, y2, col=3, pch=3, cex=0.5)
lines(d2$x, ip2, col=2, lwd=2)

i1 = which((d1$x>15)&(d1$x<20)); i2 = which((d2$x>8)&(d2$x<12)) 
plot(d1$x[i1], pr1[i1], las=1, xlab='x', ylab='CDF')
points(d1$x[i1], y1[i1], col=3, pch=3, cex=0.5)
curve(pgamma(x,length(lambda1),1), add=TRUE, col=2, lwd=2)
## Hypoexp
plot(d2$x[i2], pr2[i2], las=1, xlab='x', ylab='CDF')
points(d2$x[i2], y2[i2], col=3, pch=3, cex=0.5)
lines(d2$x[i2], ip2[i2], col=2, lwd=2)