|| || || || || || || ||

---

FitWeibullMarginal.h

Go to the documentation of this file.

#ifndef Impala_Core_Histogram_FitWeibullMarginal_h
#define Impala_Core_Histogram_FitWeibullMarginal_h

#include <cmath>

namespace Impala
{
namespace Core
{
namespace Histogram
{


/* constants: pi, Euler's constant, and log(2) */
#ifndef M_PI
#define M_PI 3.1415926535897932384626433832795029
#endif
#ifndef MMGAMMA
#define MMGAMMA 0.5772156649015328606065120900824024
#endif
#ifndef M_LN2
#define M_LN2 0.6931471805599453094172321214581766
#endif


// HistT refers to the (elementary) data type used to store the histogram data,
// not to a "histogram" type
template<class HistT>
class FitWeibullMarginal
{
public:

    /* input: histogram of values minrange..maxrange */
    FitWeibullMarginal(HistT* hist, int bins, double minrange, double maxrange)
    {
        mHist = hist;
        mBins = bins;
        mMinval = minrange;
        mDx = (maxrange-minrange)/bins;
        mPrecision = 0.0;
        DoFit(hist, bins);
    }

    ~FitWeibullMarginal()
    {
        if (mHist)
            delete [] mHist;
    }

    double
    Mu()
    {
        return mMu;
    }

    double
    Beta()
    {
        return mBeta;
    }

    double
    Gamma()
    {
        return mGamma;
    }

    double
    A2()
    {
        return AndersonDarling();
    }

private:

    void
    DoFit(HistT* orgHist, int orgBins)
    {
        //std::cout << "org mBins = " << orgBins << std::endl;
        //mMu = MuEstNew();
        mMu = MuEst();
        if (mNorm == 0) // computed by MuEst, histogram is empty
        {
            mMu = 0;
            mBeta = 0;
            mHist = 0;
            return;
        }
        //std::cout << "mMu = " << mMu << std::endl;
        //std::cout << "mDx = " << mDx << std::endl;
        //std::cout << "bin for 0.0 = " << (0.0 - mMinval)/mDx << std::endl;
        //std::cout << "double bin0 = " << (mMu - mMinval)/mDx << std::endl;
        int bin0 = (int)((mMu - mMinval)/mDx);
        //std::cout << "bin0 = " << bin0 << std::endl;
        int i;

        // fill new histogram with mirrored data
        mBins = bin0>orgBins/2 ? bin0 : orgBins-bin0;
        //std::cout << "new mBins = " << mBins << std::endl;
        mHist = new HistT [mBins];

        for (i=0; i<mBins; i++)
            mHist[i] = 0;
        for (i=0; i<orgBins; i++)
        {
            int j = bin0 > i ? bin0-i-1 : i-bin0;
            //std::cout << "  org " << i << " to " << j << std::endl;
            mHist[j] += orgHist[i];
        }
        // remove empty tail
        for (i=mBins-1; i>0 && (mHist[i] <= mPrecision); i--)
            ;
        mBins = i+1;

        // initial estimate of upper and lower bound on gamma (between 0.3 and 50)
        double gammal, gammah;
        double score;
        int its = 0;

        mGamma = 2.0;
        score = GammaEst(mGamma);

        if (score > 0.0)
        {
            while(score > 0.0)
            {
                gammah = mGamma;
                mGamma /= 2.0;
                if (mGamma < 0.3)
                    break;
                score = GammaEst(mGamma);
                if (its++ > 30)
                    break;
            }
            gammal = mGamma;
        }
        else
        {
            while(score < 0.0)
            {
                gammal = mGamma;
                mGamma *= 2.0;
                if (mGamma > 50.0)
                    break;
                score = GammaEst(mGamma);
                if (its++ > 30)
                    break;
            }
            gammah = mGamma;
        }

        // approximate gamma by newton-raphson

        //    double tol = 2.0*TOL*2.0;
        //    while ((gammah-gammal) > tol) {

        while (its < 10)
        {
            mGamma = (gammal+gammah)/2.0;
            score = GammaEst(mGamma);
            if (score >= 0.0)
                gammah = mGamma;
            else
                gammal = mGamma;
            if (its++ > 30)
            {
                std::cerr << "Weibull MLE::no convergence..." << std::endl;
                break;
            }
        }

        mGamma = (gammal+gammah)/2.0;
        mBeta = BetaEst(mGamma);

    }

    double
    MuEst()
    {
        double sum = 0;
        double x = mMinval;

        mNorm = 0;
        for (int i=0; i<mBins; i++, x += mDx)
        {
            if (mHist[i] > mPrecision)
            {
                sum += x*mHist[i];
                mNorm += mHist[i];
            }
        }

        if (mNorm == 0)
            return 0;
        return sum/mNorm;
    }

    double
    MuEstNew()
    {
        mNorm = 0;
        for (int i=0 ; i<mBins ; i++)
        {
            mNorm += mHist[i];
            //std::cout << "sum " << i << " = " << mNorm << std::endl;
        }
        if (mNorm == 0)
            return 0;

        double sum = 0;
        int i=0;
        for (i=0 ; sum < 0.5*mNorm ; i++)
            sum += mHist[i];
        int medianbin = (i>0) ? i-1 : 0;
        return mMinval + (medianbin+0.5)*mDx;
    }

    double
    BetaEst(double gamma)
    {
        double sum = 0;
        double x = mDx;

        for (int i=0; i<mBins; i++, x+=mDx)
            if (mHist[i] > mPrecision)
                sum += pow(x, gamma)*mHist[i];

        if (mNorm == 0)
            return 0;
        return pow(sum/mNorm, 1./gamma);
    }

    double
    GammaEst(double gamma)
    { 
        double sum = 0;

        mBeta = BetaEst(gamma);

        double ddx = mDx/mBeta;
        double x = mDx/mBeta;

        for (int i=0; i<mBins; i++, x += ddx)
            if (mHist[i] > mPrecision)
                sum += pow(x, gamma)*(1.0-gamma*log(x))*mHist[i];

        return -(sum/mNorm+gamma-1.0+log(gamma)+DiGamma(1./gamma));
    }


    double
    AndersonDarling() const
    { 
        double A2 = 0;
        double x;
        int i, n=0;

        double Fi = 0.0;

        for (i=0, x = mDx; i<mBins; i++, x += mDx)
        {
            if (mHist[i] > mPrecision)
            {
                double xt = pow(x/mBeta, mGamma);
                double F = exp(-xt);
                double dF = (mGamma/mBeta)*pow(x/mBeta, mGamma-1.0)*F*mDx;
            
                Fi += mHist[i]/mNorm;
                F = 1.0-F;
                if ((F>0.0) && (F<1.0))
                {
                    A2 += ((Fi-F)*(Fi-F)/(F*(1-F)))*dF;
                    n++;
                }
            }
        }

        return n*A2;
    }

    double
    gammaFunc(double xx) const
    {
        double x,y,tmp,ser;
        static double cof[6] =
            {
                76.18009172947146, -86.50532032941677,
                24.01409824083091, -1.231739572450155,
                0.1208650973866179e-2, -0.5395239384953e-5
            };
        
        int j;
        y=x=xx;
        tmp=x+5.5;
        tmp -= (x+0.5)*log(tmp);
        ser=1.000000000190015;
        for (j=0;j<=5;j++)
            ser += cof[j]/++y;
        xx = -tmp+log(2.5066282746310005*ser/x);
        return exp(xx);
    }

    double
    DiGamma(double x) const
    {
        /* force into the interval 1..3 */
        if( x < 0.0 )
            return DiGamma(1.0-x)+M_PI/tan(M_PI*(1.0-x)) ;  /* refection formula */
        else if( x < 1.0 )
            return DiGamma(1.0+x)-1.0/x ;
        else if ( x == 1.0)
            return -MMGAMMA ;
        else if ( x == 2.0)
            return 1.0-MMGAMMA ;
        else if ( x == 3.0)
            return 1.5-MMGAMMA ;
        else if ( x > 3.0)
            /* duplication formula */
            return 0.5*(DiGamma(x/2.0)+DiGamma((x+1.0)/2.0))+M_LN2 ;
        else
        {
            static double Kncoe[] =
                { .30459198558715155634315638246624251,
                  .72037977439182833573548891941219706, -.12454959243861367729528855995001087,
                  .27769457331927827002810119567456810e-1, -.67762371439822456447373550186163070e-2,
                  .17238755142247705209823876688592170e-2, -.44817699064252933515310345718960928e-3,
                  .11793660000155572716272710617753373e-3, -.31253894280980134452125172274246963e-4,
                  .83173997012173283398932708991137488e-5, -.22191427643780045431149221890172210e-5,
                  .59302266729329346291029599913617915e-6, -.15863051191470655433559920279603632e-6,
                  .42459203983193603241777510648681429e-7, -.11369129616951114238848106591780146e-7,
                  .304502217295931698401459168423403510e-8, -.81568455080753152802915013641723686e-9,
                  .21852324749975455125936715817306383e-9, -.58546491441689515680751900276454407e-10,
                  .15686348450871204869813586459513648e-10, -.42029496273143231373796179302482033e-11,
                  .11261435719264907097227520956710754e-11, -.30174353636860279765375177200637590e-12,
                  .80850955256389526647406571868193768e-13, -.21663779809421233144009565199997351e-13,
                  .58047634271339391495076374966835526e-14, -.15553767189204733561108869588173845e-14,
                  .41676108598040807753707828039353330e-15, -.11167065064221317094734023242188463e-15 } ;

            double Tn_1 = 1.0 ; /* T_{n-1}(x), started at n=1 */
            double Tn = x-2.0 ; /* T_{n}(x) , started at n=1 */
            double resul = Kncoe[0] + Kncoe[1]*Tn ;

            x -= 2.0 ;

            for(int n = 2 ; n < sizeof(Kncoe)/sizeof(double) ;n++)
            {
                const double Tn1 = 2.0 * x * Tn - Tn_1 ;    /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */
                resul += Kncoe[n]*Tn1 ;
                Tn_1 = Tn ;
                Tn = Tn1 ;
            }
            return resul ;
        }
    }


    HistT* mHist;
    int    mBins;

    double mMinval;
    double mDx;
    double mNorm;

    double mMu;
    double mBeta;
    double mGamma;
    double mPrecision;
};

} // namespace Histogram
} // namespace Core
} // namespace Impala

#endif

---