Thread Subject: normal distribution pdfs and cdf: speed and accuracy

Dear All:

I am working on a Gibbs sampler with Markov Chain Monte Carlo.

In the marginal likelihood calculation, I compute a lot of normpdf's currently using the Lightspeed package.

Technically, what I would like to compute is $\Delta\cdot \rho$, that is
Delta * normpdf(a).

Very very technically what I am computing is:
normcdf(a+Delta/2)-normcdf(a-Delta/2)

which is approximately equal to
Delta * normpdf(a)

for the small values of Delta that I am using.

Now the normcdf computation is _significantly_ slower, and I'm not honestly sure since an expansion is used to approximate it (?) that it is actually more accurate than just calculating Delta * normpdf.

I would appreciate any input on this. Specifically:
(a) is it any more accurate in the case of small Deltas (say Delta = 0.05 or 0.05/3) to use the cdf here?
(b) if so, is there a faster implementation?

Thank you
Misha

"Misha Koshelev" <mk144210@bcm.edu> wrote in message <h09h17$obm$1@fred.mathworks.com>...
> Dear All:
>
> I am working on a Gibbs sampler with Markov Chain Monte Carlo.
>
> In the marginal likelihood calculation, I compute a lot of normpdf's currently using the Lightspeed package.
>
> Technically, what I would like to compute is $\Delta\cdot \rho$, that is
> Delta * normpdf(a).
>
> Very very technically what I am computing is:
> normcdf(a+Delta/2)-normcdf(a-Delta/2)
>
> which is approximately equal to
> Delta * normpdf(a)
>
> for the small values of Delta that I am using.
>
> Now the normcdf computation is _significantly_ slower, and I'm not honestly sure since an expansion is used to approximate it (?) that it is actually more accurate than just calculating Delta * normpdf.
>
> I would appreciate any input on this. Specifically:
> (a) is it any more accurate in the case of small Deltas (say Delta = 0.05 or 0.05/3) to use the cdf here
> (b) if so, is there a faster implementation?
>

cdf is integration of pdf. You might use a high order-quadrature rule of pdf on the interval a + (-Delta,+Delta)/2.

Bruno

"Misha Koshelev" <mk144210@bcm.edu> wrote in message <h09h17$obm$1@fred.mathworks.com>...
> Dear All:
>
> I am working on a Gibbs sampler with Markov Chain Monte Carlo.
>
> In the marginal likelihood calculation, I compute a lot of normpdf's currently using the Lightspeed package.
>
> Technically, what I would like to compute is $\Delta\cdot \rho$, that is
> Delta * normpdf(a).
>
> Very very technically what I am computing is:
> normcdf(a+Delta/2)-normcdf(a-Delta/2)
>
> which is approximately equal to
> Delta * normpdf(a)
>
> for the small values of Delta that I am using.
>
> Now the normcdf computation is _significantly_ slower, and I'm not honestly sure since an expansion is used to approximate it (?) that it is actually more accurate than just calculating Delta * normpdf.
>

It IS more accurate to use the cdf here. Perhaps you
should test that fact rather than worrying about it.
You will see that it is more accurate. Do a numerical
integration using quadgk to verify the result and the
accuracy obtained.


> I would appreciate any input on this. Specifically:
> (a) is it any more accurate in the case of small Deltas (say Delta = 0.05 or 0.05/3) to use the cdf here?

Yes, the cdf is a better approximation than your
simple rectangle rule approximation.

See above.


> (b) if so, is there a faster implementation?

The cdf implementation is even reasonably fast. But
your desire for speed is apparently more than your
CPU can satisfy. There will always be someone out
there who demands more speed. It matters not how
fast it runs, someone will say "Not fast enough!"

You can always use a poor approximation to the cdf
to gain speed. There are lots of them out there. In
fact you have chosen one that is quite poor, a
rectangle rule approximation. Depending on the
variance of your distribution, this may be entirely
adequate for your needs. It appears that you are
using a standard normal, so the variance == 1.
Only you know if it is good enough for you.

Compare the results and make the decision. Only
you can make that tradeoff between speed and the
accuracy you will lose. Remember to use

format long g

to display the different results.

John

> Now the normcdf computation is _significantly_ slower, and I'm not
> honestly sure since an expansion is used to approximate it (?) that it is
> actually more accurate than just calculating Delta * normpdf.

Misha, look inside normcdf and you'll see that it's basically a call to
erfc(), with a small amount of error checking and other overhead. You might
want to call erfc yourself directly, and see if that's any faster.

-- Tom

"Tom Lane" <tlane@mathworks.com> wrote in message <h0bujk$bpt$1@fred.mathworks.com>...
> > Now the normcdf computation is _significantly_ slower, and I'm not
> > honestly sure since an expansion is used to approximate it (?) that it is
> > actually more accurate than just calculating Delta * normpdf.
>
> Misha, look inside normcdf and you'll see that it's basically a call to
> erfc(), with a small amount of error checking and other overhead. You might
> want to call erfc yourself directly, and see if that's any faster.
>
> -- Tom
>

Thank you all for your valuable advice. I will look into calling the erfc function myself.

Misha

Misha Koshelev <mk144210@bcm.edu> wrote:
> Dear All:
>
> I am working on a Gibbs sampler with Markov Chain Monte Carlo.
>
> In the marginal likelihood calculation, I compute a lot of normpdf's
currently using the Lightspeed package.
>
> Technically, what I would like to compute is $\Delta\cdot \rho$, that is
> Delta * normpdf(a).
>
> Very very technically what I am computing is:
> normcdf(a+Delta/2)-normcdf(a-Delta/2)
>
> which is approximately equal to
> Delta * normpdf(a)
>
> for the small values of Delta that I am using.
>
> Now the normcdf computation is _significantly_ slower, and I'm not
honestly sure since an expansion is used to approximate it (?) that it is
actually more accurate than just calculating Delta * normpdf.
>
> I would appreciate any input on this. Specifically:
> (a) is it any more accurate in the case of small Deltas (say Delta = 0.05
or 0.05/3) to use the cdf here?

See John's response in this.


> (b) if so, is there a faster implementation?
>

 * Marsaglia, George "Evaluating the Normal Distribution",
 * Journal of Statistical Software 11, 4 (July 2004).
 * http://www.jstatsoft.org/

You will need a few lines of C code, but as always with Marsaglia's stuff,
this is fast and excellent, and worth the effort required to understand
what he is doing.


--
Dr Tristram J. Scott
Energy Consultant

tristram.scott@ntlworld.com (Tristram Scott) wrote in message <Nx3Xl.14512$2t7.3733@newsfe16.ams2>...
> Misha Koshelev <mk144210@bcm.edu> wrote:
> > Dear All:
> >
> > I am working on a Gibbs sampler with Markov Chain Monte Carlo.
> >
> > In the marginal likelihood calculation, I compute a lot of normpdf's
> currently using the Lightspeed package.
> >
> > Technically, what I would like to compute is $\Delta\cdot \rho$, that is
> > Delta * normpdf(a).
> >
> > Very very technically what I am computing is:
> > normcdf(a+Delta/2)-normcdf(a-Delta/2)
> >
> > which is approximately equal to
> > Delta * normpdf(a)
> >
> > for the small values of Delta that I am using.
> >
> > Now the normcdf computation is _significantly_ slower, and I'm not
> honestly sure since an expansion is used to approximate it (?) that it is
> actually more accurate than just calculating Delta * normpdf.
> >
> > I would appreciate any input on this. Specifically:
> > (a) is it any more accurate in the case of small Deltas (say Delta = 0.05
> or 0.05/3) to use the cdf here?
>
> See John's response in this.
>
>
> > (b) if so, is there a faster implementation?
> >
>
> * Marsaglia, George "Evaluating the Normal Distribution",
> * Journal of Statistical Software 11, 4 (July 2004).
> * http://www.jstatsoft.org/
>
> You will need a few lines of C code, but as always with Marsaglia's stuff,
> this is fast and excellent, and worth the effort required to understand
> what he is doing.
>
>
> --
> Dr Tristram J. Scott
> Energy Consultant

Thank you. I made a c function. Once it is approved I will post link here.

Misha

"Misha Koshelev" <mk144210@bcm.edu> wrote in message <h0kfl6$b70$1@fred.mathworks.com>...
> tristram.scott@ntlworld.com (Tristram Scott) wrote in message <Nx3Xl.14512$2t7.3733@newsfe16.ams2>...
> > Misha Koshelev <mk144210@bcm.edu> wrote:
> > > Dear All:
> > >
> > > I am working on a Gibbs sampler with Markov Chain Monte Carlo.
> > >
> > > In the marginal likelihood calculation, I compute a lot of normpdf's
> > currently using the Lightspeed package.
> > >
> > > Technically, what I would like to compute is $\Delta\cdot \rho$, that is
> > > Delta * normpdf(a).
> > >
> > > Very very technically what I am computing is:
> > > normcdf(a+Delta/2)-normcdf(a-Delta/2)
> > >
> > > which is approximately equal to
> > > Delta * normpdf(a)
> > >
> > > for the small values of Delta that I am using.
> > >
> > > Now the normcdf computation is _significantly_ slower, and I'm not
> > honestly sure since an expansion is used to approximate it (?) that it is
> > actually more accurate than just calculating Delta * normpdf.
> > >
> > > I would appreciate any input on this. Specifically:
> > > (a) is it any more accurate in the case of small Deltas (say Delta = 0.05
> > or 0.05/3) to use the cdf here?
> >
> > See John's response in this.
> >
> >
> > > (b) if so, is there a faster implementation?
> > >
> >
> > * Marsaglia, George "Evaluating the Normal Distribution",
> > * Journal of Statistical Software 11, 4 (July 2004).
> > * http://www.jstatsoft.org/
> >
> > You will need a few lines of C code, but as always with Marsaglia's stuff,
> > this is fast and excellent, and worth the effort required to understand
> > what he is doing.
> >
> >
> > --
> > Dr Tristram J. Scott
> > Energy Consultant
>
> Thank you. I made a c function. Once it is approved I will post link here.
>
> Misha

Here it is:
http://www.mathworks.com/matlabcentral/fileexchange/24373

Its faster than calling erfc (I used the profile function) and seems to be approx as accurate for my uses (I get a difference over the entire log marginal likelihood of about 10^-12th of the value, which does not seem too bad).

Thank you all again.

Misha