Binomial distribution in q

thesweeheng's weblog

Recently I was using SciPy’s scipy.stats.binom.pmf(x,n,p). I though it would be great if I could have such a function in q. So a simple idea is to construct a binomial tree with probabilities attached. Recalling that a Pascal triangle is generated using n{0+':x,0}1, I modified it to get:

q)pmf:{[n;p]n{(0,y*1-x)+x*y,0}[p]/1#1f}
q)pmf[6;0.3]
0.000729 0.010206 0.059535 0.18522 0.324135 0.302526 0.117649
q)sum pmf[1000;0.3]
1f

What is great about this method is that it is stable. Compared to SciPy 0.7.0, it was more accurate too (it is a known issue that older SciPy has buggy binom.pmf):

>>> scipy.stats.binom.pmf(range(0,41),40,0.3)[-5:]
>>> array([3.33066907e-15, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.11022302e-16])
q)-5#pmf[40;0.7]
3.293487e-15 1.52594e-16 5.162955e-18 1.134715e-19 1.215767e-21

Unfortunately this method is too slow for large n. For large n, we need more sophisticated methods. For the interested reader, take a look at Catherine Loader’s Fast and Accurate Computation of Binomial Probabilities paper and an implementation of a binomial distribution in…

View original post 1 more word

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s