import pylab
import scipy.stats as ss
nrm = ss.norm
nx = nrm.pdf
# Lookup table based implementations ------------------------------------------
def init_nx_table(xlim=5,N=1001):
"""Go from -xlim to +xlim and make N entries, return us the dx and x0.
if N is made odd it is better"""
idx0 = int(N/2)
tbl = pylab.zeros(N)
x = pylab.linspace(-xlim,xlim,N)
dx = x[1] - x[0]
tbl = nx(x)
return x, tbl, idx0, dx
def nx_lookup(x,mu,tbl, idx0, dx):
"""x needs to be an array."""
sz = tbl.size
ret = pylab.zeros(x.size) #Our results
idx = (x-mu)/dx + idx0 + .5 #indexes into our table need +.5 because of rounding
idxidx = pylab.find((idx>=0) & (idx<sz)) #indexes of valid indexes
ret[idxidx] = tbl[idx[idxidx].astype('int16')]
return ret
def testnx(dotiming=False):
xtbl, tbl, idx0, dx = init_nx_table()
if dotiming:
import cProfile
x = pylab.linspace(-10,10,1000000)
cProfile.runctx('nx(x)',globals(),locals())
cProfile.runctx('nx_lookup(x, 0, tbl, idx0, dx)',globals(),locals())
cProfile.runctx('pylab.interp(x, xtbl, tbl, left=0, right=0)',globals(),locals())
else:
x = pylab.linspace(-10,10,1000)
x0 = nx(x)
x1 = nx_lookup(x, 0, tbl, idx0, dx)
x2 = pylab.interp(x, xtbl, tbl, left=0, right=0)
pylab.plot(x, x0-x1)
pylab.plot(x, x0-x2)
pylab.ylabel('Error')
if __name__ == "__main__":
testnx()
testnx(True)
