[SciPy-User] Skellam, fft and characteristic functions
josef.pktd at gmail.com
josef.pktd at gmail.com
Wed Dec 30 12:10:44 EST 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