1 #include  2 #include  3 #include  4 5 //' Update attributes and latent class probabilities  6 //'  7 //' Update attributes and latent class probabilities by sampling from full  8 //' conditional distribution.  9 //'  10 //' @param Amat A \eqn{C \times K}{C x K} matrix of latent classes.  11 //' @param Q A \eqn{N \times K}{N x K} matrix indicating which  12 //' skills are required for which items.  13 //' @param ss A \eqn{J} vector of item slipping parameters.  14 //' @param gs A \eqn{J} vector of item guessing parameters.  15 //' @param Y A \eqn{N \times J}{N x J} matrix of observed responses.  16 //' @param PIs A \eqn{C} vector of latent class probabilities.  17 //' @param ALPHAS A \eqn{N \times K}{N x K} matrix of latent attributes.  18 //' @param delta0 A \eqn{J} vector of Dirichlet prior parameters.  19 //'  20 //' @return  21 //' A \eqn{N \times K}{N x K} matrix of attributes and a C vector of  22 //' class probabilities.  23 //'  24 //' @author  25 //' Steven Andrew Culpepper  26 //'  27 //' @noRd  28 // [[Rcpp::export]]  29 1 Rcpp::List update_alpha(const arma::mat &Amat, const arma::mat &Q,  30  const arma::vec &ss, const arma::vec &gs,  31  const arma::mat &Y, const arma::vec &PIs,  32  arma::mat &ALPHAS, const arma::vec &delta0)  33 {  34 35 1  unsigned int N = Y.n_rows;  36 1  unsigned int J = Q.n_rows;  37 1  unsigned int C = Amat.n_rows;  38 1  unsigned int K = Q.n_cols;  39  double ci;  40 41 1  arma::mat alphas_new = arma::zeros(N, K);  42 1  arma::vec PYCS(C);  43 1  arma::vec CLASSES(N);  44 1  arma::vec PS;  45 1  arma::vec Ncs = arma::zeros(C);  46  double etaij;  47 48 1  for (unsigned int i = 0; i < N; i++) {  49 1  PYCS = arma::ones(C);  50 51 1  for (unsigned int c = 0; c < C; c++) {  52 53 1  for (unsigned int j = 0; j < J; j++) {  54 1  etaij = 1.0;  55 1  if (arma::dot(Amat.row(c), Q.row(j)) <  56 1  arma::dot(Q.row(j), Q.row(j))) {  57 1  etaij = 0.0;  58  }  59 1  if (etaij == 1.0 and Y(i, j) == 1.0) {  60 1  PYCS(c) = (1.0 - ss(j)) * PYCS(c);  61  }  62 1  if (etaij == 0.0 and Y(i, j) == 1.0) {  63 1  PYCS(c) = gs(j) * PYCS(c);  64  }  65 1  if (etaij == 1.0 and Y(i, j) == 0.0) {  66 1  PYCS(c) = ss(j) * PYCS(c);  67  } else {  68 1  PYCS(c) = (1.0 - gs(j)) * PYCS(c);  69  }  70  }  71  }  72 1  PS = PYCS % PIs / (arma::conv_to::from(PYCS.t() * PIs));  73 1  ci = rgen::rmultinomial(PS);  74 1  ALPHAS.row(i) = Amat.row(ci);  75 1  Ncs(ci) = 1.0 + Ncs(ci);  76 1  CLASSES(i) = ci;  77  }  78 1  return Rcpp::List::create(  79 1  Rcpp::Named("PYCS") = PYCS, Rcpp::Named("PS") = PS,  80 1  Rcpp::Named("PIs_new") = rgen::rdirichlet(Ncs + delta0),  81 1  Rcpp::Named("CLASSES") = CLASSES);  82 }  83 84 //' Update item parameters  85 //'  86 //' Update guessing and slipping parameters from full conditional distribution.  87 //'  88 //' @param Y A N by J matrix of observed responses.  89 //' @param Q A N by K matrix indicating which skills are required  90 //' for which items.  91 //' @param ALPHAS A N by K matrix of latent attributes.  92 //' @param ss_old A J vector of item slipping parameters from prior iteration.  93 //' @param as0 Slipping prior alpha parameter for Beta distribution.  94 //' @param bs0 Slipping prior beta parameter for Beta distribution.  95 //' @param ag0 Guessing prior alpha parameter for Beta distribution.  96 //' @param bg0 Guessing prior beta parameter for Beta distribution.  97 //'  98 //' @return  99 //' A list with two J vectors of guessing and slipping parameters.  100 //'  101 //' @author Steven Andrew Culpepper  102 //' @noRd  103 // [[Rcpp::export]]  104 1 Rcpp::List update_sg(const arma::mat &Y, const arma::mat &Q,  105  const arma::mat &ALPHAS, const arma::vec &ss_old,  106  double as0, double bs0, double ag0, double bg0)  107 {  108 109 1  unsigned int N = Y.n_rows;  110 1  unsigned int J = Y.n_cols;  111 112 1  arma::vec ETA;  113 1  arma::vec ss_new(J);  114 1  arma::vec gs_new(J);  115 1  arma::mat AQ = ALPHAS * Q.t();  116  double T, S, G, y_dot_eta, qq, ps, pg;  117  double ug, us;  118 119 1  for (unsigned int j = 0; j < J; j++) {  120 1  us = R::runif(0, 1);  121 1  ug = R::runif(0, 1);  122 1  ETA = arma::zeros(N);  123 1  qq = arma::conv_to::from(Q.row(j) * (Q.row(j)).t());  124 1  ETA.elem(arma::find(AQ.col(j) == qq)).fill(1.0);  125 126 1  y_dot_eta = arma::conv_to::from((Y.col(j)).t() * ETA);  127 1  T = sum(ETA);  128 1  S = T - y_dot_eta;  129 1  G = sum(Y.col(j)) - y_dot_eta;  130 131  // sample s and g as linearly truncated bivariate beta  132 133  // draw g conditoned upon s_t-1  134 1  pg = R::pbeta(1.0 - ss_old(j), G + ag0, N - T - G + bg0, 1, 0);  135 1  gs_new(j) = R::qbeta(ug * pg, G + ag0, N - T - G + bg0, 1, 0);  136  // draw s conditoned upon g  137 1  ps = R::pbeta(1.0 - gs_new(j), S + as0, T - S + bs0, 1, 0);  138 1  ss_new(j) = R::qbeta(us * ps, S + as0, T - S + bs0, 1, 0);  139  }  140 1  return Rcpp::List::create(Rcpp::Named("ss_new") = ss_new,  141 1  Rcpp::Named("gs_new") = gs_new);  142 }  143 144 //' Generate Posterior Distribution with Gibbs sampler  145 //'  146 //' Function for sampling parameters from full conditional distributions.  147 //' The function returns a list of arrays or matrices with parameter posterior  148 //' samples. Note that the output includes the posterior samples in objects.  149 //'  150 //' @param Y A \eqn{N \times J}{N x J} matrix of observed  151 //' responses.  152 //' @param Amat A \eqn{C \times K}{C x K} matrix of latent  153 //' classes.  154 //' @param Q A \eqn{N \times K}{N x K} matrix indicating  155 //' which skills are required for which items.  156 //' @param chain_length Number of MCMC iterations.  157 //'  158 //' @return  159 //' A list with samples from the posterior distribution with each  160 //' entry named:  161 //'  162 //' - CLASSES = individual attribute profiles,  163 //' - PIs = latent class proportions,  164 //' - SigS = item slipping parameters, and  165 //' - GamS = item guessing parameters.  166 //'  167 //' @author  168 //' Steven Andrew Culpepper  169 //'  170 //' @seealso  171 //' [simcdm::sim_dina_items()] and [simcdm::attribute_classes()]  172 //'  173 //' @noRd  174 // [[Rcpp::export]]  175 1 Rcpp::List DINA_Gibbs_cpp(const arma::mat &Y, const arma::mat &Q,  176  unsigned int chain_length = 10000)  177 {  178 179  // Number of Observations  180 1  unsigned int N = Y.n_rows;  181 182  // Number of Items  183 1  unsigned int J = Y.n_cols;  184 185  // Number of Attributes  186 1  unsigned int K = Q.n_cols;  187 188  // Number of Latent Classes (2^k)  189 1  unsigned int C = static_cast(pow(2.0, static_cast(K)));  190 191  // Generate the latent class alpha matrix  192 1  arma::mat Amat = simcdm::attribute_classes(K);  193 194  // Prior values for betas and Dirichlet distribution  195 1  arma::vec delta0 = arma::ones(C);  196 1  double as0 = 1.0;  197 1  double bs0 = 1.0;  198 1  double ag0 = 1.0;  199 1  double bg0 = 1.0;  200 201 1  arma::vec pil0 = arma::ones(C) / double(C); // prior probability  202 203  // Saving output  204 1  arma::mat SigS(J, chain_length);  205 1  arma::mat GamS(J, chain_length);  206 1  arma::mat US(J, chain_length);  207 1  arma::mat PIs(C, chain_length);  208 1  arma::mat CLASSES(N, chain_length);  209 1  arma::cube QS(J, K, chain_length);  210 211  // Need to initialize, alphas, ss, gs, and pis  212  // arma::mat alphas = arma::zeros(N,K); //K>1 is assumed  213 1  arma::mat alphas = arma::randu(N, K); // K>1 is assumed  214 1  alphas.elem(find(alphas > 0.5)).ones();  215 1  alphas.elem(find(alphas <= 0.5)).zeros();  216 217 1  arma::vec ss = arma::randu(J);  218 1  arma::vec gs = arma::randu(J);  219 1  arma::vec pis = arma::randu(C);  220 221  // Start Markov chain  222 1  for (unsigned int t = 0; t < chain_length; t++) {  223 224  // updata alpha and pi  225  Rcpp::List step1a =  226 1  update_alpha(Amat, Q, ss, gs, Y, pis, alphas, delta0);  227 228  // update value for pis. alphas are updated via pointer. save classes  229  // and PIs  230 1  pis = Rcpp::as(step1a[2]);  231 1  CLASSES.col(t) = Rcpp::as(step1a[3]);  232 1  PIs.col(t) = pis;  233 234  // update s and g  235 1  Rcpp::List step1b = update_sg(Y, Q, alphas, ss, as0, bs0, ag0, bg0);  236 237  // update value for ss and gs.  238 1  ss = Rcpp::as(step1b[0]);  239 1  gs = Rcpp::as(step1b[1]);  240 1  SigS.col(t) = ss;  241 1  GamS.col(t) = gs;  242  }  243 244 1  return Rcpp::List::create(  245 1  Rcpp::Named("CLASSES", CLASSES), Rcpp::Named("PIs", PIs),  246 1  Rcpp::Named("SigS", SigS), Rcpp::Named("GamS", GamS));  247 } 

Read our documentation on viewing source code .