Impala Documentation

template<class HistT>

void Impala::Core::Histogram::FitWeibull< HistT >::DoFit ( ) [inline, private]

Definition at line 55 of file FitWeibull.h.

Referenced by Impala::Core::Histogram::FitWeibull< HistT >::FitWeibull().

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         //    tol = 2.0*0.000001*gamma;
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     }

Here is the call graph for this function: