1

I'm in the middle of a code translation from Matlab to C++, and for some important reasons I must obtain the cumulative distribution function of a 'normal' function (in matlab, 'norm') with mean=0 and variance=1.

The implementation in Matlab is something like this:

map.c = cdf( 'norm', map.c, 0,1 );

Which is suppose to be the equalization of the histogram from map.c.

The problem comes at the moment of translating it into C++, due to the lack of decimals I have. I tried a lot of typical cdf implementations: such as the C++ code I found here, Cumulative Normal Distribution Function in C/C++ but I got an important lack of decimals, so I tried the boost implementation:

#include "boost/math/distributions.hpp" 

boost::math::normal_distribution<> d(0,1);

but still it is not the same implementation as Matlab's (I guess it seems to be even more precise!)

Does anyone know where I could find the original Matlab source for such a process, or which is the correct amount of decimals I should consider?

Thanks in advance!

SecretAgentMan
  • 2,856
  • 7
  • 21
  • 41
Ciller
  • 13
  • 6
  • If the boost one is more precise, shouldn't it be better for you? – Shahbaz Jun 04 '12 at 15:24
  • The thing is, I already have some results with this kind of Matlab's CDF which I should campare to, so I should use the same algorithm... :( – Ciller Jun 04 '12 at 15:45
  • "Lack of decimals" is an imprecise phrase. double-precision float is the standard precision used by MATLAB, and is the same precision you're using in C++ with the 'double' type. What results make you say that you're having a precision problem? Is it possible you're printing them without enough precision, leading you to believe there's not enough internal precision? – Peter Jun 04 '12 at 18:42
  • Maybe I just didn't explain it properly enough... I just need to implement the Matlab's normal CDF in C++. I tried boost normal distribution, but I don't know why (maybe its because of the cdf's error table too...) I'm not getting the same results that I obtain with Matlab :D – Ciller Jun 04 '12 at 20:07
  • The SO post you cited had comments saying the `pow` in the code needed to be changed to Horner's. Did you do that? Basically, the last expression in the computation of `w` can be changed into `((((a5 * K + a4) * K + a3) * K + a2) * K + a1) * K` instead. – jxh Jun 05 '12 at 00:33
  • No, no, I avoided using that function you're talking about and tried to use Boost's instead ;) Anyway, I'll keep on working on it :D! – Ciller Jun 05 '12 at 08:35
  • you have to tell boost to use double precision ! re read the most recent docs carefully – pyCthon Sep 10 '12 at 20:05

3 Answers3

1

Octave is an open-source Matlab clone. Here's the source code for Octave's implementation of normcdf: http://octave-nan.sourcearchive.com/documentation/1.0.6/normcdf_8m-source.html

It should be (almost) the same as Matlab's, if it helps you.

Daniel Hanrahan
  • 4,801
  • 7
  • 31
  • 44
1

The Gaussian CDF is an interesting function. I do not know whether my answer will interest you, but it is likely to interest others who look up your question later, so here it is.

One can compute the CDF by integrating the Taylor series of the PDF term by term. This approach works well in the body of the Gaussian bell curve, only it fails numerically in the tails. In the tails, it wants special-function techniques. The best source I have read for this is N. N. Lebedev's Special Functions and Their Applications, Ch. 2, Dover, 1972.

thb
  • 13,796
  • 3
  • 40
  • 68
0

C and C++ support long double for a more precise floating point type. You could try using that in your implementation. You can check your compiler documentation to see if it provides an even higher precision floating point type. GCC 4.3 and higher provide __float128 which has even more precision.

Shahbaz
  • 46,337
  • 19
  • 116
  • 182
jxh
  • 69,070
  • 8
  • 110
  • 193