This notebook derives the GEC (Giddings–Eyring–Carmichael) monopore characteristic function from first principles using SymPy.
Goal: Show how the CF emerges from:
Exponential egress time distribution
Poisson number of adsorptions
Independence assumption
from sympy import symbols, exp, factorial, simplify, I, oo, Sum
import sympy as sp
# Define symbols
w, n, t_bar, n_bar = symbols('w n t_bar n_bar', real=True, positive=True)# CF of single exponential random variable
phi_tau = 1 / (1 - I*w*t_bar)
print("CF of single egress time:")
phi_tauCF of single egress time:
# CF of sum of n i.i.d. exponentials
phi_sum_n = phi_tau**n
print("CF of sum of n egress times:")
phi_sum_nCF of sum of n egress times:
Step 4: Law of Total Expectation (Requires Independence!)¶
For random sum where :
Critical: This step requires independence between and !
# Poisson probability mass function
P_n = exp(-n_bar) * n_bar**n / factorial(n)
print("Poisson PMF P(r_M = n):")
P_nPoisson PMF P(r_M = n):
Step 5: Evaluate the Sum (Poisson Generating Function)¶
We need to evaluate:
This is the probability generating function of Poisson:
# Compound Poisson CF: sum over all possible n values
# φ_tS(ω) = Σ P(n) * [φ_τ(ω)]^n
summand = P_n * phi_sum_n
print("Summand P(r_M=n) * [φ_τ(ω)]^n:")
summandSummand P(r_M=n) * [φ_τ(ω)]^n:
# The sum Σ (e^(-n_bar) * n_bar^n / n!) * z^n = exp(n_bar*(z-1))
# We substitute z = phi_tau
z = symbols('z')
# Poisson generating function
pgf = exp(n_bar * (z - 1))
print("Poisson generating function E[z^r_M]:")
display(pgf)
# Substitute z = phi_tau to get GEC CF
gec_cf_derived = pgf.subs(z, phi_tau)
print("\nGEC CF (substitute z = φ_τ(ω)):")
gec_cf_derivedPoisson generating function E[z^r_M]:
GEC CF (substitute z = φ_τ(ω)):
# Simplify the exponent
exponent = simplify(phi_tau - 1)
print("Exponent φ_τ(ω) - 1:")
display(exponent)
# Full GEC CF
gec_cf_final = simplify(gec_cf_derived)
print("\nFinal GEC monopore CF:")
gec_cf_finalExponent φ_τ(ω) - 1:
Final GEC monopore CF:
from sympy import symbols, exp, simplify, I, cancel
w = symbols('w')
n1, t1 = symbols('n1, t1', positive=True, real=True)
# Original form from Dondi Eq. 43
gec_monopore_cf_original = exp(n1*(1/(1 - I*w*t1) - 1))
print("Original form (Dondi Eq. 43):")
display(gec_monopore_cf_original)
# Manually simplify the exponent to show equivalence
# 1/(1 - iωτ) - 1 = [1 - (1 - iωτ)]/(1 - iωτ) = iωτ/(1 - iωτ)
exponent_original = n1*(1/(1 - I*w*t1) - 1)
exponent_simplified = simplify(exponent_original)
print("\nOriginal exponent:")
display(exponent_original)
print("\nSimplified exponent:")
display(exponent_simplified)
# This should be: n1*I*w*t1/(1 - I*w*t1)
# Verify by manual calculation
numerator = 1 - (1 - I*w*t1)
denominator = 1 - I*w*t1
manual_simplified = n1 * numerator / denominator
manual_simplified = simplify(manual_simplified)
print("\nManual simplification [1 - (1 - iωτ)]/(1 - iωτ):")
display(manual_simplified)
# Alternative form
gec_alternative = exp(n1*I*w*t1/(1 - I*w*t1))
print("\nAlternative form exp[n̄·iω·τ̄/(1 - iω·τ̄)]:")
display(gec_alternative)
# Verify they're the same by checking the difference
print("\nDifference between original and alternative:")
display(simplify(gec_monopore_cf_original - gec_alternative))
print("\n✓ Both forms are mathematically identical!")
print(" Form 1: exp[n̄(1/(1-iωτ̄) - 1)]")
print(" Form 2: exp[n̄·iωτ̄/(1-iωτ̄)]")Original form (Dondi Eq. 43):
Original exponent:
Simplified exponent:
Manual simplification [1 - (1 - iωτ)]/(1 - iωτ):
Alternative form exp[n̄·iω·τ̄/(1 - iω·τ̄)]:
Difference between original and alternative:
✓ Both forms are mathematically identical!
Form 1: exp[n̄(1/(1-iωτ̄) - 1)]
Form 2: exp[n̄·iωτ̄/(1-iωτ̄)]
Moments from Characteristic Functions¶
The characteristic function encodes all statistical properties of a distribution. We can extract moments using derivatives:
Raw Moments (about the origin)¶
The -th raw moment is , obtained via:
: Mean (first raw moment)
: Second raw moment
Central Moments (about the mean)¶
The -th central moment is where :
: Variance
: Skewness component (normalized: )
: Kurtosis component (normalized: )
Method: Central moments can be computed from the shifted CF:
This is implemented below by dividing the CF by before differentiation.
How to use SymPy to symbolically compute moments¶
We can extract statistical moments from the characteristic function using SymPy as follows.
from sympy import symbols, exp, diff, simplify, I
w = symbols('w')
def raw_moment(cf, k):
return simplify((-I)**k * diff(cf, w, k).subs(dict(w=0)))
def central_moment(cf, k):
m = raw_moment(cf, 1)
return simplify(raw_moment(cf/exp(I*w*m),k))n1, t1 = symbols('n1, t1')
gec_monopore_cf = exp(n1*(1/(1 - I*w*t1) - 1))
gec_monopore_cfraw_moment(gec_monopore_cf, 1)central_moment(gec_monopore_cf, 2)central_moment(gec_monopore_cf, 3)