pmf.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 <math.h>
#include <assert.h>
#include <inttypes.h>

#include "jwsc/base/limits.h"
#include "jwsc/prob/pdf.h"
#include "jwsc/prob/pmf.h"


#define  INV_E_SINGLE       0.3678794411714423f

#define  INV_E_DOUBLE       0.36787944117144232159552377016146

#define  ALMOST_ONE_SINGLE  0.9999999999999999

#define  ALMOST_ONE_DOUBLE  0.99999999999999999999999999999999


float poisson_pmf_f(float lambda, uint32_t n)
{
    float p;
    int i;

    if (lambda <= 0 && lambda == n)
    {
        return 1;
    }
    else if (lambda <= 0)
    {
        return 0;
    }

    p = 1;

    if (lambda < n)
    {
        for (i = n; i > (n - lambda); i--)
        {
            p *= lambda * INV_E_SINGLE * (1 / (float)i);
        }

        for (; i > 0; i--)
        {
            p *= lambda * (1 / (float)i);
        }
    }
    else
    {
        for (i = n; i > 0; i--)
        {
            p *= lambda * INV_E_SINGLE * (1 / (float)i);
        }
        
        p *= exp(-(lambda - n));
    }

#if defined(JWSC_HAVE_ISNAN) && defined(JWSC_HAVE_ISINF)
    if (isnan(p) || isinf(p))
    {
        return 0;
    }
#endif

    return p;
}

double poisson_pmf_d(double lambda, uint32_t n)
{
    double p;
    int i;

    if (lambda <= 0 && lambda == n)
    {
        return 1;
    }
    else if (lambda <= 0)
    {
        return 0;
    }

    p = 1;

    if (lambda < n)
    {
        for (i = n; i > (n - lambda); i--)
        {
            p *= lambda * INV_E_DOUBLE * (1 / (double)i);
        }

        for (; i > 0; i--)
        {
            p *= lambda * (1 / (double)i);
        }
    }
    else
    {
        for (i = n; i > 0; i--)
        {
            p *= lambda * INV_E_DOUBLE * (1 / (double)i);
        }
        
        p *= exp(-(lambda - n));
    }

#if defined(JWSC_HAVE_ISNAN) && defined(JWSC_HAVE_ISINF)
    if (isnan(p) || isinf(p))
    {
        return 0;
    }
#endif

    return p;
}

uint32_t sample_poisson_pmf_f
(
    float lambda, 
    uint32_t a, 
    uint32_t b
)
{
    uint32_t s;
    uint32_t accept;
    float max_p, p, u;

    if (lambda <= 0 && ((lambda >= a) && (lambda <= b)))
    {
        return (int)lambda;
    }
    else if (lambda <= 0)
    {
        return (a - lambda) < (lambda - b) ? a : b;
    }

    max_p = poisson_pmf_f(lambda, floor(lambda + 0.5));
    accept = 0;

    while (!accept)
    {
        s = floor(sample_uniform_pdf_f((float)a, (float)b+1.0));
        if (s == b+1)
        {
            s = b;
        }
        p = poisson_pmf_f(lambda, s);
        u = sample_uniform_pdf_f(0, max_p);

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

    return s;
}

uint32_t sample_poisson_pmf_d
(
    double lambda, 
    uint32_t a, 
    uint32_t b
)
{
    uint32_t s;
    uint32_t accept;
    double max_p, p, u;

    if (lambda <= 0 && ((lambda >= a) && (lambda <= b)))
    {
        return (int)lambda;
    }
    else if (lambda <= 0)
    {
        return (a - lambda) < (lambda - b) ? a : b;
    }

    max_p = poisson_pmf_d(lambda, floor(lambda + 0.5));
    accept = 0;

    while (!accept)
    {
        s = floor(sample_uniform_pdf_d((double)a, (double)b+1.0));
        if (s == b+1)
        {
            s = b;
        }
        p = poisson_pmf_d(lambda, s);
        u = sample_uniform_pdf_d(0, max_p);

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

    return s;
}

float bernoulli_pmf_f(float p, uint32_t n)
{
    return n ? p : (1.0 - p);
}

double bernoulli_pmf_d(double p, uint32_t n)
{
    return n ? p : (1.0 - p);
}

uint8_t sample_bernoulli_pmf_f(float p)
{
    float u;

    assert(p >= 0.0 && p <= 1.0);

    u = sample_uniform_pdf_f(0.0, 1.0);

    if (u <= p)
    {
        return 1;
    }

    return 0;
}

uint8_t sample_bernoulli_pmf_d(double p)
{
    double u;

    assert(p >= 0.0 && p <= 1.0);

    u = sample_uniform_pdf_d(0.0, 1.0);

    if (u <= p)
    {
        return 1;
    }

    return 0;
}

float geometric_pmf_f(float p, uint32_t n)
{
    uint32_t i;
    float P;

    assert(p > 0 && p < 1);

    P = 1.0;

    for (i = 0; i < n; i++)
    {
        P *= (1.0 - p);
    }

    return P * p;
}

double geometric_pmf_d(double p, uint32_t n)
{
    uint32_t i;
    double P;

    assert(p > 0 && p < 1);

    P = 1.0;

    for (i = 0; i < n; i++)
    {
        P *= (1.0 - p);
    }

    return P * p;
}

float log_geometric_pmf_f(float p, uint32_t n)
{
    float q = 1 - p;

    assert(p > 0 && p < 1);
    if (p < JWSC_MIN_LOG_ARG)
    {
        p = JWSC_MIN_LOG_ARG;
    }
    else if (q < JWSC_MIN_LOG_ARG)
    {
        q = JWSC_MIN_LOG_ARG;
    }

    return log(p) + n*log(q);
}

double log_geometric_pmf_d(double p, uint32_t n)
{
    double q = 1 - p;

    assert(p > 0 && p < 1);
    if (p < JWSC_MIN_LOG_ARG)
    {
        p = JWSC_MIN_LOG_ARG;
    }
    else if (q < JWSC_MIN_LOG_ARG)
    {
        q = JWSC_MIN_LOG_ARG;
    }

    return log(p) + n*log(q);
}

uint32_t sample_geometric_pmf_f(float p)
{
    float u;

    assert(p > 0 && p < 1);

    u = sample_uniform_pdf_f(p, ALMOST_ONE_SINGLE);

    return (uint32_t)floor((log(1.0 - u) / log(1.0 - p)) - 1.0);
}

uint32_t sample_geometric_pmf_d(double p)
{
    double u;

    assert(p > 0 && p < 1);

    u = sample_uniform_pdf_d(p, ALMOST_ONE_DOUBLE);

    return (uint32_t)floor((log(1.0 - u) / log(1.0 - p)) - 1.0);
}