Commit 452150c9 authored by Lynn Garren's avatar Lynn Garren
Browse files

first implementation of RandSkewNormal

parent 5863f771
2011-05-27 Mark Fischler <mf@fnal.gov> and Lynn Garren <garren@fnal.gov>
* Random/RandSkewNormal.h, Random/RandSkewNormal.icc,
src/RandSkewNormal.cc: An implementation of the Azzalini
skew-normal distribution, as requested in bug #75534
==============================
11.11.10 Release CLHEP-2.1.0.1
==============================
......
......@@ -56,6 +56,8 @@ pkginclude_HEADERS = \
RandPoissonQ.icc \
RandPoissonT.h \
RandPoissonT.icc \
RandSkewNormal.h \
RandSkewNormal.icc \
RandStudentT.h \
RandStudentT.icc \
RanecuEngine.h \
......
// $Id: RandSkewNormal.h,v 1.1 2011/05/27 20:36:28 garren Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
// HEP Random
// --- RandSkewNormal ---
// class header file
// -----------------------------------------------------------------------
// RandSkewNormal ---
// returns a skew-normal distribution with shape parameter k
// To get a distribution with scale parameter b and location m:
// r = m + b * RandSkewNormal.fire(k);
// http://azzalini.stat.unipd.it/SN/
// algorithm from K. McFarlane, June 2010.
// =======================================================================
// M Fischler and L Garren - Created: 26 May 2011
// =======================================================================
#ifndef RandSkewNormal_h
#define RandSkewNormal_h 1
#include "CLHEP/Random/defs.h"
#include "CLHEP/Random/Random.h"
#include "CLHEP/Utility/memory.h"
namespace CLHEP {
/**
* @author <mf@fnal.gov>
* @ingroup random
*/
class RandSkewNormal : public HepRandom {
public:
inline RandSkewNormal ( HepRandomEngine& anEngine, double shape=0. );
inline RandSkewNormal ( HepRandomEngine* anEngine, double shape=0. );
// These constructors should be used to instantiate a RandSkewNormal
// distribution object defining a local engine for it.
// The static generator will be skipped using the non-static methods
// defined below.
// If the engine is passed by pointer the corresponding engine object
// will be deleted by the RandSkewNormal destructor.
// If the engine is passed by reference the corresponding engine object
// will not be deleted by the RandSkewNormal destructor.
virtual ~RandSkewNormal();
// Destructor
// Static methods to shoot random values using the static generator
static double shoot();
static double shoot( double shape );
static void shootArray ( const int size, double* vect,
double shape=0. );
// Static methods to shoot random values using a given engine
// by-passing the static generator.
static double shoot( HepRandomEngine* anEngine );
static double shoot( HepRandomEngine* anEngine, double shape );
static void shootArray ( HepRandomEngine* anEngine, const int size,
double* vect, double shape=0. );
// Methods using the localEngine to shoot random values, by-passing
// the static generator.
double fire();
double fire( double shape );
void fireArray ( const int size, double* vect );
void fireArray ( const int size, double* vect, double shape );
double operator()();
double operator()( double shape );
// Save and restore to/from streams
std::ostream & put ( std::ostream & os ) const;
std::istream & get ( std::istream & is );
std::string name() const;
HepRandomEngine & engine();
static std::string distributionName() {return "RandSkewNormal";}
// Provides the name of this distribution class
protected:
static double gaussianSkewNormal ( HepRandomEngine *e, double k);
double getShapeParameter() { return shapeParameter; }
inline HepRandomEngine* getLocalEngine();
private:
shared_ptr<HepRandomEngine> localEngine;
double shapeParameter;
};
} // namespace CLHEP
#ifdef ENABLE_BACKWARDS_COMPATIBILITY
// backwards compatibility will be enabled ONLY in CLHEP 1.9
using namespace CLHEP;
#endif
#include "CLHEP/Random/RandSkewNormal.icc"
#endif
// $Id: RandSkewNormal.icc,v 1.1 2011/05/27 20:36:28 garren Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
// HEP Random
// --- RandSkewNormal ---
// inlined functions implementation file
// -----------------------------------------------------------------------
// =======================================================================
// M Fischler and L Garren - Created: 26 May 2011
// =======================================================================
namespace CLHEP {
RandSkewNormal::RandSkewNormal(HepRandomEngine & anEngine, double shape )
: HepRandom(), localEngine(&anEngine, do_nothing_deleter()),
shapeParameter(shape) {}
RandSkewNormal::RandSkewNormal(HepRandomEngine * anEngine, double shape )
: HepRandom(), localEngine(anEngine), shapeParameter(shape) {}
HepRandomEngine * RandSkewNormal::getLocalEngine()
{
return localEngine.get();
}
} // namespace CLHEP
......@@ -38,6 +38,7 @@ libCLHEP_Random_@VERSION@_a_SOURCES = \
RandPoisson.cc \
RandPoissonQ.cc \
RandPoissonT.cc \
RandSkewNormal.cc \
RandStudentT.cc \
RanecuEngine.cc \
Ranlux64Engine.cc \
......
// $Id: RandSkewNormal.cc,v 1.1 2011/05/27 20:36:57 garren Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
// HEP Random
// --- RandSkewNormal ---
// class implementation file
// -----------------------------------------------------------------------
// =======================================================================
// M Fischler and L Garren - Created: 26 May 2011
// =======================================================================
#include "CLHEP/Random/defs.h"
#include "CLHEP/Random/RandSkewNormal.h"
#include "CLHEP/Random/RandGaussT.h"
#include "CLHEP/Random/DoubConv.hh"
namespace CLHEP {
std::string RandSkewNormal::name() const {return "RandSkewNormal";}
HepRandomEngine & RandSkewNormal::engine() {return *localEngine;}
RandSkewNormal::~RandSkewNormal() {}
double RandSkewNormal::operator()()
{
return fire( shapeParameter );
}
double RandSkewNormal::operator()( double shape )
{
return fire( shape );
}
//-------------
double RandSkewNormal::shoot(HepRandomEngine* anEngine)
{
// really dumb use of RandSkewNormal
double k = 1;
return gaussianSkewNormal( anEngine, k );
}
double RandSkewNormal::shoot(HepRandomEngine* anEngine, double shape)
{
return gaussianSkewNormal( anEngine, shape );
}
double RandSkewNormal::shoot()
{
// really dumb use of RandSkewNormal
HepRandomEngine* anEngine = HepRandom::getTheEngine();
double k = 1;
return gaussianSkewNormal( anEngine, k );
}
double RandSkewNormal::shoot(double shape)
{
HepRandomEngine* anEngine = HepRandom::getTheEngine();
return gaussianSkewNormal( anEngine, shape );
}
void RandSkewNormal::shootArray( const int size, double* vect,
double shape )
{
for( double* v = vect; v != vect+size; ++v )
*v = shoot(shape);
}
void RandSkewNormal::shootArray(HepRandomEngine* anEngine, const int size,
double* vect, double shape )
{
for( double* v = vect; v != vect+size; ++v )
*v = shoot(anEngine, shape);
}
//-------------
double RandSkewNormal::fire() {
return gaussianSkewNormal( getLocalEngine(), getShapeParameter() );
}
double RandSkewNormal::fire(double shape) {
return gaussianSkewNormal( getLocalEngine(), shape );
}
void RandSkewNormal::fireArray( const int size, double* vect)
{
for( double* v = vect; v != vect+size; ++v )
*v = fire( shapeParameter );
}
void RandSkewNormal::fireArray( const int size, double* vect,
double shape )
{
for( double* v = vect; v != vect+size; ++v )
*v = fire( shape );
}
//-------------
double RandSkewNormal::gaussianSkewNormal ( HepRandomEngine* e, double k)
{
// RandSkewNormal is an implementation of Azzalini's SN generator
// http://azzalini.stat.unipd.it/SN/
// K. McFarlane, June 2010.
// For location parameter m = 0., scale = 1.
// To get a distribution with scale parameter b and location m:
// r = m + b * RandSkewNormal.fire(k);
double u[2] = {0.};
RandGaussT::shootArray(e, 2, u, 0, 1);
double delta = k/sqrt(1. + k*k);
double u1 = delta*u[0] + sqrt(1 - delta*delta)*u[1];
double r = u[0] >= 0 ? u1 : -u1;
return r;
}
//-------------
std::ostream & RandSkewNormal::put ( std::ostream & os ) const {
int pr=os.precision(20);
std::vector<unsigned long> t(2);
os << " " << name() << "\n";
os << "Uvec" << "\n";
t = DoubConv::dto2longs(shapeParameter);
os << shapeParameter << " " << t[0] << " " << t[1] << "\n";
os.precision(pr);
return os;
}
std::istream & RandSkewNormal::get ( std::istream & is ) {
std::string inName;
is >> inName;
if (inName != name()) {
is.clear(std::ios::badbit | is.rdstate());
std::cerr << "Mismatch when expecting to read state of a "
<< name() << " distribution\n"
<< "Name found was " << inName
<< "\nistream is left in the badbit state\n";
return is;
}
if (possibleKeywordInput(is, "Uvec", shapeParameter)) {
std::vector<unsigned long> t(2);
is >> shapeParameter >> t[0] >> t[1]; shapeParameter = DoubConv::longs2double(t);
return is;
}
// is >> shapeParameter encompassed by possibleKeywordInput
return is;
}
} // namespace CLHEP
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment