Page MenuHomeHEPForge
Files F10881403

No OneTemporary
Actions

View Options

Index: trunk/npstat/stat/OrthoPolyGOFTest.icc
===================================================================
--- trunk/npstat/stat/OrthoPolyGOFTest.icc (revision 733)
+++ trunk/npstat/stat/OrthoPolyGOFTest.icc (revision 734)
@@ -1,140 +1,141 @@
#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 - 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 - 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 ? 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 ? j : 1U; m<=maxdeg_; ++m, ++iCol)
{
if (iCol >= iRow)
{
assert(iCol < nRows);
const long double prod_m = prod_j*importanceWeights_[m];
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 ? i : 1U; k<=maxdeg_; ++k, ++iRow)
{
unsigned iCol = 0;
for (unsigned j=0; j<=maxdeg_; ++j)
{
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 ? i : 1U; k<=maxdeg_; ++k, ++iRow)
- buf[iRow] = unnorm[i][k]*importanceWeights_[i]*importanceWeights_[k];
+ buf[iRow] = (unnorm[i][k] - (i == k ? 1.0 : 0.0))*
+ 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());
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:20 AM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4982984
Default Alt Text
(5 KB)

Event Timeline