Skip to content

Instantly share code, notes, and snippets.

@andyfaff
Created February 6, 2018 22:54
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save andyfaff/c8a94ab87333c42791e58e7d6e12a05e to your computer and use it in GitHub Desktop.
Save andyfaff/c8a94ab87333c42791e58e7d6e12a05e to your computer and use it in GitHub Desktop.
quad + stats.norm

import numpy as np from scipy.stats import norm from scipy.integrate import quad xs = np.linspace(0, 40, 40) [quad(norm.pdf, -np.inf, x)[0] for x in xs] [0.4999999999999999, 0.8474695889892241, 0.9798802563169332, 0.9989542536697613, 0.9999795701700213, 0.9999998537411905, 0.99999999962187, 0.9999999999996502, 1.0000000000000115, 1.0000000000013474, 1.0000000000000004, 1.0000000000000002, 1.0, 1.0000000000000016, 1.0000000000000002, 1.0000000000000002, 0.9999999999999999, 1.0000000000000002, 1.0000000000000004, 1.0, 1.0000000000000002, 1.8562200694967622e-10, 1.0000000000000004, 1.0000000000000007, 1.000000000001782, 1.0000000000000002, 1.0000000000000002, 1.0, 1.0, 1.0000000000000004, 1.0000000000000004, 0.9999999999999996, 0.9999999999999998, 0.9999999999999996, 1.0000000000000058, 2.4991668789183437e-11, 8.295218115538186e-15, 9.616318458943375e-19, 3.893483499692327e-23, 5.505747798875221e-28]

@chrisb83
Copy link

chrisb83 commented Feb 8, 2018

You have the same problem in R, I think both use QUADPACK. It seems like it is a known problem in general, see e.g. https://stat.ethz.ch/R-manual/R-devel/library/stats/html/integrate.html, it mentions problems when you integrate norm.pdf from 0 to 20000.

I am not aware that you can change the subdivisions of the integration interval.

Would one of the lines below help?
[0.5 + quad(norm.pdf, 0, x)[0] for x in xs]
[1 - quad(norm.pdf, x, np.inf)[0] for x in xs]

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment