1
#include <RcppArmadillo.h>
2
#include <rgen.h>
3
#include <simcdm.h>
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<arma::mat>(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<arma::vec>(C);
46
    double etaij;
47

48 1
    for (unsigned int i = 0; i < N; i++) {
49 1
        PYCS = arma::ones<arma::vec>(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<double>::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<arma::vec>(N);
123 1
        qq = arma::conv_to<double>::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<double>::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<unsigned int>(pow(2.0, static_cast<double>(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<arma::vec>(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<arma::vec>(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<arma::mat>(N,K); //K>1 is assumed
213 1
    arma::mat alphas = arma::randu<arma::mat>(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<arma::vec>(J);
218 1
    arma::vec gs = arma::randu<arma::vec>(J);
219 1
    arma::vec pis = arma::randu<arma::vec>(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<arma::vec>(step1a[2]);
231 1
        CLASSES.col(t) = Rcpp::as<arma::vec>(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<arma::vec>(step1b[0]);
239 1
        gs = Rcpp::as<arma::vec>(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 .

Loading