Last updated: 2018-08-03

\[ \left(\begin{matrix} \hat{x} \\ \hat{y} \end{matrix} \right) | \left(\begin{matrix} x \\ y \end{matrix} \right) \sim N(\left(\begin{matrix} \hat{x} \\ \hat{y} \end{matrix} \right) ; \left(\begin{matrix} x \\ y \end{matrix} \right), \left( \begin{matrix} 1 & \rho \\ \rho & 1 \end{matrix} \right)) \] \[ \left(\begin{matrix} x \\ y \end{matrix} \right) \sim \sum_{k=0}^{K} \pi_{k} N( \left(\begin{matrix} x \\ y \end{matrix} \right); 0, U_{k} ) \] \(\Rightarrow\) \[ \left(\begin{matrix} \hat{x} \\ \hat{y} \end{matrix} \right) \sim \sum_{k=0}^{K} \pi_{k} N( \left(\begin{matrix} \hat{x} \\ \hat{y} \end{matrix} \right); 0, \left( \begin{matrix} 1 & \rho \\ \rho & 1 \end{matrix} \right) + U_{k} ) \] \[ \Sigma_{k} = \left( \begin{matrix} 1 & \rho \\ \rho & 1 \end{matrix} \right) + U_{k} = \left( \begin{matrix} 1 & \rho \\ \rho & 1 \end{matrix} \right) + \left( \begin{matrix} u_{k11} & u_{k12} \\ u_{k21} & u_{k22} \end{matrix} \right) = \left( \begin{matrix} 1+u_{k11} & \rho+u_{k12} \\ \rho+u_{k21} & 1+u_{k22} \end{matrix} \right) \] Let \(\sigma_{k11} = \sqrt{1+u_{k11}}\), \(\sigma_{k22} = \sqrt{1+u_{k22}}\), \(\phi_{k}=\frac{\rho+u_{k12}}{\sigma_{k11}\sigma_{k22}}\)

library(mashr)
Loading required package: ashr
set.seed(1)
n <- 10000; pi <- c(0.8,0.2); rho <- 0.8
V <- matrix(c(1,rho,rho,1),2,2)
U0 <- matrix(0,2,2)
U1 <- diag(2)
X0 <- mvtnorm::rmvnorm(n*pi[1], sigma = V + U0)
X1 <- mvtnorm::rmvnorm(n*pi[2], sigma = V + U1)
X <- rbind(X0,X1)

U2 <- matrix(1,2,2)
U3 <- matrix(0,2,2); U3[1,1] <- 1
U4 <- matrix(0,2,2); U4[2,2] <- 1
Ulist <- list(U0 = U0, U1 = U1, U2=U2, U3 = U3, U4=U4)

MLE

The loglikelihood is (with penalty) \[ l(\rho, \pi) = \sum_{i=1}^{n} \log \sum_{k=0}^{K} \pi_{k}N(x_{i}; 0, \Sigma_{k}) + \sum_{k=0}^{K} (\lambda_{k}-1) \log \pi_{k} \]

The penalty on \(\pi\) encourages over-estimation of \(\pi_{0}\), \(\lambda_{k}\geq 1\).

\[ l(\rho, \pi) = \sum_{i=1}^{n} \log \sum_{k=0}^{K} \pi_{k}\frac{1}{2\pi\sigma_{k11}\sigma_{k22}\sqrt{1-\phi_{k}^2}} \exp\left( -\frac{1}{2(1-\phi_{k}^2)}\left[ \frac{x_{i}^2}{\sigma_{k11}^2} + \frac{y_{i}^2}{\sigma_{k22}^2} - \frac{2\phi_{k}x_{i}y_{i}}{\sigma_{k11}\sigma_{k22}}\right] \right) + \sum_{k=0}^{K} (\lambda_{k}-1) \log \pi_{k} \]

Note: This probelm is convex with respect to \(\pi\). In terms of \(\rho\), the covenxity depends on the data.

Algorithm:

Input: X, init_pi, init_rho, Ulist
Compute loglikelihood
delta = 1
while delta > tol
  Given rho, estimate pi by max loglikelihood
  Given pi, estimate rho by max loglikelihood
  Compute loglikelihood
  Update delta
source('../code/estimate_cor.R')
result <- optimize_pi_rho_times(X, Ulist)
plot(result[[1]]$loglik)

The estimated \(\rho\) is 0.7998163.

Session information

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] mashr_0.2-9 ashr_2.2-10

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.18      knitr_1.20        magrittr_1.5     
 [4] REBayes_1.3       MASS_7.3-50       doParallel_1.0.11
 [7] pscl_1.5.2        SQUAREM_2017.10-1 lattice_0.20-35  
[10] foreach_1.4.4     plyr_1.8.4        stringr_1.3.1    
[13] tools_3.5.1       parallel_3.5.1    grid_3.5.1       
[16] rmeta_3.0         htmltools_0.3.6   iterators_1.0.10 
[19] assertthat_0.2.0  yaml_2.2.0        rprojroot_1.3-2  
[22] digest_0.6.15     Matrix_1.2-14     codetools_0.2-15 
[25] evaluate_0.11     rmarkdown_1.10    stringi_1.2.4    
[28] compiler_3.5.1    Rmosek_8.0.69     backports_1.1.2  
[31] mvtnorm_1.0-8     truncnorm_1.0-8

This R Markdown site was created with workflowr