cumres.cpp
Go to the documentation of this file.
1 
11 #include <target/cumres.hpp>
12 
13 namespace target {
14 
22  cumres::cumres(const arma::vec &r,
23  const arma::mat &dr,
24  const arma::mat &ic)
25  : r(r), dr(dr), ic(ic) {
26 #ifndef ARMA_R
27  arma::arma_rng::set_seed_random();
28 #endif
29  n = r.n_elem;
30  arma::vec inp(n, 1);
31  for (unsigned i=0; i < n; i++) inp(i) = i;
32  this->ord = arma::conv_to<arma::uvec>::from(inp);
33  this->inp = inp;
34  this->order(inp);
35  this->b = arma::vec();
36  }
37 
38  void cumres::order(const arma::mat &inp, arma::vec b) {
39  unsigned p = inp.n_cols;
40  arma::umat ord(inp.n_rows, p);
41  this->inp = arma::mat(inp.n_rows, p);
42  for (unsigned j=0; j < p; j++) {
43  arma::vec inpj = inp.col(j);
44  ord.col(j) = arma::stable_sort_index(inpj);
45  this->inp.col(j) = inpj.elem(ord.col(j));
46  }
47  this->ord = ord;
48  if (p == 1 && b.n_elem == 0) {
49  eta = arma::cumsum(dr.rows(ord), 0);
50  }
51  if (b.n_elem > 0) {
52  if (b.n_elem < p) {
53  arma::vec newb = arma::vec(p);
54  for (unsigned i=0; i < p; i++) newb(i) = b(0);
55  b = newb;
56  }
57  this->b = b;
58  }
59  }
60 
64  arma::vec cumres::rnorm() {
65 #ifdef ARMA_R
66  Rcpp::RNGScope scope;
67  return Rcpp::as<arma::vec>(Rcpp::rnorm(n));
68 #else
69  return arma::randn<arma::vec>(n);
70 #endif
71  }
72 
73 
89  arma::mat cumres::obs() {
90  arma::mat res(n, inp.n_cols);
91  double scale = std::sqrt(static_cast<double>(n));
92  for (unsigned j=0; j < inp.n_cols; j++) {
93  if (b.n_elem > 0) { // Moving average
94  arma::vec ro = r.elem(ord.col(j));
95  arma::vec rj = arma::cumsum(ro);
96  unsigned pos = 0;
97  double lval = 0.0;
98  for (unsigned i=0; i < ro.n_elem; i++) {
99  if (inp(i)-b(j) < lval)
100  while (inp(pos) >= inp(i)-b(j)) pos++;
101  res(i, j) = (rj(i) - rj(pos)) / scale;
102  }
103  } else { // Standard cumulative sum
104  arma::vec rj = arma::cumsum(r.elem(ord.col(j))) / scale;
105  res.col(j) = rj;
106  }
107  }
108  return res;
109  }
110 
117  arma::mat cumres::sample(const arma::umat &idx) {
118  arma::vec g = rnorm();
119  unsigned N = n;
120  unsigned p = ord.n_cols;
121  if (!idx.is_empty()) {
122  N = idx.n_rows;
123  }
124  arma::mat res(N, p);
125  for (unsigned i=0; i < p; i++) {
126  arma::uvec curord = ord.col(i);
127  arma::vec w1 = arma::cumsum(r.elem(curord)%g);
128  if (p > 1) {
129  // Cumulative derivative of residuals (cumsum for each column):
130  eta = arma::cumsum(dr.rows(curord), 0);
131  }
132  arma::rowvec B(ic.n_cols);
133  for (unsigned j=0; j < ic.n_cols; j++) {
134  arma::vec tmp = ic.col(j);
135  B(j) = sum(tmp.elem(curord)%g);
136  }
137  arma::vec w2(N);
138  if (!idx.is_empty()) {
139  w1 = w1.elem(idx.col(i));
140  for (unsigned j=0; j < N; j++) {
141  w2(j) = arma::as_scalar(B*eta.row(idx(j, i)).t());
142  }
143  } else {
144  for (unsigned j=0; j < n; j++) {
145  w2(j) = arma::as_scalar(B*eta.row(j).t());
146  }
147  }
148  res.col(i) = (w1+w2)/std::sqrt(static_cast<double>(n));
149  }
150 
151  return res;
152  }
153 
162  arma::mat cumres::sample(unsigned R, // Number of samples
163  const arma::umat &idx,
164  bool quantiles) {
165  unsigned p = ord.n_cols;
166  arma::mat res(R, 2*p);
167  // qt.fill(0);
168  // unsigned n = this->n;
169  // arma::mat qt(n, std::ceil(R*0.05));
170  for (unsigned i=0; i < R; i++) {
171  arma::mat wi = this->sample(idx);
172  for (unsigned j=0; j < p; j++) {
173  arma::vec t0 = inp.col(j);
174  if (!idx.is_empty()) {
175  t0 = t0.elem(idx.col(j));
176  }
177  res(i, j*2) = SupTest(wi.col(j));
178  res(i, j*2+1) = L2Test(wi.col(j), t0);
179  }
180  /* if (quantiles) { // TODO: Disable for now. Capture quantiles
181  for (unsigned j=0; j<n; j++) {
182  wi = abs(wi);
183  arma::rowvec qtj = qt.row(j);
184  unsigned i = qtj.index_min();
185  if (wi(j)>qt(j,i)) qt(j,i) = wi(j);
186  }
187  } */
188  }
189  // this->qt = qt;
190  return res;
191  }
192 
193 } // namespace target
arma::mat obs()
Calculate observed cumulative residual process.
Definition: cumres.cpp:89
arma::mat ic
Influence curve.
Definition: cumres.hpp:27
arma::mat sample(const arma::umat &idx=arma::umat())
Simulate one process under the null hypothesis of a correctly specified model.
Definition: cumres.cpp:117
arma::vec rnorm()
Draw n samples from standard normal distribution.
Definition: cumres.cpp:64
arma::mat eta
Cumulative derivative of residuals.
Definition: cumres.hpp:31
Generic class for cumulative residuals.
arma::umat ord
Stores order of observations to cumulate after.
Definition: cumres.hpp:25
arma::mat dr
Derivative of residuals wrt model parameters.
Definition: cumres.hpp:26
arma::mat inp
Variable to order residuals after.
Definition: cumres.hpp:28
arma::vec r
Residuals.
Definition: cumres.hpp:24
arma::vec b
Bandwidth of moving average.
Definition: cumres.hpp:29
cumres(const arma::vec &r, const arma::mat &dr, const arma::mat &ic)
Constructor.
Definition: cumres.cpp:22
unsigned n
Sample size.
Definition: cumres.hpp:23
void order(const arma::mat &inp, arma::vec b=arma::vec())
Set variables to order after and bandwidth.
Definition: cumres.cpp:38