MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

New to MATLAB?

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

Subject: normal distribution pdfs and cdf: speed and accuracy

From: Misha Koshelev

### Misha Koshelev (view profile)

Date: 4 Jun, 2009 22:15:03

Message: 1 of 8

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

Subject: normal distribution pdfs and cdf: speed and accuracy

From: Bruno Luong

### Bruno Luong (view profile)

Date: 5 Jun, 2009 06:18:01

Message: 2 of 8

"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

Subject: normal distribution pdfs and cdf: speed and accuracy

From: John D'Errico

### John D'Errico (view profile)

Date: 5 Jun, 2009 09:15:03

Message: 3 of 8

"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

 Subject: normal distribution pdfs and cdf: speed and accuracy From: Tom Lane Date: 5 Jun, 2009 20:19:00 Message: 4 of 8 > 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

Subject: normal distribution pdfs and cdf: speed and accuracy

From: Misha Koshelev

### Misha Koshelev (view profile)

Date: 5 Jun, 2009 21:22:01

Message: 5 of 8

"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

 Subject: normal distribution pdfs and cdf: speed and accuracy Date: 8 Jun, 2009 07:59:41 Message: 6 of 8 Misha Koshelev 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

Subject: normal distribution pdfs and cdf: speed and accuracy

From: Misha Koshelev

### Misha Koshelev (view profile)

Date: 9 Jun, 2009 01:59:02

Message: 7 of 8

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 Subject: normal distribution pdfs and cdf: speed and accuracy From: Misha Koshelev ### Misha Koshelev (view profile) Date: 9 Jun, 2009 13:31:01 Message: 8 of 8 "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

## Tags for this Thread

No tags are associated with this thread.

### What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.