icp_ppm: Function that fits the multivariate independent product...

In ppmSuite: A Collection of Models that Employ Product Partition Distributions as a Prior on Partitions

Function that fits the multivariate independent product partition change point model

Description

icp_ppm is a function that fits a Bayesian product partition change point model. Each series is treated independently.

Usage

icp_ppm(ydata,
         a0, b0,
         mltypes,
         thetas,
         nburn, nskip, nsave,
         verbose = FALSE)

Arguments

ydata	An L \times n data matrix, where L is the number of time series and n, the number of time points.
a0	Vector of dimension L with shape 1 Beta parameters (see Details).
b0	Vector of dimension L with shape 2 Beta parameters (see Details).
mltypes	Type of marginal likelihood. Currently only available is: mltypes = 1. Observations within a block are conditionally independent Normal(\mu, \sigma^2) variates with mean \mu and variance \sigma^2. The desired marginal likelihood is obtained after integrating (\mu, \sigma^2) with respect to a Normal-Inverse-Gamma(\mu_0, \kappa_0, \alpha_0, \beta_0) prior.
thetas	An L \times q matrix containing hyperparameters associated with the marginal likelihood. The number of rows (L) corresponds to the number of series. The number of columns (q) depend on the marginal likelihood: If mltypes = 1, then q = 4 and thetas equals the hyperparameter (\mu_{0}, \kappa_{0}, \alpha_{0}, \beta_{0}) of the Normal-Inverse-Gamma prior.
nburn	The number of initial MCMC iterates to be discarded as burn-in.
nskip	The amount to thinning that should be applied to the MCMC chain.
nsave	Then number of MCMC iterates to be stored.
verbose	Logical indicating whether to print to screen the MCMC progression. The default value is verbose = FALSE.

Details

As described in Barry and Hartigan (1992) and Loschi and Cruz (2002), for each time series \boldsymbol{y}_{i} = (y_{i,1}, \ldots , y_{i,n})':

\boldsymbol{y}_{i} \mid \rho_{i} \sim \prod_{j = 1}^{b_{i}}\mathcal{F}(\boldsymbol{y}_{i,j} \mid \boldsymbol{\theta}_{i})

\rho_{i} \mid p_{i} \sim p_{i}^{b_{i} - 1}(1 - p_{i})^{n - b_{i}}

p_{i} \sim Beta(a_{i,0}, b_{i,0}).

Here, \rho_{i} = \{S_{i,1}, \ldots , S_{i,b_{i}}\} is a partition of the set \{1, \ldots , n\} into b_{i} contiguous blocks, and \boldsymbol{y}_{i,j} = (y_{i,t} : t \in S_{i,j})'. Also, \mathcal{F}( \cdot \mid \boldsymbol{\theta}_{i}) is a marginal likelihood function which depends on the nature of \boldsymbol{y}_{i}, indexed by a hyperparameter \boldsymbol{\theta}_{i}. Notice that p_{i} is the probability of observing a change point in series i, at each time t \in \{2, \ldots , n\}.

Value

The function returns a list containing arrays filled with MCMC iterates corresponding to model parameters. In order to provide more detail, in what follows let M be the number of MCMC iterates collected. The output list contains the following:

• C. An M \times \{L(n - 1)\} matrix containing MCMC iterates associated with each series indicators of a change point. The mth row in C is divided into L blocks; the first (n - 1) change point indicators for time series 1, the next (n - 1) change point indicators for time series 2, and so on.

• P. An M \times \{L(n - 1)\} matrix containing MCMC iterates associated with each series probability of a change point. The mth row in P is divided into L blocks; the first (n - 1) change point probabilities for time series 1, the next (n - 1) change point probabilities for time series 2, and so on.

Examples

# Generate data that has two series, each with 100 observations
y1 <- replicate(25, rnorm(4, c(-1, 0, 1, 2), c(0.1, 0.25, 0.5, 0.75)))
y2 <- replicate(25, rnorm(4, c(2, 1, 0, -2), c(0.1, 0.25, 0.5, 0.75)))
y <- rbind(c(t(y1)), c(t(y2)))
n <- ncol(y)
# Marginal likelihood parameters
thetas <- matrix(1, nrow = 2, ncol = 4)
thetas[1,] <- c(0, 1, 2, 1)
thetas[2,] <- c(0, 1, 2, 1)

# Fit the Bayesian ppm change point model
fit <- icp_ppm(ydata = y,
               a0 = c(1, 1),
               b0 = c(1, 1),
               mltypes = c(1, 1),
               thetas = thetas,
               nburn = 1000, nskip = 1, nsave = 1000)

cpprobsL <- matrix(apply(fit$C,2,mean), nrow=n-1, byrow=FALSE)

ppmSuite documentation built on July 26, 2023, 5:22 p.m.