### R functions ### rm(list=ls()) em_pcd_exchangeable <- function(x, n, g, p, class, tol, restriction, constrain, iteration){ pi<-double(g) alpha<-double(1) mu <- double(p*g) sigma <- double(p*g) loglik <- double(1) convergence<-integer(1) ez<-double(g^p*n) j<-integer(1) id<-integer(p) pi <- rep(0.7/g, g) alpha <- 0.3/(g^p-g) mu<-rep(c(-2,0,2), p) sigma<-rep(c(1.5,1,1.5),p) id<-rep(0,p) j<-1 results<-.C("em_pdim_pcd_exchange", as.integer(n), as.integer(g), as.integer(p),as.integer(j), as.integer(id), as.double(ez), as.double(pi),as.double(alpha), as.double(x), as.double(mu), as.double(sigma), loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) pi<-results[[7]] alpha<-results[[8]] mu<-results[[10]] sigma<-results[[11]] loglik<-results[[12]] tol<-results[[13]] restriction<-results[[14]] constrain<-results[[15]] iteration<-results[[16]] convergence<-results[[17]] par.save.temp<-c(pi,alpha,mu,sigma,loglik,tol,restriction,constrain,iteration,convergence) return(par.save.temp) } em_pcd_combination <- function(x, n, g, p, class, tol, restriction, constrain, iteration){ groupnum<-factorial(p+g-1)/(factorial(p)*factorial(g-1))-3 j<-1 pi_group<-1 rr<-rep(0,groupnum) id_group<-rep(0,(g^p)) id_group_diag<-rep(0,(g^p)) id<-rep(0,p) counter<-rep(0,groupnum*g) counter[1]<-2 counter[2]<-1 preresults<-.C("linking_matrix",as.integer(n),as.integer(g),as.integer(p), as.integer(j), as.integer(pi_group), as.integer(rr), as.integer(id), as.integer(counter), as.integer(id_group),as.integer(id_group_diag)) pi_group<-preresults[[5]] rr<-preresults[[6]] id_group<-preresults[[9]] pi<-double(g) alpha<-double(pi_group) mu <- double(p*g) sigma <- double(p*g) loglik <- double(1) convergence<-integer(1) ez<-double(g^p*n) j<-integer(1) id<-integer(p) pi <- rep(0.7/g, g) alpha <- rep(0.3/(g^p-g), pi_group) mu<-rep(c(-2,0,2), p) sigma<-rep(c(1.5,1,1.5),p) id<-rep(0,p) j<-1 results<-.C("em_pdim_pcd_new", as.integer(n), as.integer(g), as.integer(p),as.integer(j), as.integer(pi_group), as.integer(id), as.integer(id_group), as.integer(rr), as.double(ez), as.double(pi),as.double(alpha), as.double(x), as.double(mu), as.double(sigma),loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) pi<-results[[10]] alpha<-results[[11]] mu<-results[[13]] sigma<-results[[14]] loglik<-results[[15]] tol<-results[[16]] restriction<-results[[17]] constrain<-results[[18]] iteration<-results[[19]] convergence<-results[[20]] par.save.temp<-c(pi,alpha,mu,sigma,loglik,tol,restriction,constrain,iteration,convergence) return(par.save.temp) } ### Running samples ### #simulate z-scores set.seed(123) #true values N <- 1000 g <- 3 p <- 3 pi_true.vec<-c(0.100,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.200,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.100) mu1_true <- c(-2,0,2) mu2_true <- c(-2,0,2) mu3_true <- c(-2,0,2) sigmasq1_true <- c(1.5,1,1.5) sigmasq2_true <- c(1.5,1,1.5) sigmasq3_true <- c(1.5,1,1.5) x1 <- rep(0,N) x2 <- rep(0,N) x3 <- rep(0,N) for (k in 1:N){ u<-runif(1) if (u<=pi_true.vec[1]){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>pi_true.vec[1]&&u<=sum(pi_true.vec[1:2])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:2])&&u<=sum(pi_true.vec[1:3])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:3])&&u<=sum(pi_true.vec[1:4])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:4])&&u<=sum(pi_true.vec[1:5])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:5])&&u<=sum(pi_true.vec[1:6])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:6])&&u<=sum(pi_true.vec[1:7])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:7])&&u<=sum(pi_true.vec[1:8])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:8])&&u<=sum(pi_true.vec[1:9])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:9])&&u<=sum(pi_true.vec[1:10])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:10])&&u<=sum(pi_true.vec[1:11])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:11])&&u<=sum(pi_true.vec[1:12])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:12])&&u<=sum(pi_true.vec[1:13])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:13])&&u<=sum(pi_true.vec[1:14])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:14])&&u<=sum(pi_true.vec[1:15])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:15])&&u<=sum(pi_true.vec[1:16])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:16])&&u<=sum(pi_true.vec[1:17])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:17])&&u<=sum(pi_true.vec[1:18])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:18])&&u<=sum(pi_true.vec[1:19])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:19])&&u<=sum(pi_true.vec[1:20])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:20])&&u<=sum(pi_true.vec[1:21])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:21])&&u<=sum(pi_true.vec[1:22])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:22])&&u<=sum(pi_true.vec[1:23])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:23])&&u<=sum(pi_true.vec[1:24])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:24])&&u<=sum(pi_true.vec[1:25])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:25])&&u<=sum(pi_true.vec[1:26])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:26])&&u<=sum(pi_true.vec[1:27])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } } #some preliminaries for p=3 { x<-rep(0,3*N) for (l in 1:N){ x[3*l-2]<-x1[l] x[3*l-1]<-x2[l] x[3*l]<-x3[l] } class<-cbind(rep(1,N), rep(1,N), rep(1,N)) class[x1>0+1,1] <- 2 class[x1<0-1,1] <- 0 class[x2>0+1,2] <- 2 class[x2<0-1,2] <- 0 class[x3>0+1,3] <- 2 class[x3<0-1,3] <- 0 tol<-0.0001 constrain<-c(-1,0,1,-1,0,1,-1,0,1) iteration<-1 restriction<-1 g <- 3 p <- 3 } #need appropriate directory #load library #dyn.load("pdimpcdexchangeable.dll") dyn.load("pdimpcdcombination.dll") em_pcd_combination(x, N, g, p, class, tol, restriction, constrain, iteration) ### Running samples ### #simulate z-scores set.seed(321) #true values N <- 1000 g <- 3 p <- 3 pi_true.vec<-c(0.100,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.200,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.100) mu1_true <- c(-2,0,2) mu2_true <- c(-2,0,2) mu3_true <- c(-2,0,2) sigmasq1_true <- c(1.5,1,1.5) sigmasq2_true <- c(1.5,1,1.5) sigmasq3_true <- c(1.5,1,1.5) x1 <- rep(0,N) x2 <- rep(0,N) x3 <- rep(0,N) for (k in 1:N){ u<-runif(1) if (u<=pi_true.vec[1]){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>pi_true.vec[1]&&u<=sum(pi_true.vec[1:2])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:2])&&u<=sum(pi_true.vec[1:3])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:3])&&u<=sum(pi_true.vec[1:4])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:4])&&u<=sum(pi_true.vec[1:5])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:5])&&u<=sum(pi_true.vec[1:6])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:6])&&u<=sum(pi_true.vec[1:7])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:7])&&u<=sum(pi_true.vec[1:8])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:8])&&u<=sum(pi_true.vec[1:9])){ x1[k]<-rnorm(1,mean=mu1_true[1], sd=sqrt(sigmasq1_true[1])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:9])&&u<=sum(pi_true.vec[1:10])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:10])&&u<=sum(pi_true.vec[1:11])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:11])&&u<=sum(pi_true.vec[1:12])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:12])&&u<=sum(pi_true.vec[1:13])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:13])&&u<=sum(pi_true.vec[1:14])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:14])&&u<=sum(pi_true.vec[1:15])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:15])&&u<=sum(pi_true.vec[1:16])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:16])&&u<=sum(pi_true.vec[1:17])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:17])&&u<=sum(pi_true.vec[1:18])){ x1[k]<-rnorm(1,mean=mu1_true[2], sd=sqrt(sigmasq1_true[2])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:18])&&u<=sum(pi_true.vec[1:19])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:19])&&u<=sum(pi_true.vec[1:20])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:20])&&u<=sum(pi_true.vec[1:21])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[1], sd=sqrt(sigmasq2_true[1])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:21])&&u<=sum(pi_true.vec[1:22])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:22])&&u<=sum(pi_true.vec[1:23])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:23])&&u<=sum(pi_true.vec[1:24])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[2], sd=sqrt(sigmasq2_true[2])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } if (u>sum(pi_true.vec[1:24])&&u<=sum(pi_true.vec[1:25])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[1], sd=sqrt(sigmasq3_true[1])) } if (u>sum(pi_true.vec[1:25])&&u<=sum(pi_true.vec[1:26])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[2], sd=sqrt(sigmasq3_true[2])) } if (u>sum(pi_true.vec[1:26])&&u<=sum(pi_true.vec[1:27])){ x1[k]<-rnorm(1,mean=mu1_true[3], sd=sqrt(sigmasq1_true[3])) x2[k]<-rnorm(1,mean=mu2_true[3], sd=sqrt(sigmasq2_true[3])) x3[k]<-rnorm(1,mean=mu3_true[3], sd=sqrt(sigmasq3_true[3])) } } #some preliminaries for p=3 { x<-rep(0,3*N) for (l in 1:N){ x[3*l-2]<-x1[l] x[3*l-1]<-x2[l] x[3*l]<-x3[l] } class<-cbind(rep(1,N), rep(1,N), rep(1,N)) class[x1>0+1,1] <- 2 class[x1<0-1,1] <- 0 class[x2>0+1,2] <- 2 class[x2<0-1,2] <- 0 class[x3>0+1,3] <- 2 class[x3<0-1,3] <- 0 tol<-0.0001 constrain<-c(-1,0,1,-1,0,1,-1,0,1) iteration<-1 restriction<-1 g <- 3 p <- 3 } #need appropriate directory #load library dyn.load("pdimpcdexchangeable.dll") #dyn.load("pdimpcdcombination.dll") em_pcd_exchangeable(x, N, g, p, class, tol, restriction, constrain, iteration)