ReClaM: VarianceEstimator.cpp Source File

00001 //===========================================================================
00062 //===========================================================================
00063 
00064 
00065 #include <SharkDefs.h>
00066 #include <Array/ArrayOp.h>
00067 #include <LinAlg/LinAlg.h>
00068 #include <ReClaM/VarianceEstimator.h>
00069 
00070 
00071 //===========================================================================
00123 void VarianceEstimator::estimateFisherInformation
00124 (
00125     Model& model,
00126     const Array< double >& inputA,
00127     const Array< double >& outputA,
00128     Array< double >&       infMatA
00129 )
00130 {
00131     SIZE_CHECK(inputA.ndim() == 2)
00132 
00133     double          mseL = 0;
00134     Array< double > dwL;
00135     Array< double > modelOutputL(outputA.dim(1));
00136 
00137     for (unsigned i = 0; i < inputA.dim(0); ++i)
00138     {
00139         model.modelDerivative(inputA[ i ], modelOutputL, dwL);
00140         mseL += sqrDistance(outputA[ i ], modelOutputL);
00141 
00142         dwL.resize(dwL.nelem(), true);
00143 
00144         if (i == 0)
00145         {
00146             infMatA = outerProduct(dwL, dwL);
00147         }
00148         else
00149         {
00150             infMatA += outerProduct(dwL, dwL);
00151         }
00152     }
00153 
00154     mseL    /= inputA.dim(0);
00155     infMatA /= mseL;
00156 }
00157 
00158 
00159 //===========================================================================
00210 void VarianceEstimator::estimateInvFisher
00211 (
00212     Model& model,
00213     const Array< double >& inputA,
00214     const Array< double >& outputA,
00215     Array< double >&       invInfMatA,
00216     Array< double >&       transInvInfMatA,
00217     double&                s2A
00218 )
00219 {
00220     SIZE_CHECK(inputA.ndim() == 2)
00221 
00222     Array< double > infMatL;
00223     Array< double > dwL;
00224     Array< double > dwMeanL;
00225     Array< double > modelOutputL(outputA.dim(1));
00226 
00227     s2A = 0;
00228     for (unsigned i = 0; i < inputA.dim(0); ++i)
00229     {
00230         model.modelDerivative(inputA[ i ], modelOutputL, dwL);
00231         s2A += sqrDistance(outputA[ i ], modelOutputL);
00232 
00233         dwL.resize(dwL.nelem(), true);
00234 
00235         if (i == 0)
00236         {
00237             infMatL = outerProduct(dwL, dwL);
00238             dwMeanL = dwL;
00239         }
00240         else
00241         {
00242             infMatL += outerProduct(dwL, dwL);
00243             dwMeanL += dwL;
00244         }
00245     }
00246 
00247     s2A     /= inputA.dim(0);
00248     infMatL /= s2A;
00249 
00250     invInfMatA = invert(infMatL);
00251 
00252     //
00253     // innerProduct yields the same results as matrixProduct (not
00254     // implemented yet) since all matrices are symmetric
00255     //
00256     transInvInfMatA =
00257         innerProduct(invInfMatA,
00258                      innerProduct(outerProduct(dwMeanL, dwMeanL),
00259                                   invInfMatA));
00260 }
00261 
00262 
00263 //===========================================================================
00314 double VarianceEstimator::overallVariance
00315 (
00316     Model& model,
00317     const Array< double >& inputA,
00318     const Array< double >& outputA
00319 )
00320 {
00321     SIZE_CHECK(inputA.ndim() == 2)
00322 
00323     double          mseL = 0;
00324     Array< double > infMatL;
00325     Array< double > invInfMatL;
00326     Array< double > invAdwMeanL;
00327     Array< double > dwL;
00328     Array< double > dwMeanL;
00329     Array< double > modelOutputL(outputA.dim(1));
00330 
00331     for (unsigned i = 0; i < inputA.dim(0); ++i)
00332     {
00333         model.modelDerivative(inputA[ i ], modelOutputL, dwL);
00334         mseL += sqrDistance(outputA[ i ], modelOutputL);
00335 
00336         dwL.resize(dwL.nelem(), true);
00337 
00338         if (i == 0)
00339         {
00340             infMatL = outerProduct(dwL, dwL);
00341             dwMeanL = dwL;
00342         }
00343         else
00344         {
00345             infMatL += outerProduct(dwL, dwL);
00346             dwMeanL += dwL;
00347         }
00348     }
00349 
00350     mseL       /= inputA.dim(0);
00351     infMatL    /= mseL;
00352     invInfMatL  = invert(infMatL);
00353     invAdwMeanL = innerProduct(invInfMatL, outerProduct(dwMeanL, dwMeanL));
00354 
00355     return scalarProduct(dwMeanL, innerProduct(invInfMatL, dwMeanL))
00356            + trace(invAdwMeanL);
00357 }
00358 
00359 
00360 //===========================================================================
00394 double VarianceEstimator::estimateVariance
00395 (
00396     Model& model,
00397     const Array< double >& inputA,
00398     const Array< double >& invInfMatA
00399 )
00400 {
00401     Array< double > dwL;
00402     model.modelDerivative(inputA, dwL);
00403     dwL.resize(dwL.nelem(), true);
00404     return scalarProduct(dwL, innerProduct(invInfMatA, dwL));
00405 }
00406 
00407 
00408 //===========================================================================
00450 double VarianceEstimator::estimateVarianceChange
00451 (
00452     Model& model,
00453     const Array< double >& inputA,
00454     const Array< double >& invInfMatA,
00455     const Array< double >& transInvInfMatA,
00456     double                 s2A
00457 )
00458 {
00459     Array< double > dwL;
00460     model.modelDerivative(inputA, dwL);
00461     dwL.resize(dwL.nelem(), true);
00462     return scalarProduct(dwL, innerProduct(transInvInfMatA, dwL))
00463            / (s2A + scalarProduct(dwL, innerProduct(invInfMatA, dwL)));
00464 }
00465