#----------------------
# Install the package
#----------------------

library("NLPmix")

#-----------------------------------------------------------
# Missepecified Model. Paper: On choosing mixture components 
# via Non-local priors
#-----------------------------------------------------------

set.seed(123456)
library(mvtnorm)
n.sample.miss<-600
sigma.s1<-matrix(c(1,-0.5,-0.5,1),2,2)
sigma.s2<-matrix(c(1,-0.5,-0.5,1),2,2)
sigma.s3<-matrix(c(1,-0.5,-0.5,1),2,2)
mu1.t<-c(-1,1)
mu2.t<-c(1,-1)
mu3.t<-c(6,6)
y.s<-rbind(rmvt(n.sample.miss/3,delta=mu1.t,sigma.s1, df = 4),
rmvt(n.sample.miss/3,delta=mu2.t,sigma.s2, df = 4),
rmvt(n.sample.miss/3,delta=mu3.t,sigma.s3, df = 4)
)
y.s2<-as.data.frame(y.s)
plot(y.s2)

#----------------------------------------
# Applying the EM algorithm under NLPs
#----------------------------------------

#--------------------------
# Number of components
#--------------------------

number.c.NLPs.cov<-3

#----------------------------------------------
# Initial values
#----------------------------------------------

library(mclust)
bic.cov2 <- mclustBIC(y.s,G=number.c.NLPs.cov,modelNames = c("EEE")) 
mclust.cov2 <- mclustModel(y.s,BICvalues=bic.cov2 ) 
pars.cov <- mclust.cov2$parameters
mixEMNLPs.cov=list(k=number.c.NLPs.cov,Mu0 = t(pars.cov$mean),
Sigma = pars.cov$variance$Sigma,
Ps = pars.cov$pro
, g=1.87, m=matrix(0,1,dim(y.s)[2]),
nu=5+(dim(y.s)[2]-1), S=diag(dim(y.s)[2])*((5+(dim(y.s)[2]-1))), alfa=16.5, penalize=4)

#----------------------------------------------
# Running the function
#----------------------------------------------

EM.results2.cov<-emNLP(mixEMNLPs.cov,y.s,0.001, equalcol=TRUE)

#--------------------------------------------------------------
# Plot using the EM package
#--------------------------------------------------------------

parametersBIC2.cov<-mclust.cov2$parameters
parametersBIC2.cov$pro<-EM.results2.cov$P
parametersBIC2.cov$mean<-t(EM.results2.cov$Mu)
parametersBIC2.cov$variance$Sigma<-EM.results2.cov$Sigma
parametersBIC2.cov$variance$sigma[,,1]<-EM.results2.cov$Sigma
parametersBIC2.cov$variance$sigma[,,2]<-EM.results2.cov$Sigma
parametersBIC2.cov$variance$sigma[,,3]<-EM.results2.cov$Sigma
parametersBIC2.cov$variance$cholSigma<-chol(EM.results2.cov$Sigma)
parametersBIC2.cov$z<-EM.results2.cov$z
parametersBIC2miss.cov2<-parametersBIC2.cov
mclust2Dplot(data = y.s2, what = "classification",parameters = parametersBIC2.cov, 
z = parametersBIC2miss.cov2$z,col=8, main= "NLPs")





#------------------------------------------------------------
# 2) Computing the Integrated Likelihood
#------------------------------------------------------------

#------------------------------------------------------
# MCMC: number of iterations for the MCMC
# burn: burning period
# Precision parameter g11 (for LPs) and g21 (for NLPs)
#-------------------------------------------------------

MCMC<-10
burn<-5
g11<-2.97
g21<-1.87


#------------------------------------------------------------
# Case: a common variance-covariance matrix for LPs for 
# the missespecified model with a 5 component mixture
#------------------------------------------------------------

#---------------
# Initial values
#---------------

library(mclust)
bic5 <- mclustBIC(y.s,G=5,modelNames = c("EEE")) 
mclust5 <- mclustModel(y.s,BICvalues=bic5) 
pars5 <- mclust5$parameters

#-----------------------------
# Integrated likelihood
#-----------------------------

mixMCMC5LP5=list(k=5,Mu0 = t(pars5$mean),
Sigma0 = pars5$variance$Sigma,
Ps0 = pars5$pro
, g=g11 , m=matrix(0,1,dim(y.s)[2]),
nu=5+(dim(y.s)[2]-1), S=diag(dim(y.s)[2])*((5+(dim(y.s)[2]-1))), alfa=1, penalize=1)
simu_LP5<-Integratemix(mixMCMC5LP5,equalcol=TRUE,y.s,MCMC,burn)
LP5<-simu_LP5$local.integrate
LP5

#------------------------------------------------------------
# Case: a common variance-covariance matrix for NLPs for 
# the missespecified model with a 3 component mixture
#------------------------------------------------------------ 

#---------------
# Initial values
#---------------

bic3 <- mclustBIC(y.s,G=3,modelNames = c("EEE")) 
mclust3 <- mclustModel(y.s,BICvalues=bic3) 
pars3 <- mclust3$parameters

#-----------------------------
# Integrated likelihood
#-----------------------------

mixMCMC3NLP3=list(k=3,Mu0 = t(pars3$mean),
Sigma0 = pars3$variance$Sigma,
Ps0 = pars3$pro
, g=g21 , m=matrix(0,1,dim(y.s)[2]),
nu=5+(dim(y.s)[2]-1), S=diag(dim(y.s)[2])*((5+(dim(y.s)[2]-1))), alfa=16.5, penalize=96)
simu_NLP3<-Integratemix(mixMCMC3NLP3,equalcol=TRUE,y.s,MCMC,burn)
NLP3<-simu_NLP3$non.local.integrate
NLP3