####go to https://r-forge.r-project.org/R/?group_id=85

###download .zip file if you are using windows 
###or .tar.gz file if you are using unix-based systems
###then install from local file

###or if you are using unix-based systems directly run
### install.packages("RHmm", repos="http://R-Forge.R-project.org")


library("RHmm")

m1 <- c(1,1)
m2 <- c(-2, -2)
cov1 <- matrix(c(1, 1, 1, 4)/4, nrow=2)
cov2 <- matrix(c(1, -1, -1, 9)/4, nrow=2)
n_2d_2s <- distributionSet("NORMAL", mean=list(m1, m2),
                                        cov=list(cov1, cov2))

initProb2 <- rep(1,2)/2
transMat2 <- rbind(c(0.9, 0.1), c(0.2, 0.8))

hmm1 <- HMMSet(initProb2, transMat2, n_2d_2s)


x.sim = HMMSim(1000, hmm1)

plot(x.sim$obs[,1],x.sim$obs[,2],xlab="x",ylab="y")

bic = c()
for(n in 2:10){
	 tmp1 = HMMFit(x.sim$obs,nStates=n,control=list(nInit=10))
	bic = c(bic,tmp1$BIC)
	}

plot(bic)

hmmfit = HMMFit(x.sim$obs,nStates=2,control=list(nInit=10))
##look at the states
hmmfit$HMM

state.pred = viterbi(hmmfit, x.sim$obs)

plot(state.pred$states-x.sim$states,type="h")


#### real data

load("cnv.RData") ### you need to put the data under the working directory to load the data
plot(y)

hmmfit = HMMFit(y,nStates=4,control=list(nInit=10))
state.pred = viterbi(hmmfit, y)
plot(y)
lines(hmmfit$HMM$distribution$mean[state.pred$states],col="red")




hmmfit = HMMFit(y,nStates=5,control=list(nInit=10))
state.pred = viterbi(hmmfit, y)
plot(y)
lines(hmmfit$HMM$distribution$mean[state.pred$states],col="red")