Distribution.cpp

Go to the documentation of this file.

/* -*- Mode: C++; -*- */
// VER: $Id$
// copyright (c) 2004 by Christos Dimitrakakis <dimitrak@idiap.ch>
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <learning/Distribution.h>

static const double PI = 3.14159265358979323846;

void setRandomSeed(unsigned int seed)
{
    srand(seed);
}
real urandom()
{
    real x;
    do {
        x = ((real) rand())/((real) (RAND_MAX));
    } while (x==1.0);
    return x;
}

real urandom(real min, real max)
{
    return min + ((max-min)*urandom());
}

real UniformDistribution::generate()
{
    return m + (urandom()-0.5)*s;
}

real UniformDistribution::pdf(real x)
{
    real dx=x-m;
    real sh=.5*s;
    if ((dx > -sh)&&(dx <= sh))
        return 1.0/s;
    return 0.0;
}



// Taken from numerical recipes in C
real NormalDistribution::generate()
{
    if(!cache) {
        normal_x = urandom();
        normal_y = urandom();
        normal_rho = sqrt(-2.0 * log(1.0 - normal_y));
        cache = true;
    } else {
        cache = false;
    }
    
    if (cache) {
        return normal_rho * cos(2.0 * PI * normal_x) * s + m;
    } else {
        return normal_rho * sin(2.0 * PI * normal_x) * s + m;   
    }
}

real NormalDistribution::pdf(real x)
{
    real d = (m-x)/s;
    return exp(-0.5 * d*d)/(sqrt(2.0 * PI) * s);
}

real LaplacianDistribution::generate()
{
    real x = urandom(-1.0, 1.0);
    real absx = fabs (x);
    real sgnx;
    if (x>0.0) {
        sgnx = 1.0;
    } else {
        sgnx = -1.0;
    }
    
    return m + sgnx * log(1.0 - absx) / l;

}

real LaplacianDistribution::pdf(real x)
{
    return 0.5*l * exp (-l*fabs(x-m));
}

real ExponentialDistribution::generate()
{
    real x = urandom();
    return - log (1.0 - x) / l;
}

real ExponentialDistribution::pdf(real x)
{
    real d = x - m;
    if (d>0.0) {
        return l * exp (-l*d);
    }
    return 0.0;
}


DiscreteDistribution::DiscreteDistribution() {
    p = NULL; n_outcomes=0;
}

DiscreteDistribution::DiscreteDistribution(int N) {
    p = NULL;
    n_outcomes = 0;
    p = (real*) malloc (sizeof(real) * N);
    n_outcomes = N;
    real invN = 1.0/((real) N);
    for (int i=0; i<N; i++) {
        p[i] = invN;
    }
}

DiscreteDistribution::~DiscreteDistribution() {
    free (p);
}

real DiscreteDistribution::generate()
{
    real d=urandom();
    real sum = 0.0;
    for (int i=0; i<n_outcomes; i++) {
        sum += p[i];
        if (d < sum) {
            return (real) i;
        }
    }
    return 0.0;
}

real DiscreteDistribution::pdf(real x)
{
    int i=(int) floor(x);
    if ((i>=0)&&(i<n_outcomes)) {
        return p[i];
    } 
    return 0.0;
}