pava.cpp
Go to the documentation of this file.
1 
9 #include <cmath>
10 #include <target/utils.hpp>
11 
12 
13 namespace target {
14 
22  arma::mat pava(arma::vec y,
23  const arma::vec &x,
24  arma::vec w) {
25  unsigned n = y.n_elem;
26  unsigned nb = n; // Number of blocks
27  if (!x.is_empty()) {
28  if (x.n_elem != n)
29  throw std::range_error("Wrong length of predictor variable 'x'");
30  }
31  if (w.is_empty()) {
32  w.resize(n);
33  for (unsigned i=0; i < n; i++) w[i]=1;
34  } else {
35  if (w.n_elem != n)
36  throw std::range_error("Wrong length of weights variable 'weights'");
37  }
38  std::vector<unsigned> poolEnd(n);
39  // Initialize with n pools (each observation defines a block)
40  for (unsigned i=0; i < n; i++) poolEnd[i] = i;
41 
42  unsigned i1, i2;
43  double w0;
44  while (true) {
45  bool stable = true;
46  unsigned i = 0;
47  unsigned nviolators = 0;
48  while (i < (nb-1)) {
49  unsigned pos = i + nviolators;
50  poolEnd[i] = poolEnd[pos];
51  poolEnd[i+1] = poolEnd[pos+1];
52  i1 = poolEnd[i];
53  i2 = poolEnd[i+1];
54  if (y[i1] >= y[i2]) { // Violator => pool new os. with current block
55  stable = false;
56  w0 = w[i1]+w[i2];
57  y[i2] = (w[i1]*y[i1]+w[i2]*y[i2])/w0;
58  w[i2] = w0;
59  poolEnd[i] = poolEnd[i+1];
60  nviolators++;
61  nb--;
62  }
63  i++;
64  }
65  poolEnd[nb-1] = n-1;
66  if (stable) break;
67  }
68 
69  arma::mat res(nb, 2); // 1. column: value, 2. column: index
70  for (unsigned i=0; i < nb; i++) {
71  res(i, 0) = y[poolEnd[i]];
72  }
73  res(0, 1) = 0; // index
74  for (unsigned i=0; i < (nb-1); i++) {
75  res(i+1, 1) = poolEnd[i]+1; // Beginning of each pool right after
76  // previous pool ends
77  }
78  return res;
79  }
80 
81 } // namespace target
Various utility functions and constants.
arma::mat pava(arma::vec y, const arma::vec &x=arma::vec(), arma::vec w=arma::vec())
Weighted Pooled Adjacent Violator Algorithm.
Definition: pava.cpp:22