NHMM (version 3.11)

NHMM_MVN: Bayesian Non-homogeneous Markov Model with Multivariate Normal emission distribution (NHMMMVN)

Description

NHMM_MVN calculates an NHMM for multiple sequences of data. The sequences can actually be short sets of equal length sequences (subseq). The traditional input variables (X) influence the non-homogenous transition probabiities of the model. An additional set of input variables (W) can be included to influence the mean of the emission distribution. All parameters are sampled via Gibbs steps (latent variables such that no tuning is needed.) The X variable coefficients are sampled through an unordered Multinomial logit model Polya-Gamma formulation (Polson, Scott, Windle). The hidden states are sampled through a blocked Gibbs sampler.

Usage

NHMM_MVN(
  y,
  subseq = NULL,
  X = NULL,
  betapriorm = NULL,
  betapriorp = NULL,
  K = 2,
  iters = 1000,
  burnin = 200,
  W = NULL,
  psipriorm = NULL,
  psipriorp = NULL,
  priors1 = NULL,
  priors2 = NULL,
  outdir = NULL,
  ymiss = FALSE,
  yrep = 0,
  ypred = 0,
  Xp = NULL,
  Wp = NULL,
  yhold = NULL
)

Arguments

y

T by J matrix of data (J=1 is sufficient) - missing data is denoted with NA

subseq

[optional] if y is actually a set of subsequences then give the length of those sequences (122 for JJAS) (365 is not it!). Default is subseq=T.

X

B by T matrix for the transition input data (B different inputs) Missing values are not allowed. If there are no Xs then use HMM function.

betapriorm

[optional] default=NULL which is reference prior. Or a [K+B by K matrix] for the mean of the Normal prior for the beta coefficients.

betapriorp

[optional] default=NULL which is reference prior. Or a [K+B by K+B by K array] for the precision prior(1/sig^2) of the Normal prior for the beta coefficients.

K

number of states (default=2)

iters

number of iterations to keep after burn in (default=1000)

burnin

the number of burn in (default=200)

W

[optional] is an A by T by J array of emission input data (A different inputs), missing values are not allowed, do not include an intercept term here The mean function depends on W.

psipriorm

[optional] default=NULL which is reference prior. Or a [K+A by J matrix] for the mean of the Normal prior for the beta coefficients.

psipriorp

[optional] default=NULL which is reference prior. Or a [K+A by J matrix] for the precision prior(1/sig^2) of the Normal prior for the beta coefficients.

priors1

[optional] scale parameter for each state (vector of length K) of the Wishart prior for Sigma priors1=NULL is set to 1

priors2

[optional] Covariance matrix [J by J by K] for the MVN for each state, second parameter of the Wishart for the Sigma matrix. priors2=NULL will be a diagonal matrix with diagonal set to 1

outdir

[optional] can output each set of parameters to output files in a directory use this with larger dimension data sets or with large number of iterations the output will be written line by line and not overburden the memory limit outdir needs to end with a slash or double slash depending on OS.

ymiss

[optional-TRUE/FALSE] if outdir is specified then draws for any missing data points will be saved to ymiss-J*.txt. There will be one data point for each sequence of length Tmiss (the amount missing from that sequence). Each row of the output file will be an iteration of the algorithm after burnin is removed.

yrep

[optional] number (ie. 100,200,or 500) of output replicate data sets to print to outdir. The replicates will be the same dimension as y and start after burn in period. default is zero. These replicate data sets are generated from the same input variables. - must be shorter than iters

ypred

number of predictive sets (ie. 0,100, 200, 500) - must be shorter than iters Predicted chains: will produce ypred set of predictives [pT by J] to print to outdir. yrep uses the same inputs to make replications but ypred uses new input values of X and W to make predictions over a different time span (pT) Also outputs a set of predictive z values or length (pT by ypred)

Xp

predictive set of Xs of length pT [B by pT] -missing values is not allowed -ensure that the Xp inputs are in the same order as X

Wp

predictive set of Ws of length pT [A by pT by J] -missing valus are not allowed -ensure that the Wp inputs are in the same order as W

yhold

[optioinal] a sequence of y observed values [pT by J], held out data that is of length ypred that is used to compute the predictive log score (PLS) which is a metric like BIC (ie. hold out last 10 missing data values are filled in with mean PLS value. PLS is the average PLS across sequences.

Value

my.nhmm object

Examples

# NOT RUN {
## Multivariate Normal data
# }
# NOT RUN {
data(NHMMdata)
attach(NHMMdata)

my.nhmm1=NHMM_MVN(y=ymvn,  X=tX, W=tW2, K=3, iters=50, burnin=10, 
                 priors1=rep(2,3))
OBIC(my.nhmm1)
zz=Oz(my.nhmm1)  #compare with the truth zgamma
qq=OQQ(my.nhmm1)
bb=OXcoef(my.nhmm1)
pp=OWcoef(my.nhmm1,FALSE)
tt=Oemparams(my.nhmm1,FALSE)  #just Sigma matrix for MVN, returns mean of Sigma
 
 
 #filelocation="C:\Users\iamrandom\Desktop\here\"
 #my.nhmm6=NHMM_MVN(y=ymvn[1:1800,],   X=matrix(tX[,1:1800],1,1800), 
 #          W=array(tW2[,1:1800,],dim=c(2,1800,15)), K=3, iters=50, 
 #            burnin=10,outdir=filelocation, ymiss=TRUE, yrep=10, 
 #            Xp=matrix(tX[,1801:2000],1,200), 
 #           Wp=array(tW2[,1801:2000,],dim=c(2,200,15)), ypred=10, 
 #           yhold=ymvn[1801:2000,])
 #OBIC(my.nhmm6)
 
 #Could try it with K=1, to compare K=1 to K=3
 
# }