File: higher_order_markov_chains.Rmd

package info (click to toggle)
r-cran-markovchain 0.8.5-4-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 2,060 kB
  • sloc: cpp: 2,854; sh: 13; makefile: 2
file content (295 lines) | stat: -rw-r--r-- 12,459 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
---
title:
  plain:     "Higher, possibly multivariate, Order Markov Chains in markovchain package"
  formatted: "Higher, possibly multivariate, Order Markov Chains in \\pkg{markovchain} package"
  short:     "Higher order (multivariate) Markov chains"
pagetitle: "Higher order (multivariate) Markov chains"
author:
  - name: "Deepak Yadav"
    affiliation: B-Tech student, Computer Science and Engineering
    address: >
      Indian Institute of Technology, Varanasi 
      Uttar Pradesh - 221 005, India
    email: \email{deepakyadav.iitbhu@gmail.com}
  - name: "Tae Seung Kang"
    affiliation: Ph.D student, Computer \& Information Science \& Engineering
    address: >
      University of Florida
      Gainesville, FL, USA
    email: \email{tskang3@gmail.com}
  - name: "Giorgio Alfredo Spedicato"
    affiliation: Ph.D C.Stat ACAS, UnipolSai R\&D
    address: >
      Via Firenze 11
      Paderno Dugnano 20037 Italy
    email: \email{spedicato\_giorgio@yahoo.it}
preamble: >
  \author{Deepak Yadav, Tae Seung Kang, Giorgio Alfredo Spedicato}
  \usepackage{graphicx}
  \usepackage{amsmath}
  \usepackage{tabularx}
  \usepackage{longtable}
  \usepackage{booktabs}
  \setkeys{Gin}{width=0.8\textwidth}
abstract: |
  The \pkg{markovchain} package contains functions to fit higher (possibly) multivariate order Markov chains. The functions are shown as well as simple exmaples
output: if (rmarkdown::pandoc_version() < 2.2) function(...) { rmarkdown::pdf_document(template = "./template.tex", ...) } else function(...) { bookdown::pdf_book(base_format = rticles::jss_article, ...) }
vignette: >
  %\VignetteIndexEntry{Higher order markov chains} 
  %\VignetteEngine{knitr::rmarkdown}  
  %VignetteEncoding{UTF-8}
keywords:
  plain: [Higher order Markov chains]
  formatted: [Higher order Markov chains]
documentclass: jss
classoption: nojss
bibliography: markovchainBiblio.bib
pkgdown:
  as_is: true
  extension: pdf
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width=8.5, fig.height=6, out.width = "70%")
set.seed(123)
```

# Higher Order Markov Chains

Continuous time Markov chains are discussed in the CTMC vignette which is a part of the package.

An experimental `fitHigherOrder` function has been written in order to fit a higher order Markov chain ([@ching2013higher, @ching2008higher]). 
`fitHigherOrder` takes two inputs

  1. sequence: a categorical data sequence.
  2. order: order of Markov chain to fit with default value 2.

The output will be a `list` which consists of

  1. lambda: model parameter(s).
  2. Q: a list of transition matrices. $Q_i$ is the $ith$ step transition matrix stored column-wise.
  3. X: frequency probability vector of the given sequence.

Its quadratic programming problem is solved using `solnp` function of the Rsolnp package [@pkg:Rsolnp]. 

```{r load, results='hide', warning=FALSE, message=FALSE}
require(markovchain)
library(Rsolnp)
```

```{r higherOrder}
data(rain)
fitHigherOrder(rain$rain, 2)
fitHigherOrder(rain$rain, 3)
```

# Higher Order Multivariate Markov Chains

## Introduction

HOMMC model is used for modeling behaviour of multiple categorical sequences generated by similar sources. The main reference is [@ching2008higher]. Assume that there are s categorical sequences and each has possible states in M. In nth order MMC the state probability distribution of the jth sequence at time $t = r + 1$ depend on the state probability distribution of all the sequences (including itself) at times $t = r, r - 1, ..., r - n + 1$.

\[
x_{r+1}^{(j)} = \sum_{k=1}^{s}\sum_{h=1}^{n}\lambda_{jk}^{(h)}P_{h}^{(jk)}x_{r-h+1}^{(k)}, j = 1, 2, ..., s, r = n-1, n, ...
\]

with initial distribution $x_{0}^{(k)}, x_{1}^{(k)}, ... , x_{n-1}^{(k)} (k = 1, 2, ... , s)$. Here

\[
\lambda _{jk}^{(h)} \geq 0, 1\leq j, k\leq s, 1\leq h\leq  n \enspace and \enspace \sum_{k=1}^{s}\sum_{h=1}^{n} \lambda_{jk}^{(h)} = 1, j = 1, 2, 3, ... , s.
\]

Now we will see the simpler representation of the model which will help us understand the result of `fitHighOrderMultivarMC` method.

\vspace{5mm}

Let $X_{r}^{(j)} = ((x_{r}^{(j)})^{T}, (x_{r-1}^{(j)})^{T}, ..., (x_{r-n+1}^{(j)})^{T})^{T} for \enspace j = 1, 2, 3, ... , s.$ Then

\vspace{5mm}

\[
\begin{pmatrix}
X_{r+1}^{(1)}\\ 
X_{r+1}^{(2)}\\ 
.\\ 
.\\ 
.\\ 
X_{r+1}^{(s)}
\end{pmatrix} = \begin{pmatrix}
 B^{11}&  B^{12}&  .&  .&  B^{1s}& \\ 
 B^{21}&  B^{22}&  .&  .&  B^{2s}& \\ 
 .&  .&  .&  .&  .& \\ 
 .&  .&  .&  .&  .& \\ 
 .&  .&  .&  .&  .& \\ 
 B^{s1}&  B^{s2}&  .&  .&  B^{ss}& \\ 
\end{pmatrix} \begin{pmatrix}
X_{r}^{(1)}\\ 
X_{r}^{(2)}\\ 
.\\ 
.\\ 
.\\ 
X_{r}^{(s)}
\end{pmatrix} \textrm{where}
\]

\[B^{ii} = \begin{pmatrix}
 \lambda _{ii}^{(1)}P_{1}^{(ii)}&  \lambda _{ii}^{(2)}P_{2}^{(ii)}&  .&  .&  \lambda _{ii}^{(n)}P_{n}^{(ii)}& \\ 
 I&  0&  .&  .&  0& \\ 
 0&  I&  .&  .&  0& \\ 
 .&  .&  .&  .&  .& \\ 
 .&  .&  .&  .&  .& \\ 
 0&  .&  .&  I&  0& 
\end{pmatrix}_{mn*mn} \textrm{and}
\]

\vspace{5mm}

\[
B^{ij} = \begin{pmatrix}
 \lambda _{ij}^{(1)}P_{1}^{(ij)}&  \lambda _{ij}^{(2)}P_{2}^{(ij)}&  .&  .&  \lambda _{ij}^{(n)}P_{n}^{(ij)}& \\ 
 0&  0&  .&  .&  0& \\ 
 0&  0&  .&  .&  0& \\ 
 .&  .&  .&  .&  .& \\ 
 .&  .&  .&  .&  .& \\ 
 0&  .&  .&  0&  0& 
\end{pmatrix}_{mn*mn} \textrm{when } i\neq j.
\]

\vspace{5mm}

## Representation of parameters in the code

$P_{h}^{(ij)}$ is represented as $Ph(i,j)$ and $\lambda _{ij}^{(h)}$ as Lambdah(i,j). 
For example: $P_{2}^{(13)}$ as $P2(1,3)$ and $\lambda _{45}^{(3)}$ as Lambda3(4,5).

## Definition of HOMMC class

```{r hommcObject}
showClass("hommc")
```

Any element of `hommc` class is comprised by following slots:

  1. states: a character vector, listing the states for which transition probabilities are defined.
  2. byrow: a logical element, indicating whether transition probabilities are shown by row or by column.
  3. order: order of Multivariate Markov chain.
  4. P: an array of all transition matrices.
  5. Lambda: a vector to store the weightage of each transition matrix.
  6. name: optional character element to name the HOMMC

## How to create an object of class HOMMC

```{r hommcCreate}
states <- c('a', 'b')
P <- array(dim = c(2, 2, 4), dimnames = list(states, states))
P[ , , 1] <- matrix(c(1/3, 2/3, 1, 0), byrow = FALSE, nrow = 2, ncol = 2)

P[ , , 2] <- matrix(c(0, 1, 1, 0), byrow = FALSE, nrow = 2, ncol = 2)

P[ , , 3] <- matrix(c(2/3, 1/3, 0, 1), byrow = FALSE, nrow = 2, ncol = 2)

P[ , , 4] <- matrix(c(1/2, 1/2, 1/2, 1/2), byrow = FALSE, nrow = 2, ncol = 2)

Lambda <- c(.8, .2, .3, .7)

hob <- new("hommc", order = 1, Lambda = Lambda, P = P, states = states, 
           byrow = FALSE, name = "FOMMC")
hob
```

## Fit HOMMC

`fitHighOrderMultivarMC` method is available to fit HOMMC. Below are the 3 parameters of this method.

  1. seqMat: a character matrix or a data frame, each column represents a categorical sequence.
  2. order: order of Multivariate Markov chain. Default is 2.
  3. Norm: Norm to be used. Default is 2.

# A Marketing Example

We tried to replicate the example found in [@ching2008higher] for an application of HOMMC. A soft-drink company in Hong Kong is facing an in-house problem of production planning and inventory control. A pressing issue is the storage space of its central warehouse, which often finds itself in the state of overflow or near capacity. The company is thus in urgent needs to study the interplay between the storage space requirement and the overall growing sales demand. The product can be classified into six possible states (1, 2, 3, 4, 5, 6) according to their sales volumes. All products are labeled as 1 = no sales volume, 2 = very slow-moving (very low sales volume), 3 = slow-moving, 4 = standard, 5 = fast-moving or 6 = very fast-moving (very high sales volume). Such labels are useful from both marketing and production planning points of view. The data is cointaind in `sales` object.

```{r hommsales}
data(sales)
head(sales)
```

The company would also like to predict sales demand for an important customer in order to minimize its inventory build-up. More importantly, the company can understand the sales pattern of this customer and then develop a marketing strategy to deal with this customer. Customer's sales demand sequences of five important products of the company for a year. We expect sales demand sequences generated by the same customer to be correlated to each other. Therefore by exploring these relationships, one can obtain a better higher-order multivariate Markov model for such demand sequences, hence obtain better prediction rules.

In [@ching2008higher] application, they choose the order arbitrarily to be eight, i.e., n = 8. We first estimate all the transition probability matrices $P_{h}^{ij}$ and we also have the estimates of the stationary probability distributions of the five products:.

$\widehat{\boldsymbol{x}}^{(1)} = \begin{pmatrix}
0.0818&  0.4052&  0.0483&  0.0335& 0.0037& 0.4275 
\end{pmatrix}^{\boldsymbol{T}}$

$\widehat{\boldsymbol{x}}^{(2)} = \begin{pmatrix}
0.3680&  0.1970&  0.0335&  0.0000& 0.0037& 0.3978 
\end{pmatrix}^{\boldsymbol{T}}$

$\widehat{\boldsymbol{x}}^{(3)} = \begin{pmatrix}
0.1450&  0.2045&  0.0186&  0.0000& 0.0037& 0.6283 
\end{pmatrix}^{\boldsymbol{T}}$

$\widehat{\boldsymbol{x}}^{(4)} = \begin{pmatrix}
0.0000&  0.3569&  0.1338&  0.1896& 0.0632& 0.2565
\end{pmatrix}^{\boldsymbol{T}}$

$\widehat{\boldsymbol{x}}^{(5)} = \begin{pmatrix}
0.0000&  0.3569&  0.1227&  0.2268& 0.0520& 0.2416
\end{pmatrix}^{\boldsymbol{T}}$

By solving the corresponding linear programming problems, we obtain the following higher-order multivariate Markov chain model:

\vspace{3mm}

$\boldsymbol{x}_{r+1}^{(1)} = \boldsymbol{P}_{1}^{(12)}\boldsymbol{x}_{r}^{(2)}$

$\boldsymbol{x}_{r+1}^{(2)} = 0.6364\boldsymbol{P}_{1}^{(22)}\boldsymbol{x}_{r}^{(2)} + 
0.3636\boldsymbol{P}_{3}^{(22)}\boldsymbol{x}_{r}^{(2)}$

$\boldsymbol{x}_{r+1}^{(3)} = \boldsymbol{P}_{1}^{(35)}\boldsymbol{x}_{r}^{(5)}$

$\boldsymbol{x}_{r+1}^{(4)} = 0.2994\boldsymbol{P}_{8}^{(42)}\boldsymbol{x}_{r}^{(2)} + 
0.4324\boldsymbol{P}_{1}^{(45)}\boldsymbol{x}_{r}^{(5)} + 0.2681\boldsymbol{P}_{2}^{(45)}\boldsymbol{x}_{r}^{(5)}$

$\boldsymbol{x}_{r+1}^{(5)} = 0.2718\boldsymbol{P}_{8}^{(52)}\boldsymbol{x}_{r}^{(2)} + 
0.6738\boldsymbol{P}_{1}^{(54)}\boldsymbol{x}_{r}^{(4)} + 0.0544\boldsymbol{P}_{2}^{(55)}\boldsymbol{x}_{r}^{(5)}$

\vspace{3mm}

According to the constructed 8th order multivariate Markov model, Products A and B are closely related. In particular, the sales demand of Product A depends strongly on Product B. The main reason is that the chemical nature of Products A and B is the same, but they have different packaging for marketing purposes. Moreover, Products B, C, D and E are closely related. Similarly, products C and E have the same product flavor, but different packaging. In this model, it is interesting to note that both Product D and E quite depend on Product B at order of 8, this relationship is hardly to be obtained in conventional Markov model owing to huge amount of parameters. The results show that higher-order multivariate Markov model is quite significant to analyze the relationship of sales demand.

```{r hommcFit, warning = FALSE, message = FALSE}
# fit 8th order multivariate markov chain
object <- fitHighOrderMultivarMC(sales, order = 8, Norm = 2)
```

We choose to show only results shown in the paper. We see that $\lambda$ values are quite close, but not equal, to those shown in the original paper.

```{r result, echo = FALSE}
i <- c(1, 2, 2, 3, 4, 4, 4, 5, 5, 5)
j <- c(2, 2, 2, 5, 2, 5, 5, 2, 4, 5)
k <- c(1, 1, 3, 1, 8, 1, 2, 8, 1, 2)

if(object@byrow == TRUE) {
    direction <- "(by rows)" 
} else {
    direction <- "(by cols)" 
}

cat("Order of multivariate markov chain =", object@order, "\n")
cat("states =", object@states, "\n")

cat("\n")
cat("List of Lambda's and the corresponding transition matrix", direction,":\n")

for(p in 1:10) {
  t <- 8*5*(i[p]-1) + (j[p]-1)*8
  cat("Lambda", k[p], "(", i[p], ",", j[p], ") : ", object@Lambda[t+k[p]],"\n", sep = "")
  cat("P", k[p], "(", i[p], ",", j[p], ") : \n", sep = "")
  print(object@P[, , t+k[p]])
  cat("\n")
}
```

# References