math.h - xcode - cdf

[dirtypeanut]

Hi everybody!

I am new in C++ and I am trying to find out how to use the cumulative distribution function for normal random variables. I am using the math.h in Xcode.

Any help would be really appreciated as I am really stuck.

Best,
Peanut

[Stewbond]

I don't think there is a cumulative distribution function in <math.h> or <cmath>.

This is a list of everything in there
http://cplusplus.com/reference/clibrary/cmath/

[chipp]

do you mean like updated version of cmath?

[Stewbond]

If you're looking for the CDF of a normal or gaussian distribution you can use this:

#include <iostream>
#include <cmath>
using namespace std;

// Returns the erf() of a value (not super precice, but ok)
double erf(double x)
{  
 double y = 1.0 / ( 1.0 + 0.3275911 * x);   
 return 1 - (((((
        + 1.061405429  * y
        - 1.453152027) * y
        + 1.421413741) * y
        - 0.284496736) * y 
        + 0.254829592) * y) 
        * exp (-x * x);      
}

// Returns the probability of x, given the distribution described by mu and sigma.
double pdf(double x, double mu, double sigma)
{
  //Constants
  static const double pi = 3.14159265; 
  return exp( -1 * (x - mu) * (x - mu) / (2 * sigma * sigma)) / (sigma * sqrt(2 * pi));
}

// Returns the probability of [-inf,x] of a gaussian distribution
double cdf(double x, double mu, double sigma)
{
	return 0.5 * (1 + erf((x - mu) / (sigma * sqrt(2.))));
}

int main()
{
	
	double x, mu, sigma;
	cout << "x, mu, sigma: ";
	cin >> x >> mu >> sigma;

	cout << "PDF of x is: " << pdf(x,mu,sigma) << endl;
	cout << "CDF of x is: " << cdf(x,mu,sigma) << endl;

	return 0;
}

[Stewbond]

Okay, I decided to have some fun and make a generic function to calculate the integral of f(x) from -inf to x.

In this code, you can customize the pdf function to whatever you like and the CDF will still calculate it. The only restriction is that the pdf function must accept x plus two other parameters.

#include <iostream>
#include <cmath>
using namespace std;

// Returns the probability of x, given the distribution described by mu and sigma.
double pdf(double x, double mu, double sigma)
{
	static const double pi = 3.14159265; 
	return exp( -1 * (x - mu) * (x - mu) / (2 * sigma * sigma)) / (sigma * sqrt(2 * pi));
}

// Returns the integral from -inf to x of any function that accepts x and 2 other parameters
double cdf(double x, double arg1, double arg2, double(*pPDF)(double,double,double))
{

	double sum = 0;
	double ninf = -1e3; // Negative infinity, just use something small
	double n = 1e7; // The number of "buckets" that we'll calculate, more is more accurate but takes more time;

	for (double k = 1.; k < n-1; k++)
		sum += pPDF( x + k*(x-ninf)/n ,arg1,arg2);

	return ((x - ninf) / n) * ((pPDF(x,arg1,arg2) + pPDF(ninf,arg1,arg2))/2 + sum);
}

int main()
{
	double x, mu, sigma;
	cout << "x, mu, sigma: ";
	cin >> x >> mu >> sigma;

	cout << "PDF of x is: " << pdf(x,mu,sigma) << endl;
	cout << "CDF of x is: " << cdf(x,mu,sigma,pdf) << endl; // Note we added a pointer to your pdf function

	return 0;
}

[Stewbond]

Okay one more. I've made it a little more c++ and optimized it a bit.

Based on this generic function which will approximate an integral of f(x) between a and b:

double integral(double a, double b, double(*f)(double))
{
	double sum = 0;
	double n = 1e5; // Tune for speed/accuracy
	double c = (b-a)/n;

	for (double k = 1.; k < n-1; k+=1.)
		sum += f(a + k*c);

	return c * ((f(a) + f(b)) / 2 + sum);
}

We can make a class with pdf and cdf functions:

const double pi = 3.14159265;

class NormalDistribution
{
public:
    NormalDistribution(double _mu, double _sigma) : mu(_mu), sigma(_sigma) {}
    inline double pdf(double x);
    double cdf(double x);
private:
    double mu;
    double sigma;
};

inline double NormalDistribution::pdf(double x)
{  //You will change this function if you want another type of distribution
    return exp( -1 * (x - mu) * (x - mu) / (2 * sigma * sigma)) / (sigma * sqrt(2 * pi));
}

double NormalDistribution::cdf(double x)
{
    // Integral from a to x;
    const double ninf = mu - 10 * sigma; // Dependent on the type of distribution, tune as appropriate
    double sum = 0;
    double n = 1e5; // tune for speed/accuracy
    double c = (x - ninf) / n;

    for (double k = 1.; k < n-1; k++)
        sum += pdf( ninf + k*c);

    return c * ((pdf(x) + pdf(ninf))/2 + sum);
}

Now to use it is easy:

    NormalDistribution Gauss(0. , 1.); // Make a normal distribution (mu = 0, sigma = 1)
    double x = 0.;                     // This is what we'll get the cdf of.

    std::cout << "Integral of f(x) from a to b is: " << Gauss.cdf(x) << std::endl;