C++ mat::i方法代码示例

kalmanFilter::kalmanFilter(arma::vec state, arma::mat transistion, arma::mat meas, arma::mat mNoise, arma::mat pNoise) {
	stateEstimation = state;
	transistionMatrix = transistion;
	transistionMatrixT = transistion.i();
	measurementMatrix = meas;
	transistionMatrixT = meas.i();
	measNoise = mNoise;
	processNoise = pNoise;
	sanityChecks();

	std::cout << "Kalman filter object is created.\n";
}

void QUIC_SVD::ExtractSVD(arma::mat& u,
                          arma::mat& v,
                          arma::mat& sigma)
{
  // Calculate A * V_hat, necessary for further calculations.
  arma::mat projectedMat;
  if (dataset.n_cols > dataset.n_rows)
    projectedMat = dataset.t() * basis;
  else
    projectedMat = dataset * basis;

  // Calculate the squared projected matrix.
  arma::mat projectedMatSquared = projectedMat.t() * projectedMat;

  // Calculate the SVD of the above matrix.
  arma::mat uBar, vBar;
  arma::vec sigmaBar;
  arma::svd(uBar, sigmaBar, vBar, projectedMatSquared);

  // Calculate the approximate SVD of the original matrix, using the SVD of the
  // squared projected matrix.
  v = basis * vBar;
  sigma = arma::sqrt(diagmat(sigmaBar));
  u = projectedMat * vBar * sigma.i();

  // Since columns are sampled, the unitary matrices have to be exchanged, if
  // the transposed matrix is not passed.
  if (dataset.n_cols > dataset.n_rows)
  {
    arma::mat tempMat = u;
    u = v;
    v = tempMat;
  }
}

// [[Rcpp::export]]
arma::vec Mahalanobis(arma::mat x, arma::rowvec center, arma::mat cov){
    int n = x.n_rows;
    arma::mat x_cen;
    x_cen.copy_size(x);
    for (int i=0; i < n; i++) {
        x_cen.row(i) = x.row(i) - center;
    }
    return sum((x_cen * cov.i()) % x_cen, 1);
}

 double mahalanobis(const arma::rowvec& x, const arma::rowvec& mu, const arma::mat& sigma) {
   const arma::rowvec err = x - mu;
   return arma::as_scalar(err * sigma.i() * err.t());
 }

// [[Rcpp::export]]
List PCspatregMCMC(const arma::vec& y, const arma::mat& X, const arma::mat& K, const arma::uvec& sampleIndex,
const arma::vec& betaMean, const arma::mat& betaSig, const double& phiScale, const int& burn, const int& iter) {
  int na = K.n_cols;
  int nsite = X.n_rows;
  int nb = X.n_cols;
  int ns = y.n_elem;
  arma::uvec idx = sampleIndex-1;
  arma::mat Xsmp = X.rows(idx);
  arma::mat Ksmp = K.rows(idx);
  
  arma::mat betaStor(iter, nb, fill::zeros);
  arma::mat alphaStor(iter, na, fill::zeros);
  arma::vec phiStor(iter, fill::zeros);
  arma::vec sigStor(iter, fill::zeros);
  arma::mat predStor(iter, nsite);
  
  // Beta items
  arma::mat betaPrec = betaSig.i();
  arma::mat V_beta_inv(nb, nb);
  arma::vec v_beta(nb);
  // Alpha items
  arma::mat I(na,na, fill::eye);
  arma::mat V_alpha_inv(na, na);
  arma::vec v_alpha(na);
  // phi items
  arma::vec Kbar(ns);
  double V_phi_inv;
  double v_phi;
  
  // Initial values
  arma::vec beta(nb);
  beta = solve(Xsmp.t() * Xsmp, Xsmp.t()*y);
  double tau = (ns-1)/as_scalar((y - Xsmp*beta).t()*(y - Xsmp*beta));
  double phi = 0;
  arma::vec alpha(na, fill::zeros);
  
  for(int i=0; i<iter+burn; i++){
    
    // update beta
    V_beta_inv = tau*Xsmp.t()*Xsmp + betaPrec;
    v_beta = tau * Xsmp.t() * (y - phi*Ksmp*alpha) + betaPrec*betaMean;
    beta = GCN(V_beta_inv, v_beta);
    if(i>=burn){betaStor.row(i-burn) = beta.t();}
    
    // update alpha
    V_alpha_inv = phi*phi*tau*Ksmp.t()*Ksmp + I; 
    v_alpha = phi*tau*Ksmp.t()*(y-Xsmp*beta);
    alpha = GCN(V_alpha_inv, v_alpha);
    if(i>=burn){alphaStor.row(i-burn) = phi*alpha.t();}
    
    // update phi
    Kbar = Ksmp*alpha;
    V_phi_inv = tau*as_scalar(Kbar.t()*Kbar) + 1/(phiScale*phiScale);
    v_phi = tau*dot(Kbar, y-Xsmp*beta);
    phi = as<double>(rnorm(1, v_phi/V_phi_inv, 1/sqrt(V_phi_inv)));
    if(i>=burn) phiStor(i-burn) = phi; 
    
    // update tau
    tau = as<double>(rgamma(1, ns/2, as_scalar((y-Xsmp*beta-phi*Ksmp*alpha).t()*(y-Xsmp*beta-phi*Ksmp*alpha))/2));
    if(i>=burn) sigStor(i-burn) = 1/sqrt(tau);
    
    // make prediction
    if(i>=burn){
      predStor.row(i-burn) = (X*beta + phi*K*alpha).t();
    }
    
    
  }
  
  return Rcpp::List::create(
    Rcpp::Named("beta") = betaStor, 
    Rcpp::Named("alpha") = alphaStor, 
    Rcpp::Named("phi")=phiStor,
    Rcpp::Named("sigma")=sigStor,
    Rcpp::Named("pred")=predStor
    );  
}