Why is my application of Gaussian quadrature to the Gaussian curve so off?

Question

I am using $52$ points and weights extracted from this website and simply pasting them on an Excel spreadsheet, trying to reproduce the value of the exponential part of a standard normal distribution $\exp(-1/2 x^2)$ between $x=-1$ and $x=1,$ which using Wolfram alpha turns out to be $1.711.$ By calculating the equation on the points in the list, and weighing them by the corresponding values the sum is $2.389.$

Save for a bad oversight, this is widely off the mark. What am I missing?

Answer 1

Not worth an answer, really, but still mind-blowingly precise... When done right! The mistake is simple: It turns out that the spreadsheet I was using interprets the input =B2*exp(-C2^2/2) as something different to $e^{-\frac{x^2}2}$. When entering the query in the cell as =B2*exp(-0.5*C2^2) the approximation is uncanny:

Even more amazing is that with just $2$ nodes the value is so close to the Wolfram Alpha computation of $$1.7112487837842976063466092405600901224528285701829945206404685631941007769337068369063650663943116305674373309410539791803060979978578671006175305064740831375299713796653168516207212436314557192881717234968718536490769308263847260997086025142320547263060192774029898329388574373167823550401261985881379404339028250455018105953910652604458332469390533651006390219427744143838301439717850503901152978005$$ at $1.692963449781230,$ already matching the first $5$ decimal places with just $5$ nodes $(1.711249393524460).$

---

Same goes obviously if I contrast the standard normal (Gaussian) distribution with from the "tables" at say $2$ standard deviations, i.e. $0.977249868051821$ with the Legendre-Gaussian quadrature numerical integration say from $x=-10$ (taken as "minus infinity) to $x=2.$

With a simple linear transformation of the limits of integration (originally in LG quadrature from $-1$ to $1$) via the system $-10=m(-1) + c$ and $2= m(1) +c,$ the slope of the transformation is $m=11,$ and the constant $c=-9$ with the Jacobian $dt = 11 dx.$ Therefore $x=11t-9.$ Since the pdf to integrate is $\frac{1}{\sqrt 2}e^{\frac 1 2 x^2},$ the actual Excel formula for the first row (there is a header, so the first row is the second line): =(1/(sqrt(2 *PI())))*B2*11*exp(-0.5*(11*C2-9)^2) or $\frac{1}{\sqrt{2\pi}}w_2 11 \exp(-\frac 1 2(x(t))^2)\mathrm dt$ and the result with $52$ nodes / weights: $0.9772498681.$ Amazing!

A bit surprising is that actually doing the transformation $x= 2 + \frac{t-1}{t+1},$ which would tend to $-\infty$ as $t\to-1,$ and would tend to $2$ as $t\to 1,$ and logically introducing the formula as $\frac{1}{\sqrt{2\pi}}w_2 \left(\frac{2}{(t+1)^2}\right)\exp\left(-\frac 1 2\left(x= 2 + \frac{t-1}{t+1}\right)^2\right)\mathrm dt$ results in a less approximate value of $0.9772498669.$

This can also be solved with a simple Monte Carlo simulation: mean(dnorm(runif(5e5,-10,2))) * (2-(-10)) which is the average of $5$ million values of the pdf of the standard normal calculated for a corresponding $5$ million random samples (uniformly distributed from $-10$ standard deviations to $2$ standard deviations) times the interval from max to min $x$ value, since for a given interval $[a,b]$ the mean of the function $\hat f(x)=\frac{1}{b-a}\int_a^b f(x) dx.$

Alternatively, just seeing the proportion of random points that fall below the curve in randomly generated points in the rectangle inscribing the part of the curve to integrate works:

set.seed(2019)
n = 1e6
minx = -10
maxx = 2
miny = 0
maxy = dnorm(0)
xax = runif(n,minx,maxx)
yax = runif(n,miny,maxy)
true = dnorm(xax)
mean(yax < true) * (maxy-miny) * (maxx-minx)
# 0.978367644225595 # Estimated value
pnorm(2)
# 0.977249868051821 # True value