NC::SpectrumInterpolatorFancy Class Reference

Given a set of reference histograms, interpolate between them. More...

#include <NCSpectrumInterpolator.h>

Inheritance diagram for NC::SpectrumInterpolatorFancy:
NC::ISpectrumInterpolator

List of all members.

Protected Member Functions

int intintpow (int a, int b) const
 Private utility function. Calculates a^b using integer arithmetic.
double doubleintpow (double a, int b) const
 Private utility function. Calculates a^b where b is an integer.
virtual TMatrixD FindCoeffs (const std::vector< std::vector< double > > &xs, const std::vector< double > &x) const
 Do the bulk of the work.

Detailed Description

Given a set of reference histograms, interpolate between them.

Consider a curve $y(x)$ fit to a set of $N$ points ${(y_n, x_n)}$ by a chi-squared statistic

\[ \chi^2(x)=\sum_n{(y(x_n)-y_n)^2\over\sigma(x, x_n)^2} \]

where the weight accorded to any point varies with $x$.

Let $y$ be a linear combination of some basis functions $L_j(x)$

\[ y(x)=\sum_ja_jL_j(x) \]

We want to solve for the $ \vec a $ that gives the best fit ie $ {{\rm d}\chi^2\over{\rm d}a_i}=0 $ for all $i$

\[ {1\over2}{{\rm d}\chi^2\over{\rm d}a_i}= \sum_n{(\sum_ka_kL_k(x_n)-y_n)L_i(x_n)\over\sigma_n^2}=0 \]

\[ \to \sum_n{\sum_ka_kL_k(x_n)L_i(x_n)\over\sigma_n^2}= \sum_n{L_i(x_n)y_n\over\sigma_n^2} \]

Naming the terms we have

\[ \sum_kM_{ik}a_k=\sum_nV_{in}y_n \]

and solving for $a$ gives

\[ a_k=(M^{-1})_{ki}V_{in}y_n \]

The final result $y(x) = \sum_k a_kL_k(x)$ can now be computed.

\[ y=(M^{-1})_{ki}V_{in}y_nL_k(x)=L^T(x)M^{-1}V\vec y \]

So the coefficients of all the $y_n$'s ($L^T(x)M^{-1}V$) can be precomputed as they are functions only of $x,x_n,y_n$.

The calculation is general and works for multidimensional $x$. Because the final combination of the $y_n$'s is linear and the coefficients can be precomputed we can operate directly over entire histograms.

The $L$ functions are chosen to be one-term polynomials, making up the full $y(x)$ polynomial. So, we need to compute the full set of possible combinations of powers to raise things to so we know what terms are possible.

\[ L_k(x) = \prod_d(x^d)^{P_{kd}} \]

Where $x$ is the input vector in D-dimensional space $\vec x = [x^1,...,x^D]$ and $P$ is the matrix of powers. above, pows in the code.

The calculation of the coefficients is achieved in ISpectrumInterpolator::FindCoeffs. $L$ is called basis, $M$ is mat, $V$ is vecPart, $\sigma^2$ is denom and $P$ is pows.

The function chosen for $\sigma$ is $\sigma(x,x_n)=|x-x_n|$

This has the feature that $1\over\sigma^2$ diverges whenever $x$ is any of the $x_n$'s

This is good in that it forces the fit function to go through that point by giving it infinite weight, but bad in that it causes divide-by-zeros. We detect when $x$ is very close to an $x_n$ and in that case construct the coefficient vector by hand so that it picks out the corresponding $y_n$.

Todo:
Try multiplying up by all the weights, so that, instead of the weight going infinite, it goes to zero on all the other terms. This should prevent divide-by-zero but we may still have problems with the matrix being non-invertible.

Definition at line 183 of file NCSpectrumInterpolator.h.


Member Function Documentation

double NC::SpectrumInterpolatorFancy::doubleintpow ( double  a,
int  b 
) const [protected]

Private utility function. Calculates a^b where b is an integer.

Definition at line 124 of file NCSpectrumInterpolator.cxx.

00125 {
00126   double ret = 1;
00127   while(b--) ret *= a;
00128   return ret;
00129 }

virtual TMatrixD NC::SpectrumInterpolatorFancy::FindCoeffs ( const std::vector< std::vector< double > > &  xs,
const std::vector< double > &  x 
) const [protected, virtual]

Do the bulk of the work.

Implements NC::ISpectrumInterpolator.

int NC::SpectrumInterpolatorFancy::intintpow ( int  a,
int  b 
) const [protected]

Private utility function. Calculates a^b using integer arithmetic.

Definition at line 116 of file NCSpectrumInterpolator.cxx.

00117 {
00118   int ret = 1;
00119   while(b--) ret *= a;
00120   return ret;
00121 }


The documentation for this class was generated from the following files:

Generated on 17 Jun 2018 for loon by  doxygen 1.6.1