Impala Documentation

template<class HistT>
void Impala::Core::Histogram::FitWeibullMarginal< HistT >::DoFit	(	HistT *	orgHist,
		int	orgBins
	)			`[inline, private]`
Definition at line 77 of file FitWeibullMarginal.h.
Referenced by Impala::Core::Histogram::FitWeibullMarginal< HistT >::FitWeibullMarginal().
00078     {
00079         //std::cout << "org mBins = " << orgBins << std::endl;
00080         //mMu = MuEstNew();
00081         mMu = MuEst();
00082         if (mNorm == 0) // computed by MuEst, histogram is empty
00083         {
00084             mMu = 0;
00085             mBeta = 0;
00086             mHist = 0;
00087             return;
00088         }
00089         //std::cout << "mMu = " << mMu << std::endl;
00090         //std::cout << "mDx = " << mDx << std::endl;
00091         //std::cout << "bin for 0.0 = " << (0.0 - mMinval)/mDx << std::endl;
00092         //std::cout << "double bin0 = " << (mMu - mMinval)/mDx << std::endl;
00093         int bin0 = (int)((mMu - mMinval)/mDx);
00094         //std::cout << "bin0 = " << bin0 << std::endl;
00095         int i;
00096 
00097         // fill new histogram with mirrored data
00098         mBins = bin0>orgBins/2 ? bin0 : orgBins-bin0;
00099         //std::cout << "new mBins = " << mBins << std::endl;
00100         mHist = new HistT [mBins];
00101 
00102         for (i=0; i<mBins; i++)
00103             mHist[i] = 0;
00104         for (i=0; i<orgBins; i++)
00105         {
00106             int j = bin0 > i ? bin0-i-1 : i-bin0;
00107             //std::cout << "  org " << i << " to " << j << std::endl;
00108             mHist[j] += orgHist[i];
00109         }
00110         // remove empty tail
00111         for (i=mBins-1; i>0 && (mHist[i] <= mPrecision); i--)
00112             ;
00113         mBins = i+1;
00114 
00115         // initial estimate of upper and lower bound on gamma (between 0.3 and 50)
00116         double gammal, gammah;
00117         double score;
00118         int its = 0;
00119 
00120         mGamma = 2.0;
00121         score = GammaEst(mGamma);
00122 
00123         if (score > 0.0)
00124         {
00125             while(score > 0.0)
00126             {
00127                 gammah = mGamma;
00128                 mGamma /= 2.0;
00129                 if (mGamma < 0.3)
00130                     break;
00131                 score = GammaEst(mGamma);
00132                 if (its++ > 30)
00133                     break;
00134             }
00135             gammal = mGamma;
00136         }
00137         else
00138         {
00139             while(score < 0.0)
00140             {
00141                 gammal = mGamma;
00142                 mGamma *= 2.0;
00143                 if (mGamma > 50.0)
00144                     break;
00145                 score = GammaEst(mGamma);
00146                 if (its++ > 30)
00147                     break;
00148             }
00149             gammah = mGamma;
00150         }
00151 
00152         // approximate gamma by newton-raphson
00153 
00154         //    double tol = 2.0*TOL*2.0;
00155         //    while ((gammah-gammal) > tol) {
00156 
00157         while (its < 10)
00158         {
00159             mGamma = (gammal+gammah)/2.0;
00160             score = GammaEst(mGamma);
00161             if (score >= 0.0)
00162                 gammah = mGamma;
00163             else
00164                 gammal = mGamma;
00165             if (its++ > 30)
00166             {
00167                 std::cerr << "Weibull MLE::no convergence..." << std::endl;
00168                 break;
00169             }
00170         }
00171 
00172         mGamma = (gammal+gammah)/2.0;
00173         mBeta = BetaEst(mGamma);
00174 
00175     }
Here is the call graph for this function: