---
title: "Beware of autocorrelation"
author: ""
date: ""
output: html_document
---

```{r, setup, include=FALSE}
require(mosaic)   # Load additional packages here 

# Some customization.  You can alter or delete as desired (if you know what you are doing).
trellis.par.set(theme=theme.mosaic()) # change default color scheme for lattice
knitr::opts_chunk$set(
  tidy=FALSE,     # display code as typed
  size="small")   # slightly smaller font for code
```

- First we initialize some constants:
    + Number of samples(m) and samplesize(n)
    + Mean(mu) and standard deviation(si)
    + Correlation(ro)
```{r}
library(qcc)
m <- 50
n <- 8
mu <- 10
si <- 1
ro <- 0.9
```

- Our first observation(x) is drawn from a normal distribution
  with mean mu and standard deviation si:
```{r}
x <- rnorm(1, mu, si)
```

- We simulate a series of measurements `x` where `x[i+1] = mu * (1-ro) + ro*x[i] + e[i]`
    + The next value is a weighted average of the present value and the mean plus an error term `e[i]`.
    + All measurements are in-control, i.e. they have mean mu and standard deviation si.
    + But they are NOT independent unless ro=0. The correlation between successive measurements is ro.
    
```{r}
for(i in 1:(n*m-1)){
  x <- c(x, mu*(1-ro)+ro*x[i]+rnorm(1,0,si*sqrt(1-ro^2)))
}
```

- We format `data` as an m by n matrix, where each sample(row) has n consecutive values from `x`:
```{r}
data <- matrix(x, ncol = n, byrow = TRUE)
```

- We calculate and plot the standard `xbar` chart.
```{r}
qd <- qcc(data, type = "xbar", std.dev = "UWAVE-SD")
```

- What is your comment?

- Next, change the sign of the correlation: `ro <- -0.9`.

- Redo the simulation of `data` and qcc plot.

- What is your comment?

- Independence between observations in a sample is crucial.  
  We may check this assumption with the `acf` function.  
  We get a plot of the empirical correlation between `x[i]` and
  `x[i+1]`(Lag=1), the empirical correlation between
  `x[i]` and `x[i+2]`(Lag=2), etc.  
  These correlations should be below the dotted line in case of
  independence.
  The `pistonrings` data slightly violates the independence
  assumption.
```{r}
ax <- acf(x,lag=5)
data(pistonrings)
acf(pistonrings$diameter,lag=5)
```