pdf.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

#include "jwsc/base/limits.h"
#include "jwsc/math/constants.h"
#include "jwsc/math/slatec.h"
#include "jwsc/prob/pdf.h"


#define  SQRT_2_PI_SINGLE      2.5066282746310004f

#define  SQRT_2_PI_DOUBLE      2.50662827463100050241576528481103

#define  INV_SQRT_2_PI_SINGLE  0.3989422804014326f

#define  INV_SQRT_2_PI_DOUBLE  0.39894228040143267793994605993438

#define  LOG_SQRT_2_PI_SINGLE  0.9189385332046727

#define  LOG_SQRT_2_PI_DOUBLE  .91893853320467274178032973640561


float standard_gaussian_pdf_f(float x)
{
    return INV_SQRT_2_PI_SINGLE * exp(-0.5 * x * x);
}

double standard_gaussian_pdf_d(double x)
{
    return INV_SQRT_2_PI_DOUBLE * exp(-0.5 * x * x);
}

float gaussian_pdf_f(float mu, float sigma, float x)
{
    assert(sigma > 1.0e-8);

    return INV_SQRT_2_PI_SINGLE * (1 / sigma) * 
        exp(-0.5 * (x - mu)*(x - mu) / (sigma*sigma));
}

double gaussian_pdf_d(double mu, double sigma, double x)
{
    assert(sigma > 1.0e-16);

    return INV_SQRT_2_PI_DOUBLE * (1 / sigma) * 
        exp(-0.5 * (x - mu)*(x - mu) / (sigma*sigma));
}

float log_standard_gaussian_pdf_f(float x)
{
    return - LOG_SQRT_2_PI_SINGLE - 0.5*x*x;
}

double log_standard_gaussian_pdf_d(double x)
{
    return - LOG_SQRT_2_PI_DOUBLE - 0.5*x*x;
}

float log_gaussian_pdf_f(float mu, float sigma, float x)
{
    float log_sigma;

    assert(sigma > 1.0e-8);

    log_sigma = log(sigma);

    return - LOG_SQRT_2_PI_SINGLE - log_sigma - 
        (0.5*(x - mu)*(x - mu)) / (sigma*sigma);
}

double log_gaussian_pdf_d(double mu, double sigma, double x)
{
    double log_sigma;

    assert(sigma > 1.0e-16);

    log_sigma = log(sigma);

    return - LOG_SQRT_2_PI_DOUBLE - log_sigma - 
        (0.5*(x - mu)*(x - mu)) / (sigma*sigma);
}

#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF

float integrate_gaussian_pdf_f
(
    float mu, 
    float sigma, 
    float a, 
    float b
)
{
    float x_1;
    float x_2;
    float div;
    float result;

    div = sqrt(2) * sigma;

    x_1 = (mu - b) / div;
    x_2 = (mu - a) / div;

#if defined JWSC_HAVE_SLATEC
    result = 0.5 * (slatec_derf(x_2) - slatec_derf(x_1));
#elif defined JWSC_HAVE_ERF
    result = 0.5 * (erf(x_2) - erf(x_1));
#endif

    return result;
}

double integrate_gaussian_pdf_d
(
    double mu, 
    double sigma, 
    double a, 
    double b
)
{
    double x_1;
    double x_2;
    double div;
    double result;

    div = sqrt(2) * sigma;

    x_1 = (mu - b) / div;
    x_2 = (mu - a) / div;

#if defined JWSC_HAVE_SLATEC
    result = 0.5 * (slatec_derf(x_2) - slatec_derf(x_1));
#elif defined JWSC_HAVE_ERF
    result = 0.5 * (erf(x_2) - erf(x_1));
#endif

    return result;
}

#endif


float sample_standard_gaussian_pdf_f(float a, float b)
{
    int accept;
    float s, l, u;

    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_f(a, b);
        l = log_standard_gaussian_pdf_f(s);
        u = log(sample_uniform_pdf_f(0, 1));

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}

double sample_standard_gaussian_pdf_d(double a, double b)
{
    int accept;
    double s, l, u;

    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_d(a, b);
        l = log_standard_gaussian_pdf_d(s);
        u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, 1));

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}

float sample_gaussian_pdf_f
(
    float mu, 
    float sigma, 
    float a, 
    float b
)
{
    int accept;
    float s, l, u;
    float max;

    max = gaussian_pdf_f(mu, sigma, mu);
    assert(JWSC_MIN_LOG_ARG <= max);
    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_f(a, b);
        l = log_gaussian_pdf_f(mu, sigma, s);
        u = log(sample_uniform_pdf_f(JWSC_MIN_LOG_ARG, max));

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}

double sample_gaussian_pdf_d
(
    double mu, 
    double sigma, 
    double a, 
    double b
)
{
    int accept;
    double s, l, u;
    double max;

    max = gaussian_pdf_d(mu, sigma, mu);
    assert(JWSC_MIN_LOG_ARG <= max);
    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_d(a, b);
        l = log_gaussian_pdf_d(mu, sigma, s);
        u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, max));

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}

#if defined JWSC_HAVE_TGAMMA
extern double tgamma(double x);

float gamma_pdf_f(float alpha, float beta, float x)
{
    assert(x > 1.0e-8);
    assert(alpha > 1.0e-8);
    assert(beta > 1.0e-8);

    return pow(beta, alpha) / tgamma(alpha) * pow(x, alpha - 1) * exp(-beta*x);
}

double gamma_pdf_d(double alpha, double beta, double x)
{
    assert(x > 1.0e-16);
    assert(alpha > 1.0e-16);
    assert(beta > 1.0e-16);

    return pow(beta, alpha) / tgamma(alpha) * pow(x, alpha - 1) * exp(-beta*x);
}

#endif


#if defined JWSC_HAVE_LGAMMA

float log_gamma_pdf_f(float alpha, float beta, float x)
{
    assert(x > 1.0e-8);
    assert(alpha > 1.0e-8);
    assert(beta > 1.0e-8);

    return alpha*log(beta) - lgamma(alpha) + (alpha - 1)*log(x) - beta*x;
}

double log_gamma_pdf_d(double alpha, double beta, double x)
{
    assert(x > 1.0e-16);
    assert(alpha > 1.0e-16);
    assert(beta > 1.0e-16);

    return alpha*log(beta) - lgamma(alpha) + (alpha - 1)*log(x) - beta*x;
}

#endif


#if defined JWSC_HAVE_TGAMMA

float sample_gamma_pdf_f
(
    float alpha, 
    float beta, 
    float a, 
    float b
)
{
    int accept;
    float s, l, u;
    float max;


    if(alpha < 20 && beta < 20)
    {
        max = (alpha > 1) ? gamma_pdf_f(alpha, beta, (alpha - 1)/beta) : 1;
    }
    else
    {
        // for large alpha and beta, intermediate computations are 
        // more robust in the log space.
        max = exp(log_gamma_pdf_f(alpha, beta, (alpha - 1)/beta));
    }

    assert(JWSC_MIN_LOG_ARG <= max);
    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_f(a, b);
        l = log_gamma_pdf_f(alpha, beta, s);
        u = log(sample_uniform_pdf_f(JWSC_MIN_LOG_ARG, max));

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}

double sample_gamma_pdf_d
(
    double alpha, 
    double beta, 
    double a, 
    double b
)
{
    int accept;
    double s, l, u;
    double max;

    if(alpha < 20 && beta < 20)
    {
        max = (alpha > 1) ? gamma_pdf_d(alpha, beta, (alpha - 1)/beta) : 1;
    }
    else
    {
        // for large alpha and beta, intermediate computations are 
        // more robust in the log space.
        max = exp(log_gamma_pdf_d(alpha, beta, (alpha - 1)/beta));
    }

    assert(JWSC_MIN_LOG_ARG <= max);
    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_d(a, b);
        l = log_gamma_pdf_d(alpha, beta, s);
        u = log(sample_uniform_pdf_d(JWSC_MIN_LOG_ARG, max));

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}

#endif


float sample_uniform_pdf_f(float a, float b)
{
    if (a == b)
    {
        return a;
    }

    return (b - a)*(rand()/(float)RAND_MAX) + a;
}

double sample_uniform_pdf_d(double a, double b)
{
    if (a == b)
    {
        return a;
    }

    return (b - a)*(rand()/(double)RAND_MAX) + a;
}

float sample_sine_pdf_f()
{
    int accept;
    float s, l, u;

    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_f(0, JWSC_PI_2);
        l = sin(s);
        u = sample_uniform_pdf_f(0, 1.0);

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}

double sample_sine_pdf_d()
{
    int accept;
    double s, l, u;

    accept = 0;

    while (!accept)
    {
        s = sample_uniform_pdf_d(0, JWSC_PI_2);
        l = sin(s);
        u = sample_uniform_pdf_d(0, 1.0);

        if (u <= l)
        {
            accept = 1;
        }
    }

    return s;
}