00001 #ifndef Impala_Core_Histogram_FitWeibull_h
00002 #define Impala_Core_Histogram_FitWeibull_h
00003
00004 #include <cmath>
00005
00006 namespace Impala
00007 {
00008 namespace Core
00009 {
00010 namespace Histogram
00011 {
00012
00013
00014 const double CXFITWEIBULL_TOL=0.001;
00015
00016
00017
00018
00019 template<class HistT>
00020 class FitWeibull
00021 {
00022 public:
00023
00024
00025 FitWeibull(HistT* hist, int bins, double range)
00026 {
00027 mHist = hist;
00028 mBins = bins;
00029 mDx = range/mBins;
00030 mPrecision = 0.0;
00031 DoFit();
00032 }
00033
00034 double
00035 Beta()
00036 {
00037 return mBeta;
00038 }
00039
00040 double
00041 Gamma()
00042 {
00043 return mGamma;
00044 }
00045
00046 double
00047 A2()
00048 {
00049 return AndersonDarling();
00050 }
00051
00052 private:
00053
00054 void
00055 DoFit()
00056 {
00057 double sumy(0), sumysq(0);
00058 double gammal, gammah;
00059 double xgm, gfm, tol;
00060 int i;
00061
00062 mNorm = 0;
00063 double x;
00064 for (i=0, x = mDx*0.5; i<mBins; i++, x += mDx)
00065 {
00066 if (mHist[i] > mPrecision)
00067 {
00068 double logy = log(x);
00069 sumy += mHist[i]*logy;
00070 sumysq += mHist[i]*logy*logy;
00071 mNorm += mHist[i];
00072 }
00073 }
00074
00075 sumy /= mNorm;
00076 sumysq /= mNorm;
00077 if( sumysq-sumy*sumy <= CXFITWEIBULL_TOL )
00078 {
00079 mGamma = 2.0;
00080 mBeta = 0.5*mDx;
00081 return;
00082 }
00083
00084 mGamma = 1.28/sqrt((sumysq-sumy*sumy)*mBins/(mBins-1.0));
00085 xgm = exp(sumy);
00086
00087 tol = 2.0*CXFITWEIBULL_TOL*mGamma;
00088 gfm = GFunct(mGamma,xgm);
00089
00090 int its = 0;
00091
00092 if (gfm >= 0.0)
00093 {
00094 while (gfm > 0.0)
00095 {
00096 gammah = mGamma;
00097 mGamma /= 2.0;
00098 if (mGamma < tol)
00099 break;
00100 gfm = GFunct(mGamma,xgm);
00101 if (its++ > 30)
00102 break;
00103 }
00104 gammal = mGamma;
00105 }
00106 else
00107 {
00108 while (gfm < 0.0)
00109 {
00110 gammal = mGamma;
00111 mGamma *= 2.0;
00112 gfm = GFunct(mGamma,xgm);
00113 if (its++ > 30)
00114 break;
00115 }
00116 gammah = mGamma;
00117 }
00118
00119 while ((gammah-gammal) > tol)
00120 {
00121 mGamma = (gammal+gammah)/2.0;
00122 gfm = GFunct(mGamma,xgm);
00123 if (gfm >= 0.0)
00124 gammah = mGamma;
00125 else
00126 gammal = mGamma;
00127 if (its++ > 30) {
00128 std::cerr << "Weibull MLE::no convergence..." << std::endl;
00129 break;
00130 }
00131 }
00132
00133 mGamma = (gammal+gammah)/2.0;
00134 mBeta = BetaEst(mGamma,xgm);
00135 }
00136
00137 double
00138 BetaEst(double gamma, double xgm) const
00139 {
00140 double sum = 0;
00141 double ddx = mDx/xgm;
00142 double x = 0.5*ddx;
00143
00144 for (int i=0; i<mBins; i++, x += ddx)
00145 if (mHist[i] > mPrecision)
00146 sum += pow(x, gamma)*mHist[i];
00147
00148 return xgm*pow(sum/mNorm, 1./gamma);
00149 }
00150
00151 double
00152 GFunct(double gamma, double xgm) const
00153 {
00154 double sum = 0;
00155 double beta = BetaEst(gamma, xgm);
00156 double ddx = mDx/beta;
00157 double x = 0.5*ddx;
00158
00159 for (int i=0; i<mBins; i++, x += ddx)
00160 if (mHist[i] > mPrecision)
00161 sum += log(x)*(pow(x, gamma)-1.0)*mHist[i];
00162
00163 return sum/mNorm-1./gamma;
00164 }
00165
00166 double
00167 AndersonDarling() const
00168 {
00169 double A2 = 0;
00170 double x;
00171 int i, n=0;
00172
00173 double Fi = 0.0;
00174
00175 for (i=0, x = mDx; i<mBins; i++, x += mDx)
00176 {
00177 if (mHist[i] > mPrecision)
00178 {
00179 double xt = pow(x/mBeta, mGamma);
00180 double F = exp(-xt);
00181 double dF = (mGamma/mBeta)*pow(x/mBeta, mGamma-1.0)*F*mDx;
00182
00183 Fi += mHist[i]/mNorm;
00184 F = 1.0-F;
00185 if ((F>0.0) && (F<1.0)) {
00186 A2 += ( (Fi-F)*(Fi-F)/(F*(1-F)))*dF;
00187 n++;
00188 }
00189 }
00190 }
00191
00192 return n*A2;
00193 }
00194
00195 HistT* mHist;
00196 int mBins;
00197
00198 double mDx;
00199 double mNorm;
00200
00201 double mBeta;
00202 double mGamma;
00203 double mPrecision;
00204 };
00205
00206 }
00207 }
00208 }
00209
00210 #endif