00001
00002
00004 #ifdef WIN32
00005 #pragma warning(disable:4786)
00006 #endif
00007
00008 #include "FittingFunction.h"
00009 #include "FileWriter.h"
00010 #include <iostream>
00012
00014
00015
00016
00017 FittingFunction::FittingFunction()
00018 {
00019 }
00020
00021
00022 FittingFunction::~FittingFunction()
00023 {
00024
00025 }
00026
00027
00028
00029 double FittingFunction::getParameter(std::string sParameterName) const{
00030
00031 STRING2INT::const_iterator MyIterator;
00032 if((MyIterator = _mFunctionMap.find(sParameterName)) != _mFunctionMap.end()) {
00033 return _pParameters[(*MyIterator).second];
00034 }
00035 else return 0;
00036 }
00037
00038
00039 double FittingFunction::getParameterError(std::string sParameterName) const {
00040
00041 STRING2INT::const_iterator MyIterator;
00042 if((MyIterator = _mFunctionMap.find(sParameterName)) != _mFunctionMap.end()) {
00043 return _pParametersErrors[(*MyIterator).second];
00044 }
00045 else return 0;
00046 }
00047
00048
00049 double FittingFunction::getParameter(int iIndiceParameter) const{
00050
00051 if ((0 <= iIndiceParameter) && (iIndiceParameter < _iNbParameters)) {
00052 return _pParameters[iIndiceParameter];
00053 }
00054 return 0;
00055 }
00056
00057
00058 double* FittingFunction::getAllParameters() const{
00059
00060 return _pParameters;
00061 }
00062
00063
00064 double FittingFunction::getParameterError(int iIndiceParameter) const {
00065
00066 if ((0 <= iIndiceParameter) && (iIndiceParameter < _iNbParameters)) {
00067 return _pParametersErrors[iIndiceParameter];
00068 }
00069 return 0;
00070 }
00071
00072
00073 double* FittingFunction::getAllParametersErrors() const {
00074 return _pParametersErrors;
00075 }
00076
00077
00078 double FittingFunction::getStandardDeviation()
00079 {
00080 return _dStandardDeviation;
00081 }
00082
00083
00084 double FittingFunction::computeStandardDeviation()
00085 {
00086
00087 size_t lNbdata = _cData->n;
00088
00089 double dSum = 0.0;
00090 double dSumSquare = 0.0;
00091 double dXValue;
00092 for (unsigned int i=0;i<lNbdata;i++)
00093 {
00094 dXValue = *(_cData->x+i);
00095 dSum += dXValue;
00096 dSumSquare += dXValue*dXValue;
00097 }
00098
00099 double dInvNbData = 1/lNbdata;
00100
00101 double dSqrt = dInvNbData*dSumSquare - dInvNbData*dInvNbData*dSum*dSum;
00102
00103 _dStandardDeviation = sqrt(dSqrt);
00104
00105 return getStandardDeviation();
00106
00107 }
00108
00109 double FittingFunction::computeMean(gsl_vector* vVector)
00110 {
00111 double dMean = 0.0;
00112 double dSum = 0.0;
00113 size_t dNbData = vVector->size;
00114
00115 for (unsigned int i=0;i<dNbData;i++)
00116 {
00117 dSum+= gsl_vector_get(vVector,i);
00118 }
00119 dMean = dSum / dNbData;
00120 return dMean;
00121 }
00122
00123 double FittingFunction::computeVariance(gsl_vector* vVector)
00124 {
00125 double dMean = computeMean(vVector);
00126 double dSum = 0.0;
00127 double dDiff = 0.0;
00128 double dDeviation = 0.0;
00129 size_t dNbData = vVector->size;
00130
00131 for (unsigned int i=0;i<dNbData;i++)
00132 {
00133 dDiff=gsl_vector_get(vVector,i) - dMean;
00134 dSum+= dDiff*dDiff;
00135 }
00136 dDeviation = dSum;
00137 return dDeviation;
00138 }
00139
00140 double FittingFunction::computeResidual(gsl_vector* vVectorExp,gsl_vector* vVectorPred)
00141 {
00142 double dSum = 0.0;
00143 double dDiff = 0.0;
00144 double dResidual = 0.0;
00145 size_t dNbData = vVectorExp->size;
00146
00147 for (unsigned int i=0;i<dNbData;i++)
00148 {
00149 dDiff=gsl_vector_get(vVectorExp,i) - gsl_vector_get(vVectorPred,i);
00150 dSum+= dDiff*dDiff;
00151 }
00152 dResidual = dSum;
00153 return dResidual;
00154 }
00155
00156 double FittingFunction::computeDeterminationCoefficient(gsl_vector* vExperimentalDataY,gsl_vector* vFittedDataY)
00157 {
00158 double dDeviation = computeVariance(vExperimentalDataY);
00159 double dResidual = computeResidual(vExperimentalDataY,vFittedDataY);
00160
00161 double dRSQ = 1 - dResidual/dDeviation;
00162 return dRSQ;
00163 }
00164
00165
00166 double FittingFunction::computeObservedFStatistic(gsl_vector* vExperimentalDataY,gsl_vector* vFittedDataY)
00167 {
00168 double dDeviation = computeVariance(vExperimentalDataY);
00169 double dResidual = computeResidual(vExperimentalDataY,vFittedDataY);
00170 size_t lNbData = vExperimentalDataY->size;
00171
00172 double dSSR = dDeviation - dResidual;
00173 double dNbParameters = getNbParameters();
00174 double dDegreeOfFreedom = lNbData - 1 - dNbParameters;
00175
00176 double dMSE = dResidual/dDegreeOfFreedom;
00177 double dMSR = dSSR/dNbParameters;
00178 double dFStatistic = dMSR/dMSE;
00179
00180 return dFStatistic;
00181 }
00182
00183
00184 int FittingFunction::getIndexValue(double dWantedValue,gsl_vector* vVector)
00185 {
00186 size_t lVectorSize = vVector->size;
00187 double lHalfVectorSize = 0.5*lVectorSize;
00188
00189 size_t iIndexMax = gsl_vector_max_index(vVector);
00190 double dMaxValue = gsl_vector_get(vVector, iIndexMax);
00191
00192 double dCurrentValue = dMaxValue;
00193 double dLowValue = dMaxValue;
00194 double dHighValue = dMaxValue;
00195
00196 size_t iIndexHigh = iIndexMax;
00197 size_t iIndexLow = iIndexMax;
00198
00199 size_t iCurrentIndex = iIndexMax;
00200 size_t iGoodIndex = iIndexMax;
00201
00202 double dLowStep = dLowValue;
00203 double dHighStep = dHighValue;
00204
00205
00206
00207 if (iIndexMax < lHalfVectorSize)
00208 {
00209
00210 while(iCurrentIndex < lVectorSize - 1 && dCurrentValue > dWantedValue)
00211 {
00212 dCurrentValue = gsl_vector_get(vVector, ++iCurrentIndex);
00213 iIndexHigh = iCurrentIndex;
00214 iIndexLow = iCurrentIndex - 1;
00215 dLowValue = gsl_vector_get(vVector, iIndexLow);
00216 dHighValue = gsl_vector_get(vVector, iIndexHigh);
00217
00218
00219
00220 }
00221 dLowStep = ::fabs(dLowValue - dWantedValue);
00222 dHighStep = ::fabs(dHighValue - dWantedValue);
00223 if (dLowStep < dHighStep) iGoodIndex = iIndexLow;
00224 else iGoodIndex = iIndexHigh;
00225 }
00226 else
00227 {
00228
00229 while( iCurrentIndex > 0 && dCurrentValue > dWantedValue)
00230 {
00231 dCurrentValue = gsl_vector_get(vVector, --iCurrentIndex);
00232 iIndexLow = iCurrentIndex;
00233 iIndexHigh = iCurrentIndex + 1;
00234 dLowValue = gsl_vector_get(vVector, iIndexLow);
00235 dHighValue = gsl_vector_get(vVector, iIndexHigh);
00236
00237 }
00238 dLowStep = ::fabs(dLowValue - dWantedValue);
00239 dHighStep = ::fabs(dHighValue - dWantedValue);
00240 if (dLowStep < dHighStep) iGoodIndex = iIndexLow;
00241 else iGoodIndex = iIndexHigh;
00242
00243 }
00244
00245 if (iCurrentIndex != iIndexMax)
00246 {
00247
00248 return (int)iGoodIndex;
00249
00250 }
00251 else
00252 {
00253
00254 cout << "\nFind nearest index for " << dWantedValue << " --> " << -1 << endl;
00255 return -1;
00256 }
00257
00258
00259
00260 return -1;
00261 }
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453