nb.cpp
1 
10 #include <cmath>
11 #include <target/nb.hpp>
12 
13 namespace target {
14 
15  std::vector<raggedArray> nb(arma::vec y,
16  arma::mat x,
17  arma::uvec xlev,
18  arma::vec ylev,
19  arma::vec weights,
20  double laplacesmooth) {
21  if (ylev.n_elem == 0) {
22  ylev = arma::unique(y);
23  }
24  unsigned K = ylev.n_elem;
25  unsigned n = y.n_elem;
26  unsigned p = x.n_cols;
27  if (xlev.n_elem < p) {
28  xlev = arma::uvec(p);
29  for (unsigned j=0; j < p; j++) xlev(j) = j;
30  }
31  if (weights.n_elem < n) {
32  weights = arma::vec(n);
33  for (unsigned i=0; i < n; i++) weights(i) = 1;
34  }
35 
36  std::vector<raggedArray> res(K);
37  for (unsigned k=0; k < K; k++) {
38  unsigned N = 0;
39  for (unsigned i=0; i < n; i++) {
40  if (y[i] == ylev[k]) N++;
41  }
42  arma::uvec idx(N);
43  N = 0;
44  for (unsigned i=0; i < n; i++) {
45  if (y[i] == ylev[k]) {
46  idx(N) = i;
47  N++;
48  }
49  }
50  res[k] = pcond(idx, x, xlev, weights, laplacesmooth);
51  }
52  return res;
53  }
54 
55 
56  raggedArray pcond(const arma::uvec &idx,
57  const arma::mat &x,
58  const arma::uvec &xlev,
59  const arma::vec &weights,
60  double laplacesmooth) {
61  unsigned p = x.n_cols;
62  raggedArray val(p);
63  double sw = 0;
64  for (unsigned i=0; i < idx.n_elem; i++) sw += weights(i);
65 
66  for (unsigned j=0; j < p; j++) {
67  if (xlev(j) > 0) { // factor
68  arma::vec est(xlev[j]);
69  for (unsigned i=0; i < idx.n_elem; i++) {
70  double xval = x(idx(i), j);
71  est(static_cast<int>(xval)) += weights(idx(i));
72  }
73  double s = 0;
74  if (laplacesmooth > 0) {
75  for (unsigned i=0; i < est.n_elem; i++) est(i) += laplacesmooth;
76  }
77  for (unsigned i=0; i < est.n_elem; i++) s += est(i);
78  for (unsigned i=0; i < est.n_elem; i++) est(i) /= s;
79  val[j] = est;
80  } else {
81  double m = 0, ss = 0;
82  arma::vec est(2);
83  for (unsigned i=0; i < idx.n_elem; i++) {
84  double xval = x(idx(i), j);
85  double w = weights(idx(i))/sw;
86  double xw = xval*w;
87  m += xw;
88  ss += xval*xw;
89  }
90  if (idx.n_elem == 1) {
91  est(1) = 0;
92  } else {
93  est(1) = std::sqrt((ss-m*m)*sw/(sw-1));
94  }
95  est(0) = m;
96  val[j] = est;
97  }
98  }
99  return val;
100  }
101 
102 
103  arma::mat prednb(const arma::mat &X,
104  const raggedArray &condprob,
105  const raggedArray &xord,
106  arma::uvec multinomial,
107  arma::vec prior,
108  double threshold) {
109  return X;
110  }
111 
112 // unsigned p = X.n_cols;
113 // unsigned n = X.n_rows;
114 // unsigned nclasses = condprob.size()/p;
115 // arma::mat lposterior(n,nclasses);
116 // arma::vec px(n); px.fill(0);
117 
118 // for (unsigned k=0; k<nclasses; k++) {
119 // arma::vec lpcond(n); lpcond.fill(0);
120 // for (unsigned j=0; j<p; j++) {
121 // unsigned pos = k*p + j;
122 // NumericVector estx = condprob[pos];
123 // if (multinomial[j]) {
124 // arma::uvec x = arma::conv_to<arma::uvec>::from(X.col(j));
125 // NumericVector xx = xord(j);
126 // NumericVector est = clone(xx);
127 // LogicalVector estna = is_na(est);
128 // for (unsigned i=0; i<est.size(); i++)
129 // if (estna[i]) {
130 // est[i] = threshold;
131 // } else {
132 // est[i] = estx[est[i]];
133 // if (est[i]<threshold) est[i] = threshold;
134 // }
135 // double s = sum(est);
136 // arma::vec logest(est.size());
137 // for (unsigned i=0; i<est.size(); i++) {
138 // logest[i] = log(est[i]/s);
139 // }
140 // lpcond += logest.elem(x);
141 // } else { // Gaussian
142 // arma::vec x = X.col(j);
143 // LogicalVector estna = is_na(estx);
144 // if (estna[0]) estx[0] = 0;
145 // if (estna[1] | (estx[1]<1E-16)) estx[1] = 1;
146 // for (unsigned i=0; i<n; i++) {
147 // lpcond[i] += Rf_dnorm4(x[i],estx[0],estx[1],true);
148 // }
149 // }
150 // }
151 // lposterior.col(k) = lpcond + std::log(prior[k]); // log P(x,c)
152 // px = px + exp(lposterior.col(k)); // P(x)
153 // }
154 // arma::colvec lpx = -log(px);
155 // lposterior.each_col() += lpx; // log P(c|x)
156 // return lposterior;
157 // }
158 
159 
160 } // namespace target
Weighted Naive Bayes.