gn2<-function(theta, merge2) { pix1 <- theta[1:3] piy1 <- theta[4:6] f00 <- theta[7:9] ss <- theta[10:15] rho1 <- matrix(theta[16:(15 + length(merge2))], , 2) pix2 <- c(pix1, 1 - sum(pix1)) Pix2 <- diag(pix2) piy2 <- c(piy1, 1 - sum(piy1)) Piy2 <- diag(piy2) f0 <- c(f00, 1 - sum(f00)) if(sum(abs(pix2)) > 1) { return(0) } if(sum(abs(piy2)) > 1) { return(0) } if(sum(abs(f0)) > 1) { return(0) } Sx1 <- Smatrix(ss, pix2) Sy1 <- Smatrix(ss, piy2) k <- dim(merge2)[1] if(k == 1) { F1 <- Fmatrix(theta[16], theta[17], f0, (rho1[1, 1] * Sx1), (rho1[1, 2] * Sy1), Pix2, Piy2) return(F1) } else { merge <- matrix(0, k, 2) for(i in 1:k) { merge[i, ] <- rev(sort(merge2[i, ])) } me <- merge[k, ] f1 <- f0 m <- me cc <- 0 if(merge2[, 1][k] > 0) { cc <- merge[, 1][k] v <- matrix(0, k, 1) v[1] <- merge2[k, 1] cur <- 0 max <- 1 while(cur < max) { cur <- cur + 1 for(i in c(1, 2)) { if(merge2[(v[cur]), i] > 0) { max <- max + 1 v[max] <- merge2[ (v[cur]), i] } } cc <- matrix(0, max, 1) for(i in 1:max) { cc[i] <- v[i] } } } for(i in 1:k) { iii <- k - i + 1 if(iii == k) { if(any(merge2[k, 1] == cc)) { pix <- pix2 Pix <- Pix2 piy <- piy2 Piy <- Piy2 Sx2 <- Sx1 Sy2 <- Sy1 } else { piy <- pix2 Piy <- Pix2 pix <- piy2 Pix <- Piy2 Sy2 <- Sx1 Sx2 <- Sy1 } } else { if(any(iii == cc)) { pix <- pix2 Pix <- Pix2 piy <- pix2 Piy <- Pix2 Sx2 <- Sx1 Sy2 <- Sx1 } else { piy <- piy2 Piy <- Piy2 pix <- piy2 Pix <- Piy2 Sx2 <- Sy1 Sy2 <- Sy1 } } m <- rev(sort(me)) mm <- merge[m[1], ] t1 <- rho1[iii, ][1] t2 <- rho1[iii, ][2] hp <- rev(order(me)) for(j in 1:(4^(i - 1))) { if(j == 1) { F1 <- Fmatrix(t1, t2, f1[ 1:4], Sx2, Sy2, Pix, Piy) } else { F1 <- Fmatrix(t1, t2, f1[ (4 * (j - 1) + 1): (4 * (j - 1) + 4)], Sx2, Sy2, Pix, Piy) } if(j == 1) { f <- c(as.vector(F1), f1[ -1:-4]) } else { f <- c(f[1:(4^2 * (j - 1))], as.vector(F1), f1[ -1: - (4 * j)]) } } f12 <- array(f, c(rep(4, i + 1))) f112 <- aperm(f12, c(hp)) f1 <- as.vector(f112) if(m[1] > 0) t0 <- 9999 else { return(f112) } me <- c(mm, m[-1]) } } }