See

• The definition simdistribution.h
• The description Random distributions

<project-file type=“source”/> <content> /* * Copyright 2006 Junling Ma * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * This library 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 * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA */

#include “simdistribution.h” #include “genrand2.h” #include <algorithm> #include “simutil.h”

using namespace std;

#define MEAN_DRAW 10000 float SimDistribution::mean() const {

  float m = 0;
  for (int i = 0; i < MEAN_DRAW; i++)
      m += draw();
  return m / MEAN_DRAW;

}

float expRandom(float rate) {

if (rate==0)
   return SIM_INFINITY;
else if (isInfinity(rate))
   return 0;
else return -log(1-genrand2d())/rate;

}

SimExpDistribution::SimExpDistribution(float rate) : SimDistribution() {

  Rate=rate;

}

SimExpDistribution::~SimExpDistribution() { }

float SimExpDistribution::draw() const {

  return expRandom(Rate);

}

SimNormalDistribution::SimNormalDistribution(float mean, float std) : SimDistribution() {

  Mean=mean;
  Std=std;

}

SimNormalDistribution::~SimNormalDistribution() { }

float SimNormalDistribution::draw() const {

  // Using Box-Muller method 
  float X, Y, U,V,S;
  do {
      U = 2.0 * genrand2d() - 1;
      V = 2.0 * genrand2d() - 1;
      S = U * U + V * V;
  } while (S == 0.0 || S > 1.0);
  S = sqrt(-2 * log(S) / S);
  return Mean + Std * S * U;

}

SimGammaDistribution::SimGammaDistribution(float gamma, float lambda) {

  Gamma = gamma;
  Lambda = lambda;

}

SimGammaDistribution::~SimGammaDistribution() { }

float SimGammaDistribution::draw() const {

  return gammaRandom(Gamma, Lambda);

}

SimUniformDistribution::SimUniformDistribution(float u1, float u2) {

  Left = u1;
  Length = u2 - u1;

}

SimUniformDistribution::~SimUniformDistribution() { }

float SimUniformDistribution::draw() const {

  return ( Left + Length * genrand2d() );

}

SimCutOff::SimCutOff(SimDistribution* dist, float left, float right) {

  Distribution = dist;
  Left = left;
  Right = right;

}

SimCutOff::~SimCutOff() {

  if ( Distribution != NULL )
      delete Distribution;

}

float SimCutOff::draw() const {

  if ( Distribution == NULL )
      return Left;
  float r = Distribution->draw();
  while ( r < Left || r > Right )
      r = Distribution->draw();
  return r;

}

SimDiscreteDistribution::SimDiscreteDistribution() { }

SimDiscreteDistribution::~SimDiscreteDistribution() { }

SimDiscreteDistribution::Point::Point(float value, float freq) {

  Value = value;
  Frequency = freq;
  Sum = 0;

}

void SimDiscreteDistribution::set(float value, float freq) {

  std::vector<Point>::iterator i = lower_bound( Points.begin(), Points.end(), value );
  if ( i != Points.end() && *i == value )
      i->frequency() = freq ;
  else
      i = Points.insert(i, Point(value, freq));
  if ( i == Points.begin() )
      i->sum() = freq;
  else
      i->sum() = (i-1)->sum() + freq;
  i++;
  while ( i != Points.end() ) {
      i->sum() += freq;
      i++;
  }

}

float SimDiscreteDistribution::draw() const {

  int Size = Points.size();
  switch (Size) {
  case 0:
      return SIM_INFINITY; 
  case 1:
      return point(0);
  default:
      float r = genrand2d() * sum(Size-1);
      if ( r < sum(0) )
          return point(0);
      if ( r > sum(Size-2) )
          return point(Size-1);
      // find the interval with binary search
      int left = 0; 
      int right = Size - 1; 
      int mid = ( left + right ) / 2;
      float smid = sum(mid);
      while ( ( right - left ) > 2 ) {
          if ( r > smid )
              left = mid;
          else
              right = mid;
          mid = ( left + right ) / 2;
          smid = sum(mid);
      }
      if ( r > smid )
          return point(right);
      if ( r > sum(left) )
          return point(mid);
      return point(left);
  }

}

float SimDiscreteDistribution::mean() const {

  int n = Points.size();
  float m = 0;
  for (int i = 0; i < n; i++)
      m += Points[i].value() * Points[i].frequency();
  return m / sum(n);

} </content> <use name=“simdistribution.h”/> <use name=“genrand2.h”/>