simapp<-function(theta, seqLength, merge1)
{
    pix1 <- theta[1:3]
    piy1 <- theta[4:6]
    f00 <- theta[7:9]
    ss <- theta[10:15]
    rho <- theta[16]
    pix2 <- c(pix1, 1 - sum(pix1))
    Pix2 <- diag(pix2)
    piy2 <- c(piy1, 1 - sum(piy1))
    Piy2 <- diag(piy2)
    f0 <- c(f00, 1 - sum(f00))
    Sx1 <- Smatrix(ss, pix2)
    Sy1 <- rho * Smatrix(ss, piy2)
    merge2 <- merge1
    k <- dim(merge2)[1]
    height2 <- c(theta[17:length(theta)], 1)
    height <- height2
    h <- c(0, height)
    me <- merge2[k,  ]
    m <- me
    t0 <- h[k + 1]
    cc <- 0
    if(merge2[, 1][k] > 0) {
        cc <- merge2[, 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]
            }
        }
    }
    seqr <- matrix(0, seqLength, k - 1)
    fs <- matrix(0, seqLength, k + 1)
    seq1 <- cbind(sample(1:4, seqLength, prob = f0, replace = T)
        )
    for(i in 1:k) {
        iii <- k - i + 1
        if(iii == k) {
            xx <- 0
            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 <- merge2[m[1],  ]
        if(me[1] > 0) {
            t1 <- t0 - h[me[1] + 1]
        }
        else {
            t1 <- t0
        }
        if(me[2] > 0) {
            t2 <- t0 - h[me[2] + 1]
        }
        else {
            t2 <- t0
        }
        if(t1 < 0 || t1 > 1) {
            return(0)
        }
        if(t2 < 0 || t2 > 1) {
            return(0)
        }
        hp <- rev(order(me))
        if(i != 1) {
            seq1 <- seqr[, (k - i + 1)]
        }
        seq2 <- genseq3(t1, t2, seq1, Sx2, Sy2, Pix, Piy)
        if(me[1] < 0) {
            fs[, (abs(me[1]))] <- seq2[, 1]
        }
        else {
            seqr[, (abs(me[1]))] <- seq2[, 1]
        }
        if(me[2] < 0) {
            fs[, (abs(me[2]))] <- seq2[, 2]
        }
        else {
            seqr[, (abs(me[2]))] <- seq2[, 2]
        }
        if(m[1] > 0)
            t0 <- h[m[1] + 1]
        else {
            return(fs)
        }
        me <- c(mm, m[-1])
    }
}

genseq3 <- function(t1, t2, seq, Sx2, Sy2, Pix, Piy) {
    sequ <- matrix(0, nc = 2, nr = length(seq))
    seqx<-seq
    seqy<-seq
    SS <- diag(c(1, 1, 1, 1))
    p1 <- SS+Sx2%*%Pix
    p2 <- SS+Sy2%*%Piy
    print(p1)
     for(i in 1:(length(seq) * t1)) {
        km <- sample(1:length(seq), 1)
        seqx <- replace(seqx, km, (sample(1:4, 1, prob = p1[seqx[km],  ],replace=T)))
    }
    for(i in 1:(length(seq) * t2)) {
        km <- sample(1:length(seq), 1)
        seqy <- replace(seqy, km, (sample(1:4, 1, prob = p2[seqy[km],  ],replace=T)))
    }
    sequ<-cbind(seqx,seqy)
    return(sequ)
}