forked from WinVector/RcppDynProg
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathxlin_pfits.cpp
More file actions
79 lines (73 loc) · 2.28 KB
/
xlin_pfits.cpp
File metadata and controls
79 lines (73 loc) · 2.28 KB
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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include <RcppArmadillo.h>
#include <math.h>
using Rcpp::NumericVector;
using Rcpp::NumericMatrix;
using Rcpp::IntegerVector;
//' xlin_pfits
//'
//' Calculate out of sample linear fit predictions using pseudo-inverse.
//' Please see: \url{https://win-vector.com/2019/01/08/a-beautiful-2-by-2-matrix-identity/}.
//' Zero indexed.
//'
//' @param x NumericVector, explanatory variable (length>=2).
//' @param y NumericVector, values to fit.
//' @param w NumericVector, weights (positive).
//' @param i integer, first index (inclusive).
//' @param j integer, j>=i+2 last index (inclusive);
//' @return vector of predictions.
//'
//' @keywords internal
//'
//' @examples
//'
//' xlin_pfits(c(1, 2, 3, 4), c(1, 2, 2, 1), c(1, 1, 1, 1), 0, 3)
//'
//' @export
// [[Rcpp::export]]
NumericVector xlin_pfits(NumericVector x, NumericVector y, NumericVector w,
const int i, const int j) {
const int vlen = x.length();
if((i<0) || (j>=vlen) || (vlen!=y.length()) || (vlen!=w.length())) {
throw std::range_error("Inadmissible value");
}
// build fitting data
double xx_0_0 = 0;
double xx_1_0 = 0;
double xx_0_1 = 0;
double xx_1_1 = 0;
double sy_0 = 0;
double xy_1 = 0;
double sum_w = 0;
for(int k=i; k<=j; ++k) {
xx_0_0 = xx_0_0 + w(k)*1;
xx_1_0 = xx_1_0 + w(k)*x(k);
xx_0_1 = xx_0_1 + w(k)*x(k);
xx_1_1 = xx_1_1 + w(k)*x(k)*x(k);
sy_0 = sy_0 + w(k)*y(k);
xy_1 = xy_1 + w(k)*x(k)*y(k);
sum_w = sum_w + w(k);
}
NumericVector fits = NumericVector(1+j-i);
for(int k=i; k<=j; ++k) {
// pull out k'th observation
const double xxk_0_0 = xx_0_0 - w(k)*1;
const double xxk_1_0 = xx_1_0 - w(k)*x(k);
const double xxk_0_1 = xx_0_1 - w(k)*x(k);
const double xxk_1_1 = xx_1_1 - w(k)*x(k)*x(k);
const double syk_0 = sy_0 - w(k)*y(k);
const double xyk_1 = xy_1 - w(k)*x(k)*y(k);
const double det = xxk_0_0*xxk_1_1 - xxk_0_1*xxk_1_0;
if(det!=0.0) {
// solve linear system and form estimate
const double c0 = (xxk_1_1*syk_0 - xxk_0_1*xyk_1)/det;
const double c1 = (-xxk_1_0*syk_0 + xxk_0_0*xyk_1)/det;
// form estimate
const double y_est = c0 + c1*x(k);
fits(k-i) = y_est;
} else {
// estimate directly from hold-out mean
fits(k-i) = syk_0/(sum_w - w(k));
}
}
return fits;
}