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

//' @title Safe softmax function for computing weights
//'
//' @description Computes the weights \eqn{w_i = \frac{e^{p_i}}{\sum_{j=1}^k e^{p_j}}} from \eqn{p_i}, \eqn{i=1,\ldots,k}
//' in a safe way to avoid overflows and to truncate automatically to zero low values of \eqn{w_i}.
//' 
//' @param logs matrix of logarithms where each row contains a set of \eqn{p_1,\ldots,p_k} to compute the weights from.
//' @param etrunc truncation for exponential: \code{exp(x)} with \code{x <= -etrunc} is set to zero. Defaults to \code{30}.
//' @return A matrix of the size as \code{logs} containing the weights for each row.
//' @author Eduardo Garcia-Portugues (\email{[email protected]@math.ku.dk})
//' @details The \code{logs} argument must be always a matrix.
//' @examples
//' # A matrix
//' safeSoftMax(rbind(1:10, 20:11))
//' rbind(exp(1:10) / sum(exp(1:10)), exp(20:11) / sum(exp(20:11)))
//' 
//' # A row-matrix
//' safeSoftMax(rbind(-100:100), etrunc = 30)
//' exp(-100:100) / sum(exp(-100:100))
//' @export
// [[Rcpp::export]]
arma::mat safeSoftMax(arma::mat logs, double etrunc = 30) {
  
  // Maximum of logs by rows to avoid overflows
  arma::vec m = max(logs, 1);
  
  // Recenter by columns
  logs.each_col() -= m;
  
  // Ratios by columns
  logs.each_col() -= log(sum(exp(logs), 1));
  
  // Truncate exponential by using a lambda function - requires C++ 11
  logs.transform([etrunc](double val) { return (val < -etrunc) ? double(0) : double(exp(val)); });
  
  return logs;
  
}

    /**
     * Construct the exact kernel matrix.
     *
     * @param data Input data points.
     * @param transformedData Matrix to output results into.
     * @param eigval KPCA eigenvalues will be written to this vector.
     * @param eigvec KPCA eigenvectors will be written to this matrix.
     * @param rank Rank to be used for matrix approximation.
     * @param kernel Kernel to be used for computation.
     */
    static void ApplyKernelMatrix(const arma::mat& data,
                                  arma::mat& transformedData,
                                  arma::vec& eigval,
                                  arma::mat& eigvec,
                                  const size_t /* unused */,
                                  KernelType kernel = KernelType())
  {
    // Construct the kernel matrix.
    arma::mat kernelMatrix;
    // Resize the kernel matrix to the right size.
    kernelMatrix.set_size(data.n_cols, data.n_cols);

    // Note that we only need to calculate the upper triangular part of the 
    // kernel matrix, since it is symmetric. This helps minimize the number of
    // kernel evaluations.
    for (size_t i = 0; i < data.n_cols; ++i)
    {
      for (size_t j = i; j < data.n_cols; ++j)
      {
        // Evaluate the kernel on these two points.
        kernelMatrix(i, j) = kernel.Evaluate(data.unsafe_col(i),
                                             data.unsafe_col(j));
      }
    }

    // Copy to the lower triangular part of the matrix.
    for (size_t i = 1; i < data.n_cols; ++i)
      for (size_t j = 0; j < i; ++j)
        kernelMatrix(i, j) = kernelMatrix(j, i);

    // For PCA the data has to be centered, even if the data is centered. But it
    // is not guaranteed that the data, when mapped to the kernel space, is also
    // centered. Since we actually never work in the feature space we cannot
    // center the data. So, we perform a "psuedo-centering" using the kernel
    // matrix.
    arma::rowvec rowMean = arma::sum(kernelMatrix, 0) / kernelMatrix.n_cols;
    kernelMatrix.each_col() -= arma::sum(kernelMatrix, 1) / kernelMatrix.n_cols;
    kernelMatrix.each_row() -= rowMean;
    kernelMatrix += arma::sum(rowMean) / kernelMatrix.n_cols;

    // Eigendecompose the centered kernel matrix.
    arma::eig_sym(eigval, eigvec, kernelMatrix);

    // Swap the eigenvalues since they are ordered backwards (we need largest to
    // smallest).
    for (size_t i = 0; i < floor(eigval.n_elem / 2.0); ++i)
      eigval.swap_rows(i, (eigval.n_elem - 1) - i);

    // Flip the coefficients to produce the same effect.
    eigvec = arma::fliplr(eigvec);

    transformedData = eigvec.t() * kernelMatrix;
    transformedData.each_col() /= arma::sqrt(eigval);
  }

  /** Initialization
   *
   * @return Estimated state \f$\{\tilde{\mathbf{x}}^{(i)}_0,\tilde{\omega}^{(i)}\}_{i=1}^{M}\f$
   */
  CompeleteState initialize() {
    state_par_.each_col([this](arma::vec &v) {
      v = process_.template getProcess<0>().getInitialPDF().random();
    });

    unsigned long cnt = 0;
    w_.for_each([this, &cnt](double &e) {
      e = process_.template getProcess<0>().getInitialPDF().likelihood(
          state_par_.col(cnt++));
    });

    normalizeWeights();

    return std::make_tuple(state_par_, w_);
  }

 void predict(const Args &... args) {
   state_par_.each_col([this, &args...](arma::vec &v, const Args &... args) {
     v = process_.template getProcess<0>().getCPDF().random(v, args...);
   });
 }

Rcpp::List irls_binomial_cpp_fast_br(arma::mat A, arma::vec b, double maxit, double tol)
{
//Def
arma::vec x;
x.zeros(A.n_cols,1);
arma::vec xold;
arma::mat varmatrix;

double nobs;
nobs = A.n_rows;
double df;
df = A.n_cols;
double n;

double ll;
double aic;
double bic;
double mdl;



arma::vec W(nobs);
arma::vec unit(nobs);
unit.ones(nobs);
arma::vec eta(nobs);
arma::vec g(nobs);
arma::vec gprime(nobs);
arma::vec z(nobs);
//mod
arma::vec bprime(nobs);

//int k;

for (int i = 0; i < maxit; ++i) {
  eta = A * x;
  
  for (int j=0; j < nobs; ++j) {
    g[j] = 1.0 / (1.0 + exp(-1.0 * eta[j]));
    gprime[j] = exp (-1.0 * eta[j]) / ((1.0 + exp (-1.0 * eta[j])) * (1.0 + exp (-1.0 * eta[j])));
    //mod
    bprime[j] = (b[j]+(df/nobs)*(0.5))/(1+(df/nobs));
    //bprime[j] = b[j]+(sum(b)/nobs)*(0.5);
  }
     
  //z = eta+(b-g)/gprime;
  z = eta+(bprime-g)/gprime;
  
  W = (gprime % gprime);
  W /= (g % (unit-g)); 
  //W += unit;
  xold = x;
  
  //coefficients
  //x = arma::solve(A.t()*(W % A.each_col()), A.t()*(W % z), arma::solve_opts::no_approx);
  varmatrix = A.t()*(W % A.each_col());
  x = arma::solve(varmatrix, A.t()*(W % z), arma::solve_opts::no_approx);
  //k = i;
  
if(sqrt(arma::dot(x-xold,x-xold)) < tol){
 break;
}}

n = A.n_rows;

//arma::vec e;
//double ssr;
//e = (b - A*x);
//e = (bprime - A*x);
//ssr = accu(e.t()*e);

//scores

//ll = arma::accu(-arma::dot(b,log(unit + exp(-(A*x)))) - arma::dot((unit-b),log(unit + exp(A*x))));
ll = arma::accu(-arma::dot(b,log(unit + exp(-(A*x)))) - arma::dot((unit-b),log(unit + exp(A*x))));

aic = - 2 * ll + 2 * df;

bic = - 2 * ll + log(nobs) * df;

//mdl

// arma::mat xz;
// xz.zeros(size(x));
// 
// arma::vec ez;
// double ssrz;
// double ssrtot;
// double RR;
// double F;
// double mdl;
// arma::vec yaverage(n);
// 
// ez = (b - A*xz);
// ssrz = accu(ez.t()*ez);
// F = (((ssrz - ssr)/df)/(ssr/((n-(df + 1)))));
// 
// for (int j=0; j < n; ++j) {
// yaverage[j] = b[j] - arma::mean(b);
// }
// 
//.........这里部分代码省略.........

//[[Rcpp::export(.emplik_intern)]]
List emplik_intern(arma::mat z, arma::colvec mu, arma::vec lam, double eps,
	double M = 1e30, double thresh = 1e-12, int itermax = 1000){
//(arma::mat z, arma::vec mu  = vec(z.n_cols,fill::zeros), double eps = 1/z.nrows, double M = datum::inf);
// # Backtracking line search parameters [Tweak only with extreme caution.]
// # See Boyd and Vandenberghe, pp 464-466.
double ALPHA = 0.3; // seems better than 0.01 on some 2d test data (sometimes fewer iters)
double BETA  = 0.8;
// # We need  0 < ALPHA < 1/2   and 0 < BETA < 1
// # Backtrack threshold: you can miss by this much.
double BACKEPS = 0;
// # Consider replacing 0 by 1e-10 if backtracking seems to be
// # failing due to round off.
int n = z.n_rows;
int d = z.n_cols;
 z.each_row() -= trans(mu);
//  Use lam = 0 or initial lam, whichever is best
arma::vec onen = arma::vec(n,arma::fill::ones);
arma::mat init0 = mllog(onen, eps, M, 2);
arma::mat init(init0.n_rows, init0.n_cols);
if(!any(lam)){
  init = init0;
} else {
  init = mllog(onen+z*lam, eps, M, 2);
	if(sum(init0.col(0) < sum(init.col(0)))){
    lam = arma::rowvec(z.n_cols,arma::fill::zeros);
    init = init0;
  }
}
// # Initial f, g
double fold = sum(init.col(0));
arma::rowvec gold = sum(z.each_col() % init.col(1),0);
bool converged = false;
int iter = 0;
arma::mat oldvals = init;
arma::mat newvals(oldvals.n_rows, oldvals.n_cols);
arma::vec rootllpp;
arma::mat zt(z.n_rows, z.n_cols);
arma::vec yt(z.n_rows);
arma::vec wts;
double fnew; double targ;
double logelr;
arma::vec step(z.n_cols);
int s; double ndec; double gradnorm;
bool backtrack = false;
while(!converged){
  iter += 1;
  //   # Get Newton Step
  rootllpp = sqrt(oldvals.col(2));  //# sqrt 2nd deriv of -llog lik
  zt = z;
     for(int j=0; j<d; j++){
      zt.col(j) %= rootllpp;
    }
    yt   = oldvals.col(1) / rootllpp;
    step = -svdlm(zt,yt);
    backtrack = false;
    s = 1;
    while(!backtrack){
         newvals = mllog(onen + z * (lam+s*step), eps, M, 2);
        fnew = sum(newvals.col(0));
        targ = fold + ALPHA * s * sum(trans(gold) % step) + BACKEPS; // (BACKEPS for roundoff, should not be needed)
         if(fnew <= targ){ // backtracking has converged
          backtrack = true;
          oldvals = newvals;
          fold = fnew;
     			gold = sum(z.each_col() % oldvals.col(1),0);
          lam = lam + s*step;
         } else{
          s = s * BETA;
         }
       }
    //   # Newton decrement and gradient norm
      ndec     = sqrt(sum(square(step % trans(gold))));
      gradnorm = sqrt(sum(square(gold)));
      converged = (ndec * ndec <= thresh);
       if( iter > itermax )
       	break;
     }
     wts = pow(1 + z * lam, -1)/n;
		logelr = sum(mllog(onen + z * lam, eps, M, 0).col(0));
    return Rcpp::List::create(Named("logelr")=logelr, Named("lam") = lam, Named("wts") = wts,
        Named("conv") = converged, Named("niter") = iter,Named("ndec") = ndec, Named("gradnorm") = gradnorm);
     }