C:/programs/etirm/src/DiscreteNormalDist.h

Go to the documentation of this file.

/*! \file DiscreteNormalDist.h
  
  \brief
  Function for computing a discrete approximation to a normal
  distribution.

  Estimation Toolkit for Item Response Models (ETIRM)
  http://www.smallwaters.com/software/cpp/etirm.html

  Author(s): 
  Werner Wothke, maintenance (http://www.smallwaters.com)
  Brad Hanson (http://www.b-a-h.com/)
  See the file LICENSE for information on usage and redistribution.

  Copyright (C) 2008, Werner Wothke
  Copyright (c) 2000-2001, Bradley A. Hanson
*/

#ifndef ETIRM_DISCRETENORMALDIST_H_
#define ETIRM_DISCRETENORMALDIST_H_

#ifdef ETIRM_NO_DIR_PREFIX
#include "etirmtypes.h"
#else
#include "etirm/etirmtypes.h"
#endif

#include <cmath> // for exp

// for compilers which do not put C library functions in std namespace
#ifdef BOOST_NO_STDC_NAMESPACE
namespace std {using ::exp;}
#endif

namespace etirm
{

/*!
  DiscreteNormalDist
  
  \brief
  Compute points and weights for a discrete distribution
  that approximates a normal distribution
  
  Arguments
  
  nQuad   Number of discrete points in distribution.
  minPoint  Value of smallest point in distribution (must be smaller than the mean).
  maxPoint  Value of largest point in distribution (must be larger than the mean).
  ipoints   Iterator over nQuad points of distribution.
  iweights  Iterator over nQuad probabilities corresponding to the points.
  mean    Mean of normal distribution.
  sd      Standard deviation of normal distribution.
*/
template<class PI, class WI>
void DiscreteNormalDist(int nQuad, double minPoint, double maxPoint, PI ipoints,
  WI iweights, double mean = 0.0, double sd = 1.0)
{
  int i;
  
  if ((maxPoint-minPoint) <= 0.0 || minPoint > mean || maxPoint < mean) 
    throw InvalidArgument("Invalid minimum and maximum points", "DiscreteNormalDist");
  
  double incr = (maxPoint-minPoint) / (nQuad-1.0);
  
  PI points = ipoints;
  WI weights = iweights;
  double theta = minPoint;
  double sum = 0.0;
  for (i=nQuad; i--; ++points, ++weights)
  {
    *points = (theta - mean) / sd;
    *weights = std::exp(-0.5 * theta * theta);
    sum += *weights;
    theta += incr;
  }
  
  /* Standardize weights */
  weights = iweights;
  for (i=nQuad; i--; ++weights)
  {
    *weights /= sum;
  }
  
}

} // namespace etirm

#endif // ETIRM_DISCRETENORMALDIST_H_

---