src/fisher.cpp
// Computes the full covariance matrix for the MLE.
Computes the full covariance matrix for the MLE.
Computes the full covariance matrix for the MLE .
void full_MLE_covariance_matrix(Tree &T, Model &Mod, Parameters &Par, long N, Array2 &Cov) {
  long i, j;
  long npars = Mod.df*T.nedges + Mod.rdf;
  ublas::matrix<double> CC(npars, npars);
  ublas::matrix<double> II(npars, npars);
  Array2 I;
  Fisher_information(T, Mod, Par, N, I);
  for(i=0; i < npars; i++) {
    for(j=0; j < npars; j++) {
      II(i, j) = I[i][j];
    }
  }
  bool res = matrix_inverse(II, CC);
  if (!res) {
    throw std::out_of_range( "Could not invert the Fisher information matrix." );
  }
  Cov.resize(npars);
  for(i=0; i < npars; i++) {
    Cov[i].resize(npars);
    for(j=0; j < npars; j++) {
      Cov[i][j] = CC(i,j);
    }
  }
}
void full_MLE_covariance_matrix ( Tree & T , Model & Mod , Parameters & Par , long N , Array2 & Cov ) { long i , j ; long npars = Mod . df * T . nedges + Mod . rdf ; ublas :: matrix < double > CC ( npars , npars ) ; ublas :: matrix < double > II ( npars , npars ) ; Array2 I ; Fisher_information ( T , Mod , Par , N , I ) ; for ( i = 0 ; i < npars ; i ++ ) { for ( j = 0 ; j < npars ; j ++ ) { II ( i , j ) = I [ i ] [ j ] ; } } bool res = matrix_inverse ( II , CC ) ; if ( ! res ) { throw std :: out_of_range ( " " ) ; } Cov . resize ( npars ) ; for ( i = 0 ; i < npars ; i ++ ) { Cov [ i ] . resize ( npars ) ; for ( j = 0 ; j < npars ; j ++ ) { Cov [ i ] [ j ] = CC ( i , j ) ; } } }
Computes the full covariance matrix for the MLE.
Computes the full covariance matrix for the MLE .
// Need to convert from Matrix to ublas::matrix<double>