# R/equil.R In csgillespie/HMMs: Hidden Markov models to analyse protein sequences

#### Documented in equilhmm.simrdiric

```####################################################
#Private functions
####################################################
equil <-
function(P) {
ev = eigen(t(P))\$vectors[,1]
return(Re(ev/sum(ev)))
}

rdiric <-
function(n, a) {
p <- length(a)
m <- matrix(nrow=n, ncol=p)
for (i in 1:p) {
m[,i] <- rgamma(n,a[i])
}
sumvec <- m %*% rep(1,p)
m / as.vector(sumvec)
}

####################################################
#Public functions
####################################################
#' Simulates a hidden Markov model
#'
#' @aliases equil rdiric
#' @param n Length of the observed sequence
#' @param lambda The hidden sequence transition matrix
#' @param P An array of transition matrices for observed sequence
#' @param f The number of states of observed sequence
#' @author Nina Wilkinson
#' @return A list
#' @keywords character
#' @export

hmm.sim <-
function(n,lambda,P,f)
{

## first simulate the hidden states (segmentation)
r = dim(lambda)[1]
pi.lam = equil(lambda)
s = numeric(n)
s[1] = sample(1:r,1,replace=TRUE,prob=pi.lam)
for(t in 2:n){
s[t] = sample(1:r,1,replace=TRUE,prob=lambda[match(s[t-1],1:r),])
}
## then simulate the observed states (conditional on s)
pi.P = array(0,c(f,r))
for(i in 1:r){
pi.P[,i] = equil(P[,,i])
}
y = numeric(n)
y[1] <- sample(1:f,1,replace=TRUE,prob=pi.P[,s[1]])
for(t in 2:n){
y[t] <- sample(1:f,1,replace=TRUE,prob=P[y[t-1],,s[t]])
}
list(s=s,y=y)
}
```
csgillespie/HMMs documentation built on May 14, 2019, 12:11 p.m.