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

---

FitWeibullMarginalOld.h

Go to the documentation of this file.

#ifndef Impala_Core_Histogram_FitWeibullMarginalOld_h
#define Impala_Core_Histogram_FitWeibullMarginalOld_h

//#include <algorithm>
#include <cmath>

namespace Impala
{
namespace Core
{
namespace Histogram
{


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

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

    ~FitWeibullMarginalOld()
    {
        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)
    {
        // find mu for folded histogram (2 bins together)
        double x;
        double mode = 0.0;
        int bin0 = 0;
        int i;
        mMu = 0.0;
        for (i=0, x = mMinval+0.5*mDx; i<orgBins; i++, x += mDx)
            mMu += orgHist[i]*x;
        bin0 = (int)((mMu-mMinval)/mDx);

        // fill new histogram with mirrored data
        mBins = bin0>orgBins/2 ? bin0 : orgBins-bin0;
        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;
            mHist[j] += orgHist[i];
        }

        double sumy(0), sumysq(0);
        double gammal, gammah;
        double xgm, gfm, tol;

        mNorm = 0;
        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;
        mGamma = 1.28/sqrt((sumysq-sumy*sumy)*mBins/(mBins-1.0));
        xgm = exp(sumy);
        //    tol = 2.0*0.000001*gamma;
        tol = 2.0*0.001*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;
        }
        its = 0;
        while (gfm < 0.0)
        {
            gammal = mGamma;
            mGamma *= 2.0;
            gfm = GFunct(mGamma,xgm);
            if (its++ > 30)
                break;
        }
        gammah = mGamma;

        its = 0;
        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)
                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*beta)*(pow(x, gamma)-1.0)*mHist[i];

        return sum/mNorm-1./gamma;
    }

    double
    AndersonDarling() const
    { 
        double A2 = 0;
        double x;
        int i, bins = 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 += (pow(Fi-F,2.0)/(F*(1-F)))*dF;
                bins++;
            }
        }

        return bins*A2;
    }

    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

---