Page MenuHomeHEPForge
Files F9501277

No OneTemporary
Actions

View Options

Index: trunk/npstat/stat/OrthoPolyGOFTest.icc
===================================================================
--- trunk/npstat/stat/OrthoPolyGOFTest.icc (revision 732)
+++ trunk/npstat/stat/OrthoPolyGOFTest.icc (revision 733)
@@ -1,144 +1,140 @@
#include <cassert>
#include <algorithm>
#include <stdexcept>
#include <cmath>
#include "npstat/nm/truncatedInverseSqrt.hh"
#include "npstat/nm/sumOfSquares.hh"
namespace npstat {
template <class Quadrature>
OrthoPolyGOFTest::OrthoPolyGOFTest(
const Quadrature& q,
const AbsClassicalOrthoPoly1D& poly,
const unsigned i_maxdeg,
const double* importanceWeights,
const unsigned nWeights)
: poly_(poly.clone()),
importanceWeights_(i_maxdeg+1U, 1.0),
- buf_((i_maxdeg+1U)*(i_maxdeg+2U)/2U),
+ buf_((i_maxdeg+1U)*(i_maxdeg+2U)/2U - 1U),
statBuf_(2U*i_maxdeg),
maxdeg_(i_maxdeg)
{
if (nWeights)
{
if (nWeights != i_maxdeg + 1U) throw std::invalid_argument(
"In npstat::OrthoPolyGOFTest constructor: "
"inconsistent number of importance weights");
assert(importanceWeights);
for (unsigned i=0; i<nWeights; ++i)
{
if (importanceWeights[i] <= 0.0) throw std::invalid_argument(
"In npstat::OrthoPolyGOFTest constructor: "
"all importance weights must be positive");
importanceWeights_[i] = importanceWeights[i];
}
}
calculateNormalizingMatrix(q, &invsqr_);
}
template <class Quadrature>
void OrthoPolyGOFTest::calculateNormalizingMatrix(
const Quadrature& q, Matrix<lapack_double>* result) const
{
assert(result);
const unsigned maxdegP1 = maxdeg_ + 1U;
- const unsigned nRows = maxdegP1*(maxdegP1+1U)/2U;
+ const unsigned nRows = maxdegP1*(maxdegP1+1U)/2U - 1U;
Matrix<lapack_double> covmat(nRows, nRows);
unsigned iRow = 0;
for (unsigned i=0; i<=maxdeg_; ++i)
{
const long double prod_i = importanceWeights_[i];
- for (unsigned k=i; k<=maxdeg_; ++k, ++iRow)
+ for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
{
assert(iRow < nRows);
const double k_del_ik = k == i ? 1.0 : 0.0;
const long double prod_k = prod_i*importanceWeights_[k];
unsigned iCol = 0;
for (unsigned j=0; j<=maxdeg_; ++j)
{
const long double prod_j = prod_k*importanceWeights_[j];
- for (unsigned m=j; m<=maxdeg_; ++m, ++iCol)
+ for (unsigned m=j ? j : 1U; m<=maxdeg_; ++m, ++iCol)
{
if (iCol >= iRow)
{
assert(iCol < nRows);
const long double prod_m = prod_j*importanceWeights_[m];
- double cov = 0.0;
- if (!((i == 0 && k == 0) || (j == 0 && m == 0)))
- {
- const double k_del_jm = j == m ? 1.0 : 0.0;
- cov = poly_->jointAverage(q, i, k, j, m) -
- k_del_ik*k_del_jm;
- if (i == j && k == m && cov < 0.0)
- cov = 0.0;
- }
+ const double k_del_jm = j == m ? 1.0 : 0.0;
+ double cov = poly_->jointAverage(q, i, k, j, m) -
+ k_del_ik*k_del_jm;
+ if (i == j && k == m && cov < 0.0)
+ cov = 0.0;
covmat[iRow][iCol] = cov*prod_m;
}
}
}
assert(iCol == nRows);
}
}
assert(iRow == nRows);
iRow = 0;
for (unsigned i=0; i<=maxdeg_; ++i)
{
- for (unsigned k=i; k<=maxdeg_; ++k, ++iRow)
+ for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
{
unsigned iCol = 0;
for (unsigned j=0; j<=maxdeg_; ++j)
{
- for (unsigned m=j; m<=maxdeg_; ++m, ++iCol)
+ for (unsigned m=j ? j : 1U; m<=maxdeg_; ++m, ++iCol)
if (iCol < iRow)
covmat[iRow][iCol] = covmat[iCol][iRow];
}
}
}
truncatedInverseSqrt(covmat, 2U*maxdeg_, 0.0, result);
assert(result->nRows() == 2U*maxdeg_);
}
template<typename Numeric>
void OrthoPolyGOFTest::normalizedDeviations(
const Numeric* data, const unsigned long lenData,
double* deviations, const unsigned lenDeviations) const
{
const unsigned twomaxdeg = 2U*maxdeg_;
assert(lenData);
assert(data);
assert(lenDeviations >= twomaxdeg);
assert(deviations);
const Matrix<double>& unnorm = poly_-> sampleProductAverages(
data, lenData, maxdeg_);
double* buf = const_cast<double*>(&buf_[0]);
unsigned iRow = 0;
for (unsigned i=0; i<=maxdeg_; ++i)
- for (unsigned k=i; k<=maxdeg_; ++k, ++iRow)
+ for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
buf[iRow] = unnorm[i][k]*importanceWeights_[i]*importanceWeights_[k];
assert(iRow == buf_.size());
invsqr_.timesVector(buf, iRow, deviations, twomaxdeg);
const double sqrlen = std::sqrt(static_cast<double>(lenData));
for (unsigned i=0; i<twomaxdeg; ++i)
deviations[i] *= sqrlen;
}
template<typename Numeric>
double OrthoPolyGOFTest::testStatistic(
const Numeric* data, const unsigned long lenData) const
{
double* statbuf = const_cast<double*>(&statBuf_[0]);
normalizedDeviations(data, lenData, statbuf, statBuf_.size());
return sumOfSquares(statbuf, statBuf_.size());
}
}
Index: trunk/npstat/stat/PolynomialDistro1D.cc
===================================================================
--- trunk/npstat/stat/PolynomialDistro1D.cc (revision 732)
+++ trunk/npstat/stat/PolynomialDistro1D.cc (revision 733)
@@ -1,143 +1,143 @@
#include <cassert>
#include <stdexcept>
#include <sstream>
#include <climits>
#include "npstat/nm/findRootUsingBisections.hh"
#include "npstat/stat/PolynomialDistro1D.hh"
#include "npstat/stat/distributionReadError.hh"
namespace npstat {
PolynomialDistro1D::PolynomialDistro1D(
const double location, const double scale,
const double* coeffs, const unsigned maxdeg,
const unsigned nCheck)
: AbsScalableDistribution1D(location, scale),
glq_(GaussLegendreQuadrature::minimalExactRule(maxdeg))
{
setupCoeffs(coeffs, maxdeg);
checkPositive(nCheck);
}
PolynomialDistro1D::PolynomialDistro1D(
const double location, const double scale,
const std::vector<double>& coeffs,
const unsigned nCheck)
: AbsScalableDistribution1D(location, scale),
glq_(GaussLegendreQuadrature::minimalExactRule(coeffs.size()))
{
setupCoeffs(coeffs.empty() ? (const double*)0 : &coeffs[0],
coeffs.size());
checkPositive(nCheck);
}
void PolynomialDistro1D::setupCoeffs(const double* coeffs,
const unsigned maxdeg)
{
if (maxdeg)
assert(coeffs);
allCoeffs_.resize(maxdeg + 1U);
allCoeffs_[0] = 1.0;
for (unsigned i=0; i<maxdeg; ++i)
allCoeffs_[i+1U] = coeffs[i];
}
void PolynomialDistro1D::checkPositive(const unsigned nCheck) const
{
if (nCheck < 2U) throw std::invalid_argument(
"In npstat::PolynomialDistro1D::checkPositive: "
"insufficient number of density check points");
const unsigned nM1 = nCheck - 1U;
const double step = 1.0/nM1;
for (unsigned i=0; i<nCheck; ++i)
{
double x = step*i;
if (i == nM1)
x = 1.0;
if (unscaledDensity(x) <= 0.0)
{
std::ostringstream os;
os.precision(17);
os << "In npstat::PolynomialDistro1D::checkPositive: "
- << "density is not positive at x = " << x;
+ << "density is not positive at x = " << x;
throw std::invalid_argument(os.str());
}
}
}
double PolynomialDistro1D::unscaledDensity(const double x) const
{
if (x < 0.0 || x > 1.0)
return 0.0;
else
return poly_.series(&allCoeffs_[0], allCoeffs_.size()-1, x);
}
double PolynomialDistro1D::unscaledCdf(const double x) const
{
if (x <= 0.0)
return 0.0;
else if (x >= 1.0)
return 1.0;
else
return glq_.integrate(UnscaledDensityFunctor1D(*this), 0.0, x);
}
double PolynomialDistro1D::unscaledQuantile(const double r1) const
{
if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
"In npstat::PolynomialDistro1D::unscaledQuantile: "
"cdf argument outside of [0, 1] interval");
if (r1 == 0.0)
return 0.0;
if (r1 == 1.0)
return 1.0;
double q;
const bool status = findRootUsingBisections(
UnscaledCdfFunctor1D(*this), r1, 0.0, 1.0, 4.0*DBL_EPSILON, &q);
assert(status);
return q;
}
double PolynomialDistro1D::unscaledExceedance(const double x) const
{
if (x <= 0.0)
return 1.0;
else if (x >= 1.0)
return 0.0;
else
return glq_.integrate(UnscaledDensityFunctor1D(*this), x, 1.0);
}
bool PolynomialDistro1D::isEqual(const AbsDistribution1D& other) const
{
const PolynomialDistro1D& r = static_cast<const PolynomialDistro1D&>(other);
return AbsScalableDistribution1D::isEqual(r) && allCoeffs_ == r.allCoeffs_;
}
bool PolynomialDistro1D::write(std::ostream& os) const
{
AbsScalableDistribution1D::write(os);
gs::write_pod_vector(os, allCoeffs_);
return !os.fail();
}
PolynomialDistro1D* PolynomialDistro1D::read(const gs::ClassId& id, std::istream& in)
{
static const gs::ClassId current(gs::ClassId::makeId<PolynomialDistro1D>());
current.ensureSameId(id);
double location, scale;
if (AbsScalableDistribution1D::read(in, &location, &scale))
{
std::vector<double> table;
gs::read_pod_vector(in, &table);
if (!in.fail() && table.size())
return new PolynomialDistro1D(location, scale,
&table[1], table.size()-1, 2U);
}
distributionReadError(in, classname());
return 0;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:13 PM (8 h, 5 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486489
Default Alt Text
(11 KB)

Event Timeline