src/fisher.cpp
// Computes the full covariance matrix for the MLE, using observed Fisher info.
Computes the full covariance matrix for the MLE, using observed Fisher info.
Computes the full covariance matrix for the MLE using observed Fisher info .
void full_MLE_observed_covariance_matrix(Tree &T, Model &Mod, Parameters &Par, Counts &data, 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;
  Observed_Fisher_information(T, Mod, Par, data, 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_observed_covariance_matrix ( Tree & T , Model & Mod , Parameters & Par , Counts & data , 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 ; Observed_Fisher_information ( T , Mod , Par , data , 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, using observed Fisher info.
Computes the full covariance matrix for the MLE using observed Fisher info .
// Need to convert from Matrix to ublas::matrix<double>