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

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

 Subject: normal distribution pdfs and cdf: speed and accuracy From: Misha Koshelev 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 Date: 5 Jun, 2009 06:18:01 Message: 2 of 8 "Misha Koshelev" wrote in message ... > 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 Date: 5 Jun, 2009 09:15:03 Message: 3 of 8 "Misha Koshelev" wrote in message ... > 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 Date: 5 Jun, 2009 21:22:01 Message: 5 of 8 "Tom Lane" wrote in message ... > > 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 Date: 9 Jun, 2009 01:59:02 Message: 7 of 8 tristram.scott@ntlworld.com (Tristram Scott) wrote in message ... > 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 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 Date: 9 Jun, 2009 13:31:01 Message: 8 of 8 "Misha Koshelev" wrote in message ... > tristram.scott@ntlworld.com (Tristram Scott) wrote in message ... > > 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 > > 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