-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathwpca.cpp
62 lines (54 loc) · 1.82 KB
/
wpca.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
// This Cpp includes functions to impliment Weighted PCA in
// reference Bai, J., & Liao, Y. (2013). Statistical inferences using large estimated covariances for panel data and factor models. arXiv preprint arXiv:1307.2662.
// and Inferences in panel data with interactive effects using largecovariance matrices
// It considers the heterogeneous error term in approximated factor model.
#include "RcppArmadillo.h"
// [[Rcpp::depends(RcppArmadillo)]]
#include "idrsc2.h"
#define INT_MIN (-INT_MAX - 1)
using namespace Rcpp;
using namespace arma;
using namespace std;
/*
* Auxiliary
*/
//' @keywords internal
//' @noRd
//'
//[[Rcpp::export]]
VECTYPE calculateWeight(const MATTYPE& X, const int& nPCs){
int n = X.n_rows;
MATTYPE U, V;
VECTYPE s;
svd_econ(U, s, V, X);
MATTYPE hX = U.cols(0, nPCs-1) *diagmat(s.subvec(0, nPCs-1)) * trans(V.cols(0,nPCs-1));
ROWVECTYPE Lam_vec = sum((hX-X) % (hX-X))/ n;
return Lam_vec.t();
}
//[[Rcpp::export]]
Rcpp::List wpcaCpp(const MATTYPE&X, const int& nPCs, const bool& weighted=true){
MATTYPE U, V;
VECTYPE s;
MATTYPE PCs, loadings;
svd_econ(U, s, V, X);
PCs = U.cols(0, nPCs-1) *diagmat(s.subvec(0, nPCs-1));
loadings = V.cols(0, nPCs-1);
MATTYPE dX = PCs * loadings.t() - X;
ROWVECTYPE Lam_vec = mean(dX % dX);
if(weighted){
svd_econ(U, s, V, X*diagmat(1.0/ sqrt(Lam_vec)));
// vec s2 = s % s; // s; //
MATTYPE loadings_unscale = diagmat(sqrt(Lam_vec)) * V.cols(0, nPCs-1);
MATTYPE V1;
VECTYPE s1;
svd_econ(loadings, s1, V1, loadings_unscale);
PCs = U.cols(0, nPCs-1) * diagmat(s.subvec(0, nPCs-1)) * V1 * diagmat(s1);
dX = PCs * loadings.t() - X;
Lam_vec = mean(dX % dX);
}
List output = List::create(
Rcpp::Named("PCs") = PCs,
Rcpp::Named("loadings") = loadings,
Rcpp::Named("Lam_vec") = Lam_vec);
return output;
}