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

---

FitWeibull.h

Go to the documentation of this file.

#ifndef Impala_Core_Histogram_FitWeibull_h
#define Impala_Core_Histogram_FitWeibull_h

#include <cmath>

namespace Impala
{
namespace Core
{
namespace Histogram
{


const double CXFITWEIBULL_TOL=0.001;


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

    /* input: histogram data with values 0..range */
    FitWeibull(HistT* hist, int bins, double range)
    {
        mHist = hist;
        mBins = bins;
        mDx = range/mBins;
        mPrecision = 0.0;
        DoFit();
    }

    double
    Beta()
    {
        return mBeta;
    }

    double
    Gamma()
    {
        return mGamma;
    }

    double
    A2()
    {
        return AndersonDarling();
    }

private:

    void
    DoFit()
    {
        double sumy(0), sumysq(0);
        double gammal, gammah;
        double xgm, gfm, tol;
        int i;

        mNorm = 0;
        double x;
        for (i=0, x = mDx*0.5; i<mBins; i++, x += mDx)
        {
            if (mHist[i] > mPrecision)
            {
                double logy = log(x);
                sumy += mHist[i]*logy;
                sumysq += mHist[i]*logy*logy;
                mNorm += mHist[i];
            }
        }

        sumy /= mNorm;
        sumysq /= mNorm;
        if( sumysq-sumy*sumy <= CXFITWEIBULL_TOL )
        {
            mGamma = 2.0;
            mBeta = 0.5*mDx;
            return;
        }

        mGamma = 1.28/sqrt((sumysq-sumy*sumy)*mBins/(mBins-1.0));
        xgm = exp(sumy);
        //    tol = 2.0*0.000001*gamma;
        tol = 2.0*CXFITWEIBULL_TOL*mGamma;
        gfm = GFunct(mGamma,xgm);

        int its = 0;

        if (gfm >= 0.0)
        {
            while (gfm > 0.0)
            {
                gammah = mGamma;
                mGamma /= 2.0;
                if (mGamma < tol)
                    break;
                gfm = GFunct(mGamma,xgm);
                if (its++ > 30)
                    break;
            }
            gammal = mGamma;
        }
        else
        {
            while (gfm < 0.0)
            {
                gammal = mGamma;
                mGamma *= 2.0;
                gfm = GFunct(mGamma,xgm);
                if (its++ > 30)
                    break;
            }
            gammah = mGamma;
        }

        while ((gammah-gammal) > tol)
        {
            mGamma = (gammal+gammah)/2.0;
            gfm = GFunct(mGamma,xgm);
            if (gfm >= 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,xgm);
    }

    double
    BetaEst(double gamma, double xgm) const
    { 
        double sum = 0;
        double ddx = mDx/xgm; // x[i]/xgm
        double x = 0.5*ddx;

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

        return xgm*pow(sum/mNorm, 1./gamma);
    }

    double
    GFunct(double gamma, double xgm) const
    { 
        double sum = 0;
        double beta = BetaEst(gamma, xgm);
        double ddx = mDx/beta;
        double x = 0.5*ddx;

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

        return sum/mNorm-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;
    }

    HistT* mHist;
    int    mBins;

    double mDx;
    double mNorm;

    double mBeta;
    double mGamma;
    double mPrecision;
};

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

#endif

---