arma_basics.cpp

/*********************************************************************
 * 
 * arma_basics.cpp
 * Includes simple examples for learning how to use RcppArmadillo
 * Some material taken from Dirk Eidelbuettel's wabpage:
 *      http://dirk.eddelbuettel.com/code/rcpp.armadillo.html
 * Philip Barrett, Chicago, 01may2014
 * 
 *********************************************************************/


#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
mat identity_arma(int iNN) {
// Create & return an identity matrix
  mat AA = eye( iNN, iNN ) ;
  return AA ;
}

// [[Rcpp::export]]
mat invert_arma(mat XX) {
// Invert a matrix
  mat YY = inv( XX ) ;
  return YY ;
}

// [[Rcpp::export]]
colvec solve_cpp( mat mA, colvec vB ) {
// Invert a matrix
  return solve( mA, vB ) ;
}

// [[Rcpp::export]]
List lm_cpp( colvec vY, mat mX) {
// How to do linear projection in C++
  int iObs = mX.n_rows, iVars = mX.n_cols;
      // The number of variables and observations
  colvec vCoef = solve( mX, vY );      // fit model y ~ X
  colvec vResid = vY - mX * vCoef;            // residuals
  double dSig2 = as_scalar( trans( vResid ) * vResid / ( iObs - iVars ) );
      // Var-Covar matrix
  colvec vStdErrEst = sqrt( dSig2 * diagvec( inv( trans( mX ) * mX ) ) );
      // std.error of estimate
  return List::create( Named("coefficients") = vCoef, Named("stderr") = vStdErrEst ) ;
}

// [[Rcpp::export]]
mat arma_ops( mat mX ) {
// Matrix assingment. mX should be 3x3
  
  // Error checking
//  if( mX.n_rows != 3 || mX.n_cols != 3 )
  if( !( mX.n_rows == 3 && mX.n_cols == 3 ) )
    stop("Input must be a 3x3 matrix") ;
  
  // A bunch of useful operations
  mat mY ;
      // Instantiation
  mY << 1 << 2 << 3 << endr
          << 4 << 5 << 6 << endr
          << 7 << 8 << 9 << endr ;
      // Initialization
  Rcout << "mY = \n" << mY << std::endl << std::endl ;
  Rcout << "mX = \n" << mX << std::endl << std::endl ;
      // Printing
  Rcout << "mX[ 1:2, 2:3 ] = \n" << mX.submat( 0, 1, 1, 2)  << std::endl << std::endl ;
      // Submatrix
  mat mZ = zeros( 3, 3 ) ;
  mat mO = ones( 3, 3 ) ;
  mat mI = eye( 3, 3 ) ;
  Rcout << "mX = \n" << mZ << std::endl << std::endl ;
  Rcout << "mX = \n" << mO << std::endl << std::endl ;
  Rcout << "mI = \n" << mI << std::endl << std::endl ;
      // Special matricies
  Rcout << "mX %*% mY = \n" << mX * mY << std::endl << std::endl ;
  Rcout << "mX * mY = \n" << mX % mY << std::endl << std::endl ;
  Rcout << "mX / mY = \n" << mX / mY << std::endl << std::endl ;
  return mY ;
}