[SciPy-User] Skellam, fft and characteristic functions

josef.pktd@gmai... josef.pktd@gmai...
Wed Dec 30 11:10:44 CST 2009


This is just an example for some one-liners that it took me a long
time to figure out (because I don't know enough about fft, and I still
have one missing), related to the question I asked recently about
fourier transforms.

cfdiscrete_ex.txt are the results for the example file, cfdiscrete.py
has the actual functions (and more examples and tests.)

Essentially these are just some tools to work with (discrete)
distributions. I used it to verify results for the Skellam
distribution in the scipy trac ticket.

Josef
-------------- next part --------------
Using fast fourier transform to move between pmf/pdf and characteristic function
part 1: discrete distribution
================================================================================

(part 2: continuous distribution and compound poisson in preparation)


>>> import numpy as np
>>> from scipy import stats
>>> from cfdiscrete import *

#the last command imports
#__all__ = ['cf2cdfdiscrete', 'cf2cdfdiscretequad', 'cfn2pmf', 
#           'momentfromcfn', 'pmf2cf', 'Test_CFDiscrete', 'skellam']

# pmf2cf : pmf to characteristic function using fft
# cf2pmf : characteristic function to pmf using fft
# momentfromcfn : non-central moment from cf using fft differentiation
# cf2cdfdiscrete : cf to cdf using numerical integration with sum
# cf2cdfdiscretequad : cf to cdf using integrate.quad
# no fft version of cf2cdf, I didn't manage to figure out the integration

# cfn2pmf: characteristic function to probability mass function
# -------------------------------------------------------------


# Example: Skellam from scipy ticket


>>> def skellam_cf(w, mu1, mu2):
...     poisscf = stats.poisson._cf
...     return poisscf(w, mu1) * poisscf(-w, mu2)
... 
>>> mu1, mu2 = 10, 5
>>> k = np.arange(-5,20)
>>> print skellam.pmf(k, mu1, mu2)
[ 0.00324123  0.00633737  0.01155235  0.01960615  0.03094716  0.04540174
  0.06189433  0.07842461  0.09241881  0.10139793  0.10371928  0.09907658
  0.08854666  0.07418784  0.05839277  0.04326869  0.03024816  0.01999143
  0.01251688  0.00743899  0.00420459  0.00226421  0.00116371  0.00057179
  0.000269  ]


>>> skpmf, k2 = cfn2pmf(lambda w: skellam_cf(w, mu1, mu2))
>>> print skpmf[k.astype(int)]
[ 0.00324123  0.00633737  0.01155235  0.01960615  0.03094716  0.04540174
  0.06189433  0.07842461  0.09241881  0.10139793  0.10371928  0.09907658
  0.08854666  0.07418784  0.05839277  0.04326869  0.03024816  0.01999143
  0.01251688  0.00743899  0.00420459  0.00226421  0.00116371  0.00057179
  0.000269  ]
>>> print skellam.pmf(k, mu1, mu2) - skpmf[k.astype(int)]
[  1.47451495e-17   2.08166817e-17   2.77555756e-17   4.16333634e-17
   6.24500451e-17   7.63278329e-17   1.17961196e-16   1.38777878e-16
   1.38777878e-16   1.66533454e-16   1.52655666e-16   1.24900090e-16
   1.24900090e-16   6.93889390e-17   5.55111512e-17   1.38777878e-17
  -3.46944695e-18  -1.73472348e-17  -1.73472348e-17  -2.08166817e-17
  -1.82145965e-17  -2.34187669e-17  -1.97324795e-17  -2.03830008e-17
  -1.68051337e-17]


#Example 2: x+y-z , where x,y,z are Poisson distributed
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

#Similar to the story for Skellam but two teams against one,
#(independent rvs only)

>>> def poiss2versus1_cf(w, mu1, mu2, mu3):
...     poisscf = stats.poisson._cf
...     return poisscf(w, mu1) * poisscf(w, mu2) * poisscf(-w, mu3)
... 
>>> mu1a, mu2a, mu3a = 5, 5, 10
>>> cfn2v1 = lambda w: poiss2versus1_cf(w, mu1a, mu2a, mu3a)
>>> sk2pmf, k2 = cfn2pmf(cfn2v1)

>>> np.random.seed(9876345)
>>> nobs = 1000
>>> rvs = stats.poisson.rvs(mu1a, size=nobs) + stats.poisson.rvs(mu2a, size=nobs) -\
...       stats.poisson.rvs(mu3a
...       , size=nobs)
>>> rvsmin = rvs.min()

>>> freq = np.bincount(rvs-rvsmin)
>>> print (rvs == rvsmin).sum()
1
>>> kn = np.arange(rvsmin,rvs.max()+1).astype(int)


>>> print 'comparing sample frequency and theoretical frequency out of 1000'
comparing sample frequency and theoretical frequency out of 1000
>>> print np.column_stack((kn, freq, sk2pmf[kn]*1000)).astype(int)
[[-16   1   0]
 [-15   1   0]
 [-14   1   0]
 [-13   0   1]
 [-12   3   2]
 [-11   6   4]
 [-10  11   7]
 [ -9   5  11]
 [ -8  17  17]
 [ -7  22  25]
 [ -6  46  35]
 [ -5  43  47]
 [ -4  48  59]
 [ -3  83  71]
 [ -2  65  81]
 [ -1  93  87]
 [  0  91  89]
 [  1  91  87]
 [  2  75  81]
 [  3  66  71]
 [  4  62  59]
 [  5  60  47]
 [  6  29  35]
 [  7  42  25]
 [  8  13  17]
 [  9  14  11]
 [ 10   5   7]
 [ 11   5   4]
 [ 12   0   2]
 [ 13   2   1]]

#moments of the 2 versus 1 Poisson
# first seven uncenterd moments

>>> mncn = [momentfromcfn(i, cfn2v1) for i in range(7)]
>>> print mncn
[1.0, 3.1143029770981705e-014, 20.000000000022936, 2.0039346238461752e-009,
 1220.0000014978923, 0.0001034474399337382, 126020.08746095549]

>>> print 'mean sample %f, theory %f' % (rvs.mean(), mncn[1])
mean sample 0.069000, theory 0.000000
>>> print 'var sample %f, theory %f' % (rvs.var(), mncn[2]-mncn[1]**2)
var sample 20.088239, theory 20.000000


# cf2cdf: check accuracy
# ----------------------

#Poisson

>>> k = np.arange(50)
>>> cdfp = stats.poisson.cdf(k,mu1)

>>> cdfpiq = cf2cdfdiscretequad(k, lambda k:stats.poisson._cf(k,mu1))
>>> print np.max(np.abs(cdfp - cdfpiq))
1.90958360236e-014

>>> cdfpi = cf2cdfdiscrete(k, lambda k:stats.poisson._cf(k,mu1))
>>> print np.max(np.abs(cdfp - cdfpi))
2.02288774382e-012

#Skellam

>>> k = np.arange(-20,60)
>>> cdfs = skellam.cdf(k, mu1, mu2)

>>> cdfsiq = cf2cdfdiscretequad(k, lambda k_:skellam._cf(k_,mu1, mu2))
>>> print np.max(np.abs(cdfs - cdfsiq))
1.56812197777e-006

>>> cdfsi = cf2cdfdiscrete(k, lambda k_:skellam._cf(k_,mu1, mu2))
>>> print np.max(np.abs(cdfs - cdfsi))
1.56812197588e-006
>>> 

# skellam error is 1.568e-006 but that's maybe a problem with skellam._cdf
-------------- next part --------------



import numpy as np
from scipy import stats
from cfdiscrete import *

#the last command imports
#__all__ = ['cf2cdfdiscrete', 'cf2cdfdiscretequad', 'cfn2pmf', 
#           'momentfromcfn', 'pmf2cf', 'Test_CFDiscrete', 'skellam']

# pmf2cf : pmf to characteristic function using fft
# cf2pmf : characteristic function to pmf using fft
# momentfromcfn : non-central moment from cf using fft differentiation
# cf2cdfdiscrete : cf to cdf using numerical integration with sum
# cf2cdfdiscretequad : cf to cdf using integrate.quad
# no fft version of cf2cdf, I didn't manage to figure out the integration

# cfn2pmf: characteristic function to probability mass function
# -------------------------------------------------------------


# Example: Skellam from scipy ticket


def skellam_cf(w, mu1, mu2):
    # characteristic function
    poisscf = stats.poisson._cf
    return poisscf(w, mu1) * poisscf(-w, mu2)

mu1, mu2 = 10, 5
k = np.arange(-5,20)

#pmf of skellam (uses ncx2) and pmf by inversion of the characteristic function
print skellam.pmf(k, mu1, mu2)
skpmf, k2 = cfn2pmf(lambda w: skellam_cf(w, mu1, mu2))
#print skpmf[:10]
print skpmf[k.astype(int)]

#Error
print skellam.pmf(k, mu1, mu2) - skpmf[k.astype(int)]


#Example 2: x+y-z , where x,y,z are Poisson distributed
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

#Similar to the story for Skellam but two teams against one, (independent rvs only)

def poiss2versus1_cf(w, mu1, mu2, mu3):
    # characteristic function
    poisscf = stats.poisson._cf
    return poisscf(w, mu1) * poisscf(w, mu2) * poisscf(-w, mu3)

mu1a, mu2a, mu3a = 5, 5, 10
cfn2v1 = lambda w: poiss2versus1_cf(w, mu1a, mu2a, mu3a)
sk2pmf, k2 = cfn2pmf(cfn2v1)

np.random.seed(9876345)
nobs = 1000
rvs = stats.poisson.rvs(mu1a, size=nobs) + stats.poisson.rvs(mu2a, size=nobs) -\
      stats.poisson.rvs(mu3a
      , size=nobs)

rvsmin = rvs.min()
freq = np.bincount(rvs-rvsmin)
#check that bincount is correct
print (rvs == rvsmin).sum()
kn = np.arange(rvsmin,rvs.max()+1).astype(int)

print 'comparing sample frequency and theoretical frequency out of 1000'
print np.column_stack((kn, freq, sk2pmf[kn]*1000)).astype(int)

#moments of the 2 versus 1 Poisson
# first seven uncenterd moments
mncn = [momentfromcfn(i, cfn2v1) for i in range(7)]
print mncn
print 'mean sample %f, theory %f' % (rvs.mean(), mncn[1])
print 'var sample %f, theory %f' % (rvs.var(), mncn[2]-mncn[1]**2)



# cf2cdf: check accuracy
# ----------------------

k = np.arange(50)
cdfp = stats.poisson.cdf(k,mu1)
cdfpiq = cf2cdfdiscretequad(k, lambda k:stats.poisson._cf(k,mu1))
print np.max(np.abs(cdfp - cdfpiq))
cdfpi = cf2cdfdiscrete(k, lambda k:stats.poisson._cf(k,mu1))
print np.max(np.abs(cdfp - cdfpi))


k = np.arange(-20,60)
cdfs = skellam.cdf(k, mu1, mu2)
cdfsiq = cf2cdfdiscretequad(k, lambda k_:skellam._cf(k_,mu1, mu2))
print np.max(np.abs(cdfs - cdfsiq))
cdfsi = cf2cdfdiscrete(k, lambda k_:skellam._cf(k_,mu1, mu2))
print np.max(np.abs(cdfs - cdfsi))

# skellam error is 1.568e-006 but that's maybe a problem with skellam._cdf


-------------- next part --------------
# -*- coding: utf-8 -*-
"""
Created on Sun Dec 13 11:24:06 2009

Author: josef-pktd 

mostly written from scratch 
(based on definitions and with some trial and error for fft arguments)

cf2cdf is based on:

Numerical Inversion of a Characteristic Function
Author(s): R. B. Davies
Source: Biometrika, Vol. 60, No. 2 (Aug., 1973), pp. 415-417
Published by: Biometrika Trust
Stable URL: http://www.jstor.org/stable/2334555
Accessed: 13/11/2009 14:02


TODOs
-----

* create frozen distribution factory, with characteristic function, 
  npoints and bounds as inputs to constructor
* check consistency of npoints, make it consistently a choice variable
* convert arguments to use functions/distribution instead of arrays for cf,pmf,...
* possible extension choose error bounds instead of npoints for precision
* fit method, MLE or estimation based on empirical characteristic function(?)
  for MLE: fast way of getting loglike from cf for different parameters
* find characteristic functions for other discrete distributions
* add probability generating function, pdf, to distributions
* loc and scale (?): affine transformation of w in cf (check)
* helper functions for linear combinations of random variables (weighted sum)
  (following example skellam)
* connection to compound distribution (discrete mixture of distributions)
* are there any other useful tricks with fft ?
* try to get help with cf2cdf using fft
* check cf2cdf integration limits (0,pi) versus (-pi, pi). 
  Is real part always symmetric around zero?
  
* clean up analogous version for continuous distributions
* examples for sum of independent random variable
  try sum of log-normals
  for discrete: simple case X+Y-Z with each independently Poisson distributed
  (2 teams/players against 1 instead of 1on1 as in skellam

"""

import numpy as np
import scipy
from scipy import stats, integrate
import scipy.fftpack as fft

def maxabs(x,y):
    return np.max(np.abs(x-y))
    
def maxabsr(x,y):
    return np.max(np.abs(x.real-y.real))

def cfpoissmeth(self, w, mu):
    return np.exp(mu * (np.exp(1j * w) - 1.)) 
    
#patch poisson
stats.distributions.poisson_gen._cf = cfpoissmeth

# check
#print stats.poisson._cf(np.linspace(-0.5,0.5,6)[:,None],[1,10])



def pmf2cf(pmf, k):
    '''convert probability mass point to characteristic function
    
    checked only for positive k so far
    use pmf function instead of pmf array
    
    Parameters
    ----------
    pmf : array
        probabiltiy mass function
    k : array
        support points of random variable in array pmf
        required equidistant (? for fft)
        (integer, what if not integer values, e.g. 1.5?)
    
    Returns
    -------
    cf : array
        characteristic function
    w : array
        points at which characteristic function is evaluated
    
    '''
    npoints = len(k)
    w = k * 2.*np.pi/npoints  
    cf = fft.fft(pmf).conj() #no idea why I need conj to reverse sign of imag
    return cf, w


def cf2pmf(cfv, w):   
    '''convert characteristic function to probability mass point
    
    return only real part of inversion, maybe not correct for all distributions?
    checked only for positive k so far,
    ok for Skellam 
    use cf function instead of cf array
    
    Parameters
    ----------
    cfv : array
        characteristic function
    w : array
        points at which characteristic function is evaluated
        
    
    Returns
    -------
    pmf : array
        probabiltiy mass function
    k : array
        support points of random variable in array pmf
        not tested, what if there are negative integers
    
    '''
    k = w / 2./np.pi * npoints 
    pmf = fft.ifft(cfv.conj()).real #need conj to reverse sign of imag
    return pmf, k

def cfn2pmf(cfn, npoints=2**9):   
    '''convert characteristic function to probability mass point
    
    return only real part of inversion, maybe not correct for all distributions?
    checked only for positive k so far,
    ok for Skellam, but pmf for negative k are at end, recenter with fftshift
    use cf function instead of cf array
    
    Parameters
    ----------
    cfn : function
        characteristic function
    npoints : integer
        number of points for which characteristic function is evaluated
        
    
    Returns
    -------
    pmf : array
        probabiltiy mass function
    k : array
        support points of random variable in array pmf
        not tested, what if there are negative integers (those are at end)
    
    '''
    kw = np.arange(npoints)
    w = kw * 2.*np.pi/npoints  
    cfv = cfn(w)
    k = w / 2./np.pi * npoints 
    pmf = fft.ifft(cfv.conj()).real #need conj to reverse sign of imag
    return pmf, k

def cf2cdfdiscrete(k, cfn, npoints=2**9):
    '''inversion of characteristic function for cdf
    
    w needs to cover intervall (0, pi)
    
    uses simple fixed grid integration
    I didn't manage to figure out how to use fft for this
    
    Parameters
    ----------
    k array
        support points of random variable for which cdf is calculated
        (need not be equidistant, not related to precision)
    cfn : function
        characteristic function, only one argument allowed, `cfn(w)`
    npoints : integer
        number of points at which characteristic function is evaluated,
        determines the precision of the numerical integration
    
    Returns
    -------
    cdf : array
        cumulative distribution function for points of input k
        
    
    '''
    if np.ndim(k) == 1:
        k = k[:,None]
    k = k + 1  
    
    #uses symmetrie of real part around 0 to integrate over (0,pi]
    delta = np.pi/(npoints-1)
    w = (np.arange(npoints)+0.5)*delta
    #print k.shape, w.shape
    cf = cfn(w)
    #print w.min(), w.max()    
    cdf = 0.5-(cf/(1-np.exp(-1j*w))/2./np.pi
                 * np.exp(-1j*w*k)).real.sum(1)*(w[1]-w[0])*2
    return cdf


def cf2cdfdiscretequad(k, cfn):
    '''cf -> cdf inversion using quad
    
    from equation in section 3 of Davis (1973)
    
    Parameters
    ----------
    k : array
        support points of random variable for which cdf is calculated
        (need not be equidistant, not related to precision, calculations in loop)
    cfn : function
        characteristic function, only one argument allowed, `cfn(w)`
        
    
    Returns
    -------
    cdf : array
        cumulative distribution function for points of input k
        

    '''
#    def fnp(w):
#        return cfn(w)/(1-np.exp(-1j*w))/2./np.pi
#    
#    def fnintegA(m, A):
#        n = 1.*np.arange(len(A))
#        return (A*np.exp(-1j*2*np.pi*n/len(n)*k))
        
    def fn(w, k):
        w = max(w, 1e-10)
        A = cfn(w)/(1-np.exp(-1j*w))/2./np.pi
        A *= np.exp(-1j*w*k)
        return A.real
    
    cdfv = np.zeros(len(k))
    for ii,kk in enumerate(k):
        #cdfv[ii] = 0.5-integrate.quad(fn, -np.pi, np.pi, args=(kk+1,))[0]
        cdfv[ii] = 0.5-2*integrate.quad(fn, 0., np.pi, args=(kk+1,))[0]
    return cdfv


def momentfromcf(k, cf):
    '''calculates k-th non-central moment by differentiation of cf
    
    possible sign error
    
    Parameters
    ----------
    k : integer
        determines which moment is calculated.
    cf : array
        characteristic function values used by fft
    
    Returns
    -------
    non-central moment : float
    
    '''
    return np.abs(fft.diff(cf, k))[0]
    
    # without abs:
#    # no idea why odd even distinction
#    if k%2: # odd
#        return fft.diff(cf, k).imag[0]
#    else:
#        return fft.diff(cf, k).real[0]

def momentfromcfn(k, cfn, npoints = 2**9):
    '''calculates k-th non-central moment by differentiation of cf
    
    possible sign error
    
    Parameters
    ----------
    k : array
        support points of random variable for which cdf is calculated
        (need not be equidistant, not related to precision, calculations in loop)
    cfn : function
        characteristic function, only one argument allowed, `cfn(w)`
    npoints : integer
        number of points at which characteristic function is evaluated,
        determines the precision of the moment calculation by fft

        
    Returns
    -------
    non-central moment : float

    '''
    kw = np.arange(npoints)
    w = kw * 2.*np.pi/npoints  
    cfv = cfn(w)
    #print cfv[:10]
    return np.abs(fft.diff(cfv, k))[0]



######## from scipy track ticket, plus I added _cf, not the latest version

poisson = scipy.stats.distributions.poisson
ncx2 = scipy.stats.distributions.ncx2

# Skellam distribution 

class skellam_gen(scipy.stats.distributions.rv_discrete):
    def _rvs(self, mu1, mu2):
        n = self._size
        return np.random.poisson(mu1, n)-np.random.poisson(mu2, n)
    def _pmf(self, x, mu1, mu2):
        px = np.where(x < 0, ncx2.pdf(2*mu2, 2*(1-x), 2*mu1)*2,
                         ncx2.pdf(2*mu1, 2*(x+1), 2*mu2)*2)
        return px
    def _cdf(self, x, mu1, mu2):
        x = np.floor(x)
        px = np.where(x < 0, ncx2.cdf(2*mu2, -2*x, 2*mu1),
                         1-ncx2.cdf(2*mu1, 2*(x+1), 2*mu2))
        return px
        
    def _cf(self, w, mu1, mu2):
        # characteristic function
        poisscf = scipy.stats.distributions.poisson._cf
        return poisscf(w, mu1) * poisscf(-w, mu2)

    def _stats(self, mu1, mu2):
        mean = mu1 - mu2
        var = mu1 + mu2
        g1 = mean / np.sqrt((var)**3)
        g2 = 1 / var
        return mean, var, g1, g2
skellam = skellam_gen(a=-np.inf, name="skellam", longname='A Skellam',
                      shapes="mu1,mu2", extradoc="""

Skellam distribution

   Probability distribution of the difference of two correlated or
   uncorrelated Poisson random variables.

   If k1 and k2 are two Poisson variables with expected values lam1
   and lam2, then k1-k2 follows a Skellam distribution with parameters
   m1 = lam1 - rho*sqrt(lam1*lam2) and m2 = lam2 - rho*sqrt(lam1*lam2),
   rho being the correlation coefficient between k1 and k2.

   Parameters m1 and m2 must be strictly positive.

   For details see: http://en.wikipedia.org/wiki/Skellam_distribution
   
"""
                      )

###### end ticket



# and now for a test class

class Test_CFDiscrete(object):
    # class for testing, currently prints results
    def __init__(self, dist, args):
        self.dist = dist
        self.args = args
        
    def test_cf(self):
        # this fails for Skellam, need fftshift
        lambd = 10
        npoints = 2**8
        k = np.arange(npoints)
        w = k * 2.*np.pi/npoints  
        pmf = self.dist.pmf(k, *self.args)
        
        cfv = self.dist._cf(w, *self.args)
        
        k2 = np.arange(-npoints,npoints)
        w2 = k2 * 2.*np.pi/npoints  
        pmf = self.dist.pmf(k2, *self.args)
        #print 'nan in skellam.pmf', np.isnan(pmf).sum()
        pmf[np.isnan(pmf)] = 0.0
        #print 'check cdf for nans'
        #print np.isnan(self.dist.cdf(k2, lambd1, lambd2)).sum()
    
    
        #check error 5.95059410299e-015
        print 'cf error'
        cfi = pmf2cf(fft.fftshift(pmf),k2)[0][::2]
        ncf = min(len(cfv), len(cfi))  # check shape mismatch
        print maxabs(cfv[:ncf], cfi[:ncf])
    
    def test_pmf(self):
        npoints = 2**8
        k2 = np.arange(-npoints,npoints)
        w2 = k2 * 2.*np.pi/npoints  
        cfv = self.dist._cf(w, *self.args)
        print 'pmfi error'
        #Note: cf2pmf returns pmf of negative integers at end - use fftshift
        kcheck = np.arange(-20,60)
        print maxabs(fft.fftshift(pmf)[kcheck], cf2pmf(cfv, w2)[0][kcheck])
     
    def test_cdf(self):
        d=0
        #k = np.arange(d+0,d+10)
        k = np.arange(-20,60) #kcheck
        cdfi = cf2cdfdiscrete(k, lambda k_:self.dist._cf(k_,*self.args))
        cdf = self.dist.cdf(k,*self.args)
        #print 'cdfi'
        #print cdfi
        print 'cdfi error'  # 1.56812197671e-006
        print maxabs(cdf, cdfi)
        #print cdfi-skellam.cdf(k,lambd1, lambd2)
    
    def test_cdfquad(self):
        d=0
        k = np.arange(d+0,d+10)
        k = kcheck
        cdfi = cf2cdfdiscretequad(k, lambda k_:self.dist._cf(k_,*self.args))
        cdf = self.dist.cdf(k,*self.args)
        #print 'cdfi'
        #print cdfi
        print 'cdfi quad error'  # 1.56812197671e-006
        print maxabs(cdf, cdfi)
        #print cdfi-skellam.cdf(k,lambd1, lambd2)
        
    def test_moment(self):
        npoints = 2**8
        k = np.arange(npoints)
        w = k * 2.*np.pi/npoints
        cfv = self.dist._cf(w, *self.args)
        mean = momentfromcf(1, cfv)
        mnc2 = momentfromcf(2, cfv)
        m, v = self.dist.stats(*self.args)
        print 'mean', m, mean
        print 'var', v, mnc2 - mean**2
        
        
        for mi in range(7): # error in scipy.stats for mi>6
            mm = momentfromcf(mi, cfv)
            try:
                print mi, self.dist.moment(mi,*self.args), mm, mm-np.round(mm)
            except:  # catch all possible problems with scipy.stats.distributions
                print 'exception in moment calculation', mi, mm





#__all__ = [cf2cdfdiscrete, cf2cdfdiscretequad, cfn2pmf, 
#           momentfromcfn, pmf2cf, Test_CFDiscrete, skellam]
__all__ = ['cf2cdfdiscrete', 'cf2cdfdiscretequad', 'cfn2pmf', 
           'momentfromcfn', 'pmf2cf', 'Test_CFDiscrete', 'skellam']



if __name__ == '__main__':



    lambd = 10
    npoints = 2**8
    k = np.arange(npoints)
    w = k * 2.*np.pi/npoints  
    pmf = stats.poisson.pmf(k,lambd)

    cfv = stats.poisson._cf(w, lambd)

    #check error 5.95059410299e-015
    print 'cf error'
    print maxabs(cfv, pmf2cf(pmf,k)[0])

    #cf2pmf
    #------

    #check diff 5.96744875736e-016
    print 'pmfi error'
    print maxabs(pmf, cf2pmf(cfv, w)[0])


    #cf2cdfdiscrete
    #^^^^^^^^^^^^^^

    d = 0
    k2 = np.arange(d+0,d+10)

    cdfi = cf2cdfdiscrete(k2, lambda k:stats.poisson._cf(k,lambd))
    cdf = stats.poisson.cdf(k2,lambd)
    print 'cdfi'
    print cdfi
    print 'cdfi error'
    print maxabs(cdf, cdfi)
    print cdfi-stats.poisson.cdf(k2,lambd)


    #cf2cdfdiscretequad
    #^^^^^^^^^^^^^^^^^^

    d=0
    k2 = np.arange(d+0,d+10)
    cdfi2 = cf2cdfdiscretequad(k2, lambda k:stats.poisson._cf(k,lambd))
    cdf = stats.poisson.cdf(k2,lambd)
    print 'cdfi'
    print cdfi2
    print 'cdfi2 error'
    print maxabs(cdf, cdfi2)
    #print cdfi2-stats.poisson.cdf(k2,lambd)



    #momentfromcf momentfromcfn
    #^^^^^^^^^^^^^^^^^^^^^^^^^^

    #npoints = 2**8
    #k = np.arange(npoints)
    #w = k * 2.*np.pi/npoints  
    #pmf = stats.poisson.pmf(k,lambd)
    #cfv = stats.poisson._cf(w, lambd)


    mnc = [momentfromcf(i, cfv) for i in range(7)] 
    cfnpoiss = lambda w: stats.poisson._cf(w, lambd)
    mncn = [momentfromcfn(i, cfnpoiss) for i in range(7)] 


    for mi in range(7): # error in scipy.stats for mi>6
        mm = momentfromcf(mi, cfv)
        print mi, stats.poisson.moment(mi,lambd), mm, mm-np.round(mm)



    examples = ['skellam']

    if 'skellam' in examples:
        lambd1 = 10
        lambd2 = 5
        kcheck = np.arange(-20,60)
        
        lambd = 10
        npoints = 2**8
        k = np.arange(npoints)
        w = k * 2.*np.pi/npoints  
        pmf = skellam.pmf(k, lambd1, lambd2)
        
        cfv = skellam._cf(w, lambd1, lambd2)
        
        k2 = np.arange(-npoints,npoints)
        w2 = k2 * 2.*np.pi/npoints  
        pmf = skellam.pmf(k2, lambd1, lambd2)
        print 'nan in skellam.pmf', np.isnan(pmf).sum()
        pmf[np.isnan(pmf)] = 0.0
        print 'check cdf for nans'
        print np.isnan(skellam.cdf(k2, lambd1, lambd2)).sum()
        
        #check error 5.95059410299e-015
        print 'cf error'
        #cfi = pmf2cf(pmf,k2)[0]
        cfi = pmf2cf(fft.fftshift(pmf),k2)[0][::2] # I don't know why [::2]
        ncf = min(len(cfv), len(cfi))  # check shape mismatch
        print maxabs(cfv[:ncf], cfi[:ncf])
        
        print 'pmfi error'
        #Note: cf2pmf returns pmf of negative integers at end - use fftshift
        kcheck = np.arange(-20,60)
        print maxabs(fft.fftshift(pmf)[kcheck], cf2pmf(cfv, w2)[0][kcheck])
        
        d=0
        k = np.arange(d+0,d+10)
        k = kcheck
        cdfi = cf2cdfdiscrete(k, lambda k_:skellam._cf(k_,lambd1, lambd2))
        cdf = skellam.cdf(k,lambd1, lambd2)
        print 'cdfi'
        print cdfi
        print 'cdfi error'  # 1.56812197671e-006
        print maxabs(cdf, cdfi)
        print cdfi-skellam.cdf(k,lambd1, lambd2)


    #using the testclass
    #^^^^^^^^^^^^^^^^^^^

    print '\nTesting Skellam'
    td = Test_CFDiscrete(skellam, (lambd1, lambd2))

    td.test_cf()
    td.test_pmf()
    td.test_cdf()    
    # error is large 1.56812197671e-006, problem with skellam.cdf or cf2cdf ?
    td.test_cdfquad()  # why is this almost identical to td.test_cdf()
    td.test_moment()

    print '\nTesting Poisson using version including negative integers'
    td = Test_CFDiscrete(stats.poisson, (lambd1,))
    import time
    t0 = time.time()
    td.test_cf()
    t1 = time.time()
    print 'elapsed', t1-t0
    td.test_pmf()
    t2 = time.time()
    print 'elapsed', t2-t1
    td.test_cdf()
    t3 = time.time()
    print 'elapsed', t3-t2
    td.test_cdfquad()
    t4 = time.time()
    print 'elapsed', t4-t3
    td.test_moment()
    t5 = time.time()
    print 'elapsed', t5-t4

    #Note integrate.quad is very slow, up to 55 times array sum
    #times with fft don't register (0.0) without repetition



More information about the SciPy-User mailing list