mlogit.cpp
Go to the documentation of this file.
1 
9 #include <target/mlogit.hpp>
10 
12 
13 namespace target {
14 
15  using arma::vec;
16  using arma::uvec;
17  using arma::rowvec;
18  using arma::mat;
19 
20  MLogit::MLogit(const arma::uvec &choice,
21  const arma::uvec &alt,
22  const arma::uvec &id_idx,
23  const arma::mat &z1,
24  const arma::mat &z2,
25  const arma::mat &x,
26  unsigned nalt,
27  arma::vec weights) {
28  n = alt.n_elem;
29  J = nalt;
30  ncl = id_idx.n_elem;
31  if (weights.is_empty()) {
32  weights = arma::vec(n);
33  weights.fill(1);
34  }
35  updateData(choice, alt, id_idx, z1, z2, x, weights);
36 
37  if (nalt == 0) {
38  arma::uvec ualt = arma::unique(alt);
39  J = ualt.n_elem;
40  }
41  idx_z1 = arma::uvec(z1.n_cols);
42  unsigned pos = 0;
43  for (unsigned i=0; i < idx_z1.n_elem; i++) {
44  idx_z1(i) = pos;
45  pos++;
46  }
47  idx_z2 = arma::uvec(J*z2.n_cols);
48  for (unsigned i=0; i < idx_z2.n_elem; i++) {
49  idx_z2(i) = pos;
50  pos++;
51  }
52  idx_x = arma::uvec((J-1)*x.n_cols);
53  for (unsigned i=0; i < idx_x.n_elem; i++) {
54  idx_x(i) = pos;
55  pos++;
56  }
57  alt_idx = arma::uvec(J);
58  updateRef(basealt);
59  theta_z1 = arma::vec(z1.n_cols);
60  theta_x = arma::mat(J, x.n_cols);
61  theta_z2 = arma::mat(J, z2.n_cols);
62  theta = arma::vec(idx_z1.n_elem + idx_z2.n_elem + idx_x.n_elem);
63 
64  lp = arma::vec(n);
65  logpr = arma::vec(n);
66  zx = arma::sp_mat(n, theta.n_elem);
67  }
68 
69  double MLogit::loglik() {
70  return arma::as_scalar(sum(logpr%(*Choice())%(*Weights())));
71  }
72 
73  mat MLogit::score(bool update, bool indiv) {
74  if (update) MLogit::updateZX();
75  vec pr = exp(logpr);
76  mat xp = zx;
77  xp.each_col() %= pr;
78  xp = target::groupsum(xp, *cluster(), true);
79  uvec chosen = find((*Choice()));
80  mat score = zx.rows(chosen) - xp;
81  score.each_col() %= (*Weights()).elem(chosen);
82  if (!indiv) score = sum(score, 0); // Column-sums
83  return score;
84  }
85 
86  mat MLogit::hessian(bool update) {
87  if (update) MLogit::updateZX();
88  vec pr = exp(logpr);
89  mat xp = zx;
90  xp.each_col() %= pr;
91  mat tmp = xp;
92  xp = zx - target::groupsum(xp, *cluster(), false);
93  tmp.each_col() %= target::groupsum((*Weights()) %
94  (*Choice()), *cluster(), false);
95  mat hess = -xp.t()*tmp;
96  return hess;
97  }
98 
99  void MLogit::updateRef(unsigned basealt) {
100  // Sets the base/reference alternative
101  this->basealt = basealt;
102  alt_idx.fill(0); //
103  unsigned pos = 0;
104  for (unsigned j=0; j < J; j++) {
105  if (j != basealt) {
106  this->alt_idx[j] = pos;
107  pos++;
108  }
109  }
110  }
111 
112  void MLogit::updatePar(vec theta) {
113  this->theta = theta;
114  // theta = (z1_1,...,z1_n1, z2_1, ..., z2_J*n2, x1_1, ..., x_1_(J-1)*n3)
115  for (unsigned i=0; i < p_z1; i++) {
116  theta_z1(i) = theta(idx_z1[i]);
117  }
118  theta_x.fill(0);
119  theta_z2.fill(0);
120  unsigned pos = 0;
121  if (p_x > 0) {
122  for (unsigned j=0; j < J; j++)
123  if (j != basealt) {
124  for (unsigned i=0; i < p_x; i++) {
125  theta_x(j, i) = theta(idx_x[pos]);
126  pos++;
127  }
128  }
129  }
130  pos = 0;
131  if (p_z2 > 0) {
132  for (unsigned j=0; j < J; j++)
133  for (unsigned i=0; i < p_z2; i++) {
134  theta_z2(j, i) = theta(idx_z2[pos]);
135  pos++;
136  }
137  }
138  }
139 
140  void MLogit::updateZX() {
141  for (unsigned i=0; i < n; i++) {
142  unsigned pos = 0;
143  if (p_z1 > 0) {
144  for (unsigned j=0; j < p_z1; j++)
145  zx(i, j) = (*Z1())(i, j);
146  pos = p_z1;
147  }
148  if (p_z2 > 0) {
149  for (unsigned j=0; j < p_z2; j++)
150  zx(i, j + pos + p_z2*Alt(i)) = (*Z2())(i, j);
151  pos += p_z2*J;
152  }
153  if (p_x > 0) {
154  if (Alt(i) != basealt) {
155  pos += p_x*(alt_idx(Alt(i)));
156  for (unsigned j=0; j < p_x; j++)
157  zx(i, j+pos) = (*X())(i, j);
158  }
159  }
160  }
161  }
162 
163  void MLogit::updateProb() {
164  lp.fill(0); // Linear predictor
165  if (p_z1 > 0) {
166  lp = (*Z1())*theta_z1;
167  }
168  if (p_z2 > 0) {
169  for (unsigned i=0; i < n; i++) {
170  rowvec z2 = Z2(i);
171  lp[i] += as_scalar((*Z2()).row(i)*trans(theta_z2.row(Alt(i))));
172  }
173  }
174  if (p_x > 0) {
175  for (unsigned i=0; i < n; i++) {
176  if (Alt(i) != basealt) {
177  lp[i] += as_scalar((*X()).row(i)*trans(theta_x.row(Alt(i))));
178  }
179  }
180  }
181  // Calculate log-probabilities of choices
182  unsigned stop;
183  unsigned pos = 0;
184  for (unsigned i=0; i < ncl; i++) { // Iterate over individuals
185  unsigned start = this->cluster(i);
186  if (i == (ncl-1)) {
187  stop = n-1;
188  } else {
189  stop = this->cluster(i+1)-1;
190  }
191  vec a = target::softmax(lp.subvec(start, stop)); // Multinomial logit
192  for (unsigned j=0; j < a.n_elem; j++) {
193  logpr(pos) = a(j);
194  pos++;
195  }
196  }
197  }
198 
199 
200 } // namespace target
201 
Multinomial logit models.