knitr::opts_chunk$set(echo = TRUE)
Notes about code in src/ll_funcs_v7.cpp
This function calculates $\log P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}, \sigma_{1}^2, \sigma_2^2, \rho]$
$$ P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}, \sigma_{1}^2, \sigma_2^2] = (1-q) \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma) \Sigma A(\gamma)^T + S_i(\rho) \right) + \ q \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma + \eta) \Sigma A(\gamma + \eta)^T + S_i(\rho) \right) \ $$
Where $$ A(x) \Sigma A(x)^T + S_i(\rho)= \begin{pmatrix} 1 & 0 \ x & 1 \end{pmatrix} \begin{pmatrix} \sigma_{1}^2 & 0 \ 0 & \sigma_2^2 \end{pmatrix} \begin{pmatrix} 1 & x \ 0 & 1 \end{pmatrix} + S_i(\rho) \ = \begin{pmatrix} 1 & 0 \ x & 1 \end{pmatrix} \begin{pmatrix} \sigma_1^2 & x\sigma_1^2 \ 0& \sigma_2^2 \end{pmatrix} + S_i(\rho) \ = \begin{pmatrix} \sigma_1^2 & x\sigma_1^2 \ x\sigma_1^2 & x^2 \sigma_1^2 + \sigma_2^2 \end{pmatrix} + S_i (\rho)\ = \begin{pmatrix} \sigma_1^2 + s_1^2 & x\sigma_1^2 + \rho s_1 s_2 \ x\sigma_1^2 + \rho s_1 s_2 & x^2 \sigma_1^2 + \sigma_2^2 + s_2^2 \end{pmatrix} $$
Specify the variance matrix for the first part of the equation
double v11_1 = R_pow_di(sigma1,2) + R_pow_di(s1, 2) ; double v12_1 = g*R_pow_di(sigma1,2) + rho*s1*s2 ; double v22_1 = R_pow_di(sigma2,2) + R_pow_di(g*sigma1,2) + R_pow_di(s2, 2) ;
Here g
is $\gamma$. To calulate $N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma) \Sigma A(\gamma)^T + S_i(\rho) \right)$ we use
$$ N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, \begin{pmatrix} v_{11} & v_{12} \ v_{21} & v_{22} \end{pmatrix} \right) = N(\hat{\beta}{1i}; 0, v{11} ) N(\hat{\beta}{2i}; \hat{\beta}{1i}v_{21}/v_{11}, v_{22}-v_{21}v_{12}/v_{11}) $$
Calculate conditional mean and variance for $\hat{\beta}{21}$ (note that $v{12} = v_{21}$ in this case).
double m2_g1_1 = b1*v12_1/v11_1; double v2_g1_1 = v22_1 - R_pow_di(v12_1,2)/v11_1;
Calculate $\log \left( (1-q) \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma) \Sigma A(\gamma)^T + S_i(\rho) \right) \right)$
lik1_ij = log(1-q) + Rf_dnorm4(b1, 0, sqrt(v11_1), 1) + Rf_dnorm4(b2, m2_g1_1, sqrt(v2_g1_1), 1);
Calculate $\log \left( (q) \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma + \eta) \Sigma A(\gamma + \eta)^T + S_i(\rho) \right) \right)$. Here gp
is $\eta$ because $\eta$ was called $\gamma^\prime$ in an earlier version of notes.
double v11_2 = R_pow_di(sigma1,2) + R_pow_di(s1, 2) ; double v12_2 = (g + gp)*R_pow_di(sigma1,2) + rho*s1*s2 ; double v22_2 = R_pow_di(sigma2,2) + R_pow_di((g+gp)*sigma1,2) + R_pow_di(s2, 2) ; double m2_g1_2 = b1*v12_2/v11_2; double v2_g1_2 = v22_2 - R_pow_di(v12_2,2)/v11_2; double lik2_ij = log(q) + Rf_dnorm4(b1, 0, sqrt(v11_2), 1) + Rf_dnorm4(b2, m2_g1_2, sqrt(v2_g1_2), 1);
Add them together
double lik_ij = Rf_logspace_add(lik1_ij, lik2_ij); return lik_ij;
This function integrates over the distribution of $\sigma_1$ and $\sigma_2$ so we calculate
$$ P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}] = \sum_{k=1}^{K} \pi_k P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}, \sigma_{1k}^2, \sigma_{2k}^2] $$
int k = sigma1.length(); std::vector<double> lik_i(k); for(size_t j = 0; j < k; j ++){ lik_i[j] = log(pi[j]) + ll_v7_ij(rho, g, gp, q, sigma1[j], sigma2[j], b1, b2, s1, s2); } double result = Rf_logspace_sum(&lik_i[0], k); return result;
Adds up over variants
int p = beta_hat_2.size(); NumericVector tll(p); RVector<double> ll(tll); parallel_for(blocked_range<size_t>(0, p), [&](const blocked_range<size_t>& r){ for(size_t i=r.begin(); i!=r.end(); ++i){ ll[i] = ll_v7_i(rho, g, gp, q, sigma1,sigma2,pi, beta_hat_1[i], beta_hat_2[i], seb1[i], seb2[i]); }}); double result = sum(tll); return result;
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.