Monte Carlo simulation for Uncertainty propagation

Hello ROOTers,
I’ve a question which you might know how to answer.
I’ve a very complex fit with it’s parameter and the covariance matrix (all
come from TMINUIT).

To estimate the uncertainty of the total function (actually on its extrapolation)
I want to draw random parameters (from the best fit values of the parameters and taking into account their error) and make many realizations of the fitting function. Then
I’d be able to determine the uncertainty. Since the parameters are correlated
(here the covariance matrix comes into play) the parameters cannot be extracted
as independent but one would need to define a multi-variate distribution (I guess).

Has any of you done anything similar? Or know where I can find something? Any help is highly appreciated.

ciao
marco

Hi,

Just to make sure we’re on the same page: I think you want to throw correlated random numbers from the covariance matrix from TMinuit (yes, I do it “all the time.” :smiley: ). CERNLib had this functionality built in. Here’s the code (from a friend of mine). It unfortunately assumes NP as the number of parameters, but it should be pretty easy to generalize. The formatting’s a little messed up, but it should be good enough…

Cheers,
Charles


typedef std::vector< double > DVec

// take a matrix and return its sqare root
void corset( const TMatrixD& V, TMatrixD& C );

// Given the sqrt of the covariance matrix, generate correlated random numbers.  This assumes the mean of all variables is 0.  If it is not, you'll 
void corgen( const TMatrixD& C, DVec& x );

void yourRoutine (const TMatrix &covarMatrix)
{
      TMatrixD sqrtCov( NP, NP );
      corset( cov, sqrtCov );

      for (int loop = 0; loop < 10000; ++loop)
      {
          vector< double > values (0, NP);
          corgen (sqrtCov, values);
          // values assumes all means are at 0.  If they are not, you'll have too add the means to each element now.

         // use values
      } // for loop
}


void corset( const TMatrixD& V, TMatrixD& C )
{ 
  // calculate sqrt(V) as lower diagonal matrix
  for( int i = 0; i < NP; ++i ) {
    for( int j = 0; j < NP; ++j ) {
      C[i][j] = 0;
    }
  }

  for( int j = 0; j < NP; ++j ) {
    // diagonal terms first
    double Ck = 0;
    for( int k = 0; k < j; ++k ) {
      Ck += C[j][k] * C[j][k];
    } // k 
    C[j][j] = sqrt( fabs( V[j][j] - Ck ) );

    // off-diagonal terms
    for( int i = j+1; i < NP; ++i ) {
      Ck = 0;
      for( int k = 0; k < j; ++k ) {
	Ck += C[i][k] * C[j][k];
      } //k
      C[i][j] = ( V[i][j] - Ck ) / C[j][j];
    }// i 
  } // j
}

void corgen( const TMatrixD& C, DVec& x )
{
  DVec z (0, NP);

  // np random numbers from unit Gaussian
  for( int i = 0; i < NP; ++i ) {
    z[i] = gRandom->Gaus( 0.0, 1.0 );
  }

  for( int i = 0; i < NP; ++i ) {
    x[i] = 0;
    for( int j = 0; j <= i; ++j ) {
      x[i] += C[i][j] * z[j];
    } // j
  } // i

}

Hi Charles (and all),
unfortunately I believe that’s only part of the story (see also randomizePar in RooFit). I need to draw random correlated numbers from a covariance matrix.I expect these numbers to have a correlation and a dispersion (<-- important) as the observed one. In the case above I believe I would get and rms=1.0 for all the parameters (as for RooFit0. Their means and correlation factors would be correct though.
Any suggestions ?
Thanks and ciao
Marco

Hi,

with the code above you will get a mean of zero, but the RMS (standard deviation) for each variable will be given by the sqrt of the diagonal element of the given correlation matrix.

But if you are interested in getting the uncertainty in the fitting function you don’t need to throw random numbers, you can use the method
TVirtualFitter::GetConfidenceIntervals , see

root.cern.ch/root/htmldoc/TFitte … eIntervals

See also the tutorials $ROOTSYS/tutorisals/fit/ConfidenceIntervals.C

Regards

Lorenzo