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

double
HMM::forwardProcedureCached() {

  //initialisation
  alpha_.col(0) = arma::trans(pi_ % B_.row(0));

  c_(0) = arma::accu(alpha_.col(0));
  alpha_.col(0) /= arma::as_scalar(c_(0));

  //alpha_.print("alpha");
  //c_.print("scale");
  //iteration
  for(unsigned int t = 1; t < T_; ++t) {
    alpha_.col(t) = (A_.t() * alpha_.col(t-1)) % arma::trans(B_.row(t));
    c_(t) = arma::accu(alpha_.col(t));
    alpha_.col(t) /= arma::as_scalar(c_(t)); 
  }

  pprob_ = arma::accu(arma::log(c_));
  return pprob_;
}

arma::vec dmvn_arma(arma::mat x,  
                          arma::mat mean,  
                          arma::mat sigma, 
                          bool logd = false) { 
      int n = x.n_rows;
      int xdim = x.n_cols;
      arma::vec out(n);
      arma::mat rooti = arma::trans(arma::inv(trimatu(arma::chol(sigma))));
      double rootisum = arma::sum(log(rooti.diag()));
      double constants = -(static_cast<double>(xdim)/2.0) * log2pi;
      
      for (int i=0; i < n; i++) {
        arma::vec z = rooti * arma::trans( x.row(i) - mean.row(i)) ;    
        out(i)      = constants - 0.5 * arma::sum(z%z) + rootisum;     
      }  
      
      if (logd == false) {
        out = exp(out);
      }
      return(out);
    }

// [[Rcpp::export]]
arma::vec MCMCloglikelihoodNegBinomCpp_t(const arma::vec& beta, const arma::mat& sigma, double alpha, const arma::vec& sigmaType, const arma::mat& u, 
const arma::vec& df, const arma::vec& kKi, const arma::vec& kLh, const arma::vec& kLhi, const arma::vec& kY, const arma::mat& kX, const arma::mat& kZ) {
  int kM = u.n_rows;  /** Number of MCMC samples */
  arma::vec loglike(kM);
  loglike.fill(0);
  
  for (int i = 0; i < kM; i++) {
    loglike(i) = loglikelihoodNegBinomCpp_t(beta, sigma, alpha, sigmaType, u.row(i).t(), df, kKi, kLh, kLhi, kY, kX, kZ);
  }
  
  return loglike;
}

// [[Rcpp::export]]
arma::vec MCMCloglikelihoodGammaCpp_n(const arma::vec& beta, const arma::mat& sigma, double alpha, 
const arma::mat& u, const arma::vec& kY, const arma::mat& kX, const arma::mat& kZ) {
  int kM = u.n_rows;
  arma::vec loglike(kM);
  loglike.fill(0);
  
  for (int i = 0; i < kM; i++) {
    loglike(i) = loglikelihoodGammaCpp_n(beta, sigma, alpha, u.row(i).t(), kY, kX, kZ);
  }
  
  return loglike;
}

// [[Rcpp::export]]
Rcpp::List CTCRWNLL(const arma::mat& y, const  arma::mat& Hmat, 
                    const arma::vec& beta, const arma::vec& sig2, const arma::vec& delta,
                    const arma::vec& noObs,const arma::vec& active, const arma::colvec& a,
                    const arma::mat& P)
{
  // Define fixed matrices
  int N = y.n_rows;
  double detF;
  arma::mat Z(2,4, fill::zeros);
  Z(0,0) = 1; Z(1,2) = 1;
  arma::mat T(4,4);
  arma::mat Q(4,4);
  arma::mat F(2,2, fill::zeros);
  arma::mat H(2,2, fill::zeros);
  arma::mat K(4,2, fill::zeros);
  arma::mat L(4,4, fill::zeros);
  arma::colvec v(2, fill::zeros);
  arma::colvec aest(4);
  aest=a;
  arma::mat Pest(4,4);
  Pest=P;
  double ll=0;
  //Begin Kalman filter
  for(int i=0; i<N; i++){
    T = makeT(beta(i), delta(i), active(i));
    Q = makeQ(beta(i), sig2(i), delta(i), active(i));
    // prediction
    if(noObs(i)==1){
      aest = T*aest;
      Pest = T*Pest*T.t() + Q;
    } else {
      H(0,0) = Hmat(i,0); 
      H(1,1) = Hmat(i,1); 
      H(0,1) = Hmat(i,2); 
      H(1,0) = Hmat(i,2);
      v = y.row(i).t()-Z*aest;
      F = Z*Pest*Z.t() + H; 
      detF = F(0,0)*F(1,1) - F(1,0)*F(0,1);
      if(detF<=0){
        aest = T*aest;
        Pest = T*Pest*T.t() + Q;
      } else{
        ll += - (log(detF) + dot(v,solve(F,v)))/2; 
        // K = T*Pest*Z.t()*inv_sympd(F);  
        K = T*Pest*Z.t()*F.i();
        L = T - K*Z;
        aest = T*aest + K*v;
        Pest = T*Pest*L.t() + Q;
      }
    }
  }
  return Rcpp::List::create(Rcpp::Named("ll") = ll);
}

void GMM::sample(arma::mat& X){
    std::size_t k;
    arma::vec tmp;
    dist = boost::random::discrete_distribution<>(pi);
    for(std::size_t i = 0; i < X.n_rows;i++){

        k  =  dist(generator);
        A[k]*arma::randn<arma::vec>(D);
        tmp = Means[k] + A[k]*arma::randn<arma::vec>(D);
        X.row(i) = tmp.st();
    }
}

List FFBS(const arma::mat& allprobs, const arma::vec& delta,
          const arma::mat& mGamma, const int& K, const int& T) {
  arma::mat lalpha = zeros(K, T);
  arma::mat lbeta = zeros(K, T);

  arma::vec foo(K);
  double sumfoo, lscale;
  int i;

  foo = delta % allprobs.row(0).t();
  sumfoo = sum(foo);
  lscale = log(sumfoo);
  foo = foo / sumfoo;

  lalpha.col(0) = log(foo) + lscale;
  for (i = 1; i < T; i++) {
    foo = (foo.t() * mGamma).t() % allprobs.row(i).t();
    sumfoo = sum(foo);
    lscale = lscale + log(sumfoo);
    foo = foo / sumfoo;
    lalpha.col(i) = log(foo) + lscale;
  }
  for (i = 0; i < K; i++) {
    foo(i) = 1.0 / K;
  }
  lscale = log((double)K);
  for (i = T - 2; i >= 0; i--) {
    foo = mGamma * (allprobs.row(i + 1).t() % foo);
    lbeta.col(i) = log(foo) + lscale;
    sumfoo = sum(foo);
    foo = foo / sumfoo;
    lscale = lscale + log(sumfoo);
  }

  List FS;
  FS["lalpha"] = lalpha;
  FS["lbeta"] = lbeta;

  return FS;
}

// [[Rcpp::export]]
arma::mat gillespie_next(arma::vec theta, arma::vec init_state, arma::mat trans, LogicalMatrix depend, int t_end, bool info = false){
  
  int n_event = trans.n_rows;
  int j_event;
  
  Progress p(0, false); //progress to check for user interrupt
  
  //initialize trace of Monte Carlo simulation
  arma::mat trace = arma::zeros(1,init_state.n_elem);
  trace.row(0) = init_state.t();
  arma::vec current_state;
  arma::vec current_rates;
  
  //main simulation loop
  double time = 0.0;
  int i = 1;
  
  while(time <= t_end){
    
    //check for user abort
    if(Progress::check_abort()){
      Rcout << "User abort detected at time: " << time << ", i: " << i << std::endl;
      return(arma::zeros(2,2));
    }
    
    if(info){ //print simulation details
      Rcout << "time: " << time << ", i: " << i << std::endl;  
    }
    
    current_state = trace.row(i-1).t(); //extract state at beginning of time jump
    
    if(i == 1){
      current_rates = rate_wrapper(theta,current_state,seirs_demography_rate); //calculate initial rates   
    } else {
      current_rates = update_rates_seirs_demography(theta,current_rates,current_state,depend,j_event); //update rates based on dependency graph
    }
    
    arma::vec tau(n_event); //calculate vector of times to next event
    for(int j=0; j<n_event; j++){
      tau(j) = (-1 / current_rates(j)) * log(R::runif(0,1));
    }

    j_event = tau.index_min(); //fire j_event
    current_state = current_state + trans.row(j_event).t(); //update the current state
    trace.insert_rows(i,current_state.t());
    
    time = time + tau(j_event); //update time
    i++; //update iterator
  }
  
  return(trace);
}

inline void SVDIncompleteIncrementalLearning::
                              HUpdate<arma::sp_mat>(const arma::sp_mat& V,
                                                    const arma::mat& W,
                                                    arma::mat& H)
{
  arma::mat deltaH(H.n_rows, 1);
  deltaH.zeros();

  for(arma::sp_mat::const_iterator it = V.begin_col(currentUserIndex);
                                        it != V.end_col(currentUserIndex);it++)
  {
    double val = *it;
    size_t i = it.row();
    if((val = V(i, currentUserIndex)) != 0)
      deltaH += (val - arma::dot(W.row(i), H.col(currentUserIndex))) *
                                                    arma::trans(W.row(i));
  }
  if(kh != 0) deltaH -= kh * H.col(currentUserIndex);

  H.col(currentUserIndex++) += u * deltaH;
  currentUserIndex = currentUserIndex % m;
}

arma::mat hCoef(const arma::mat& weights, const arma::mat& X) {

  int p = X.n_cols;
  arma::mat hess(p * (weights.n_rows - 1), p * (weights.n_rows - 1), arma::fill::zeros);
  for (unsigned int i = 0; i < X.n_rows; i++) {
    arma::mat XX = X.row(i).t() * X.row(i);
    for (unsigned int j = 0; j < (weights.n_rows - 1); j++) {
      for (unsigned int k = 0; k < (weights.n_rows - 1); k++) {
        if (j != k) {
          hess.submat(j * p, k * p, (j + 1) * p - 1, (k + 1) * p - 1) += XX * weights(j + 1, i)
              * weights(k + 1, i);
        } else {
          hess.submat(j * p, j * p, (j + 1) * p - 1, (j + 1) * p - 1) -= XX * weights(j + 1, i)
              * (1.0 - weights(j + 1, i));

        }

      }
    }
  }
  return hess;
}

//--------------------------------------------------------------------------------------------------
List variogram( const arma::mat& Z,
                const arma::mat& X, 
                Function d ) {
  int i, j, k;
  int n = X.n_rows;
  int N = n * ( n + 1 ) / 2;
  
  std::vector< int > I( N );
  arma::rowvec x, y;
  arma::colvec V( N );
  arma::colvec D( N );
  arma::colvec S( N );
  
  List Varg;
  
  for ( i = 0; i < N; i++ ) { 
    I[i] = i;
  }
  
  k = 0;
  for ( i = 0; i < n; i++ ) { 
    for ( j = i; j < n; j++ ) {
      x = X.row( i );
      y = X.row( j );
      
      D(k) = as<double>( d( x, y ) );
      S(k) = ( Z(i,0) - Z(j,0) ) * ( Z(i,0) - Z(j,0) );
    
      k++;
    }
  }
  
  sort( I.begin(), I.end(),
       [&]( const int& k, const int& l ) {
         return ( D(k) < D(l) );
       }
  );
  
  V( 0 ) = S( I[0] );
  for ( k = 1; k < N; k++ ) {
    V( k ) = ( V( k - 1 ) * k + S( I[ k ] ) ) / ( k + 1.0 );
  }
  V = 0.5 * V;
  
  Varg[ "variogram" ] = V;
  Varg[ "distance" ] = D;
  Varg[ "sort" ] = I;
  
  return Varg;
}

RegularizedSVDFunction::RegularizedSVDFunction(const arma::mat& data,
                                               const size_t rank,
                                               const double lambda) :
    data(data),
    rank(rank),
    lambda(lambda)
{
  // Number of users and items in the data.
  numUsers = max(data.row(0)) + 1;
  numItems = max(data.row(1)) + 1;

  // Initialize the parameters.
  initialPoint.randu(rank, numUsers + numItems);
}

// [[Rcpp::export]]
arma::mat gillespie_direct(arma::vec theta, arma::vec init_state, arma::mat trans, int t_end, bool info = false){
  
  Progress p(0, false); //progress to check for user interrupt
  
  //initialize trace of Monte Carlo simulation
  arma::mat trace = arma::zeros(1,init_state.n_elem);
  trace.row(0) = init_state.t();
  arma::vec current_state;
  
  //main simulation loop
  double time = 0.0;
  int i = 1;
  
  while(time <= t_end){
    
    //check for user abort
    if(Progress::check_abort()){
      Rcout << "User abort detected at time: " << time << ", i: " << i << std::endl;
      return(arma::zeros(2,2));
    }
    
    if(info){ //print simulation details
      Rcout << "time: " << time << ", i: " << i << std::endl;  
    }
    
    current_state = trace.row(i-1).t(); //extract state at beginning of time jump
    arma::vec current_rates = rate_wrapper(theta,current_state,sir_rate); //calculate current rates
    
    double w0 = sum(current_rates); //sum of rate (propensity) functions
    double tau = 1/w0 * log(1/R::runif(0,1)); //calculate time to next reaction
    
    double r_num = R::runif(0,1);
    double w0_rxn = r_num * w0; //instantaneous event probabilities
    
    //decide which event j occured
    int j = 0;
    while(sum(current_rates.subvec(0,j)) < w0_rxn){
      j++;
    }
    
    current_state = current_state + trans.row(j).t(); //update the current state
    trace.insert_rows(i,current_state.t());
    
    time = time + tau; //update time
    i++; //update iterator
  }
  
  return(trace);
}

arma::uvec generate_categoricals(arma::uvec zlabels, arma::mat probs) {
    int ndata = zlabels.n_elem;
    int nclusters = probs.n_rows;
    
    std::vector<boost::random::discrete_distribution<> > cat_dists;
    for (int k=0; k<nclusters; k++) {
        std::vector<double> probs_k = arma::conv_to<std::vector<double> >::from(probs.row(k));
        cat_dists.push_back(boost::random::discrete_distribution<> (probs_k.begin(), probs_k.end()));
    }
    arma::uvec categories(ndata);
    for (int i=0; i<ndata; i++) {
        categories(i) = cat_dists[zlabels(i)](rng);
    }
    return categories;
}

// [[Rcpp::export]]
arma::mat gillespie_first(arma::vec theta, arma::vec init_state, arma::mat trans, int t_end, bool info = false){

  int n_event = trans.n_rows;
  
  Progress p(0, false); //progress to check for user interrupt
  
  //initialize trace of Monte Carlo simulation
  arma::mat trace = arma::zeros(1,init_state.n_elem);
  trace.row(0) = init_state.t();
  arma::vec current_state;
  
  //main simulation loop
  double time = 0.0;
  int i = 1;
  
  while(time <= t_end){
    
    //check for user abort
    if(Progress::check_abort()){
      Rcout << "User abort detected at time: " << time << ", i: " << i << std::endl;
      return(arma::zeros(2,2));
    }
    
    if(info){ //print simulation details
      Rcout << "time: " << time << ", i: " << i << std::endl;  
    }
    
    current_state = trace.row(i-1).t(); //extract state at beginning of time jump
    arma::vec current_rates = rate_wrapper(theta,current_state,sir_rate); //calculate current rates
    
    arma::vec rand; //create vector of U[0,1]
    rand = as<arma::vec>(runif(n_event));
    
    arma::vec tau(n_event); //calculate vector of times to next event
    for(int j=0; j<n_event; j++){
      tau(j) = (-1 / current_rates(j)) * log(rand(j));
    }
    
    int first_rxn = tau.index_min(); //which even occurs first
    current_state = current_state + trans.row(first_rxn).t(); //update the current state
    trace.insert_rows(i,current_state.t());
    
    time = time + tau(first_rxn); //update time
    i++; //update iterator
  }
  
  return(trace);
}