gaussian.h

/**************************************************************************
 * Copyright (c) 2002-2011 T. M. Murali                                   *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm 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 3 of the License, or      *
 * (at your option) any later version.                                    *
 *                                                                        *
 * Biorithm is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 * GNU General Public License for more details.                           *
 *                                                                        *
 * You should have received a copy of the GNU General Public License      *
 * along with Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#ifndef _GAUSSIAN_H
#define _GAUSSIAN_H

#include <math.h>

#include "constants.h"

// say the functions are INLINE!!!!!! otherwise, they get defined multiple times.
inline double gaussian_prob(double x, double mean, double sigma)
{
  if (0 == sigma)
    // this is a spike. gaussian's prob. is 1 iff x equals mean.
    return(x == mean);
  double std = (x - mean)/sigma;
  return(exp(-std*std/2)/(sigma*SQRT2PI));
}

inline double gaussian_log_prob(double x, double mean, double sigma)
{
  if (0 == sigma)
    // this gaussian is a spike. gaussian's prob. is 1 iff x equals mean.
    {
      double std = (x == mean);
      return(log(std));
    }
  double std = (x - mean)/sigma;
  return(-std*std/2 - log(sigma*SQRT2PI));
}

// return the prob. that the gaussian takes a value >= x. 
inline double gaussian_pval(double x, double mean, double sigma)
{
  return(0.5*erfc((x - mean)/(sigma*SQRT2)));
  
}


// compute the prob of a gaussian random variable falling inside an
// interval. this value is simply the integral of the gaussian between
// the two limits of the interval or the diff. between the pvalues of
// the two limits.
inline double prob_interval_gaussian(double start, double end, double mean,
                              double sigma)
  {
    return(gaussian_pval(start, mean, sigma) - 
           gaussian_pval(end, mean, sigma));
  }
  
#endif // _GAUSSIAN_H