simulateHMM <-
function(n,hmm){
  
  stateParameters = hmm@stateParameters
  if(hmm@internVariables$k!=hmm@internVariables$r){
    G = hmm@internVariables$G
    nuP = list(nu = hmm@transitionMatrix, pList = hmm@probList)  
    gamma = nuP2Gamma(nuPList=nuP,G=G)
    
  }
  else{
    gamma = hmm@transitionMatrix  
  }
  
  distribution = hmm@internVariables$distribution
  
  
  #no of states
  m=dim(gamma)[1]
  
  
  randomSample = function(n,parameter,distribution){
    if(distribution=="normal"){
      return(rnorm(n=n,mean=parameter[1],sd=parameter[2]))
    }
    
    if(distribution=="mv_normal"){
      return(rmnorm(n=n,mean=parameter$mu,varcov=parameter$sigma))
    }  
  }
  
  #reparametrisation
  stateParameters = extern2intern(stateParameters,distribution,m)
  
  #stationary distr.
  delta = stat_distr(gamma)
  
  #sample first state
  s = sample(x=1:m,size=1,replace=T,prob=delta)
  
  if(distribution == "normal")
    D=1
  if(distribution == "mv_normal")
    D = length(stateParameters[[1]]$mu)  
  
  
  #initialise data vector
  data=matrix(nrow=n,ncol=D)
  
  #sample first obs.
  data[1,] = randomSample(n=1,parameter=stateParameters[[s[1]]],distribution)
  
  #sample path
  for(i in 2:n){
    s[i] = sample(x=1:m,size=1,replace=T,prob=gamma[s[i-1],])
    data[i,] = randomSample(n=1,parameter=stateParameters[[s[i]]],distribution)
  }
  
  if(distribution == "normal")
    data = as.vector(data)
  
  return(list(data = data,s=s))
  
}
