JASTRONOMY::JAspera Struct Reference

Auxiliary data structure to fit signal strength using likelihood ratio. More...

#include <JAspera.hh>

Inheritance diagram for JASTRONOMY::JAspera:

Classes
struct	fit_type
	Result of fit. More...

Public Member Functions
void	put (const double s, const double b)
	Put signal and background to list of pre-computed N/S values.
void	put (const size_t n, const double s, const double b)
	Put signal and background to list of pre-computed N/S values.
double	getLikelihood (const double p) const
	Get likelihood for given signal strength.
double	getDerivative (const double p) const
	Get derivative of likelihood for given signal strength.
fit_type	operator() (const bool ns=false) const
	Fit signal strength.
double	getTestStatisticForUpperLimit (const double ps) const
	Get test statistic for given signal strength.
double	getSignal () const
	Get total signal strength.
void	setSignal (const double wS)
	Set signal strength.
void	addSignal (const double wS)
	Add signal strength.

Static Public Attributes
static constexpr double	EPSILON = 1.0e-3
	precision determination of signal strength

Protected Attributes
double	ws = 0.0
	total signal strength

Detailed Description

Auxiliary data structure to fit signal strength using likelihood ratio.

Definition at line 22 of file JAspera.hh.

Member Function Documentation

put() [1/2]

void JASTRONOMY::JAspera::put ( const double s, const double b )

inline

Put signal and background to list of pre-computed N/S values.

Parameters

s	signal
b	background

Definition at line 44 of file JAspera.hh.

    {
      push_back(b/s);
 
      ws += s;
    }

put() [2/2]

void JASTRONOMY::JAspera::put ( const size_t n, const double s, const double b )

inline

Put signal and background to list of pre-computed N/S values.

Parameters

n	data
s	signal
b	background

Definition at line 60 of file JAspera.hh.

    {
      for (size_t i = 0; i != n; ++i) {
        push_back(b/s);
      }
 
      ws += s;
    }

getLikelihood()

double JASTRONOMY::JAspera::getLikelihood ( const double p ) const

inline

Get likelihood for given signal strength.

Parameters

p	signal strength

Returns

likelihood

Definition at line 78 of file JAspera.hh.

    {
      double y = -p * ws;
 
      for (const double i : static_cast<const std::vector<double>&>(*this)) {
        y += log1p(p/i);
      }
 
      return y;
    }

getDerivative()

double JASTRONOMY::JAspera::getDerivative ( const double p ) const

inline

Get derivative of likelihood for given signal strength.

Parameters

p	signal strength

Returns

derivative of likelihood

Definition at line 96 of file JAspera.hh.

    {
      double y = -ws;
 
      for (const double i : static_cast<const std::vector<double>&>(*this)) {
        y += 1.0 / (p + i);
      }
 
      return y;
    }

operator()()

fit_type JASTRONOMY::JAspera::operator() ( const bool ns = false ) const

inline

Fit signal strength.

Parameters

ns	allow for negative signal

Returns

result

Definition at line 149 of file JAspera.hh.

    {
      using namespace std;
 
      if (this->empty()) {
 
        // nothing to be done
 
        return { 0.0, 0.0 };
 
      } else if (this->size() == 1 ) {
 
        // analytical solution
 
        const double x = 1.0/ws - (*this)[0];
 
        if (x > 0.0 || ns)
          return { getLikelihood(x), x };
        else
          return { 0.0, 0.0 };
      }
 
      double x1 = 0.0;
      double x2 = 0.0;
 
      double f1 = getDerivative(0.0);           // discriminator between positive and negative signal
      double f2 = 0.0;
 
      if        (f1 == 0.0) {
 
        return { 0.0, 0.0 };
 
      } else if (f1 >  0.0) {                   // positive signal
 
        x1 = 0.0;                               // lower limit corresponds to no signal
        x2 = (double) this->size() / ws;        // upper limit corresponds to no background (i.e. all N/S = 0)
 
        f2 = getDerivative(x2);
 
      } else if (ns) {                          // negative signal
 
        x2 = 0.0;                               // upper limit corresponds to no signal
        x1 = -(*this)[0];                       // lower limit corresponds to largest negated N/S ratio
 
        for (const double x : static_cast<const std::vector<double>&>(*this)) {
          if (-x > x1) {
            x1 = -x;
          }
        }
 
        x1 += 1.0 / ws;                         // offset
 
        f2 = f1;
        f1 = getDerivative(x1);
 
      } else {
 
        return { 0.0, 0.0 };
      }
 
#if   defined(BISECTION)
 
      // binary search
 
      while (x2 - x1 > EPSILON) {
 
        const double x = 0.5 * (x1 + x2);
 
        if (getDerivative(x) > 0.0)
          x1 = x;
        else
          x2 = x;
      }
 
      const double x = 0.5 * (x1 + x2);
 
      return { getLikelihood(x), x };
 
#elif defined(NEWTON_RAPHSON)
 
      // Newton-Raphson method
 
      double dx = x2 - x1;
      double x  = 0.5 * (x1 + x2);
 
      derivatives_type y = getDerivatives(x);
 
      struct {
        double dx;
      } old;
 
      old.dx = dx;
 
      for ( ; ; ) {
 
        if ((y.fp + (x2 - x)*y.fpp)*(y.fp + (x1 - x)*y.fpp) > 0.0 || fabs(2.0*y.fp) > fabs(old.dx*y.fpp)) {
 
          old.dx = dx;
 
          dx  = 0.5 * (x2 - x1);
          x   = x1 + dx;
 
        } else {
 
          old.dx = dx;
 
          dx  = y.fp/y.fpp;
          x  -= dx;
        }
 
        if (fabs(dx) <= EPSILON) {
          break;
        }
 
        y = getDerivatives(x);
 
        if (y.fp < 0.0)
          x2 = x;
        else
          x1 = x;
      }
 
      return { getLikelihood(x), x };
 
#else
 
      // Ridder's method
 
      while (x2 - x1 > EPSILON) {
 
        const double xm = 0.5 * (x1 + x2);
        const double fm = getDerivative(xm);
 
        const double s  = sqrt(fm*fm - f1*f2);
 
        if (s == 0.0) {
          break;
        }
 
        const double xn = xm + (xm - x1) * fm/s;
        const double fn = getDerivative(xn);
 
        if (fn == 0.0) {
          return { getLikelihood(xn), xn };
        }
 
        if (signbit(fn) != signbit(fm)) {
 
          x1 = xm;
          f1 = fm;
          x2 = xn;
          f2 = fn;
 
        } else {
 
          if (signbit(fn)) {
 
            x2 = xn;
            f2 = fn;
 
          } else {
 
            x1 = xn;
            f1 = fn;
          }
        }
      }
 
      const double x = 0.5 * (x1 + x2);
 
      return { getLikelihood(x), x };
#endif
    }

getTestStatisticForUpperLimit()

double JASTRONOMY::JAspera::getTestStatisticForUpperLimit ( const double ps ) const

inline

Get test statistic for given signal strength.

See formula 16 in this reference.

Parameters

ps	signal strength

Returns

test statistic

Definition at line 331 of file JAspera.hh.

    {
      const fit_type result = (*this)(true);
 
      if      (result.signal <= 0.0)
        return 0.0                -  this->getLikelihood(ps);
      else if (result.signal <= ps)
        return result.likelihood  -  this->getLikelihood(ps);
      else
        return 0.0;
    }

getSignal()

double JASTRONOMY::JAspera::getSignal ( ) const

inline

Get total signal strength.

Returns

signal strength

Definition at line 349 of file JAspera.hh.

    {
      return ws;
    }

setSignal()

void JASTRONOMY::JAspera::setSignal ( const double wS )

inline

Set signal strength.

Parameters

wS	signal strength

Definition at line 360 of file JAspera.hh.

    {
      ws = wS;
    }

addSignal()

void JASTRONOMY::JAspera::addSignal ( const double wS )

inline

Add signal strength.

Parameters

wS	signal strength

Definition at line 371 of file JAspera.hh.

    {
      ws += wS;
    }

Member Data Documentation

EPSILON

double JASTRONOMY::JAspera::EPSILON = 1.0e-3

staticconstexpr

precision determination of signal strength

Definition at line 26 of file JAspera.hh.

ws

double JASTRONOMY::JAspera::ws = 0.0

protected

total signal strength

Definition at line 378 of file JAspera.hh.

---

The documentation for this struct was generated from the following file:

• JAspera.hh