install.packages("BCPmixture_1.0.tar.gz") library(BCPmixture) #simulate data mu0 <- 0.3 sigma0 <- 16 n <- 2000 set.seed(1) starts <- sample(1:n, 9) starts <- sort(starts) ends <- c(starts, 2000) starts <- c(0, starts) nums <- ends - starts while(any(nums < 25)){ starts <- sample(1:n, 9) starts <- sort(starts) ends <- c(starts, 2000) starts <- c(0, starts) nums <- ends - starts } all <- NULL for(i in 1:6){ onedim <- NULL mus <- runif(10, -2, 2) lamdas <- rep(0.8, 10) #round(rbeta(5, 1, 1), 1) zeros <- floor(nums*lamdas) for(j in 1:10){ x <- c(rnorm(nums[j] - zeros[j] - 1, mus[j], 0.2), rep(0, zeros[j])) onedim <- c(onedim, sample(x, nums[j]-1, replace = FALSE), rnorm(1, mus[j], 0.2)) } all <- rbind(all, onedim) } data <- t(all) change.points <- cumsum(nums) ###### run BCP am <- bm <- 1 w0 <- p0 <- 0.2 test <- BCP_mixture(data, time = 50, am = am, bm = bm, w0 = w0, p0 = p0) ###### takes about 50 mins estimated.change.points <- c(1:1998)[res[50, ] > 0.5] ###### plot result plot(data[, 1], col="red", cex=0.5, pch=20) points(data[, 2], col="blue", cex=0.5, pch=20) points(data[, 3], col="green", cex=0.5, pch=20) points(data[, 4], col="yellow", cex=0.5, pch=20) points(data[, 5], col="orange", cex=0.5, pch=20) points(data[, 6], col="purple", cex=0.5, pch=20) abline(v=estimated.change.points,lwd=2, col="black")