---
title: "Statistics and electronics - lecture 1"
author: "The ASTA team"
output:
  slidy_presentation:
    fig_caption: no
    highlight: tango
    theme: cerulean
  pdf_document:
    fig_caption: no
    keep_tex: yes
    highlight: tango
    number_section: yes
    toc: yes

---

```{r, include = FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
#options(digits = 2)
## Remember to add all packages used in the code below!
missing_pkgs <- setdiff(c("mosaic","VennDiagram"), rownames(installed.packages()))
if(length(missing_pkgs)>0) install.packages(missing_pkgs)
library(VennDiagram)
library(mosaic)

```
----

## Sources of variation

We shall study 2 types of variation


- measurement variation due to random errors on a measuring device
- component variation due to random errors in the production proces

```{r, fig.width=10, echo=FALSE, fig.align='center', fig.caption="Capacitors with nominal value 47/150 nF and tolerance 1%"}
url <- "https://asta.math.aau.dk/static-files/asta/img/kondensatorer-forsog-nr2.jpg"
z <- tempfile()
download.file(url, z, mode = "wb")
grid::grid.raster(jpeg::readJPEG(z))
invisible(file.remove(z))
```

## Data from Peter Koch

Peter has done 100 independent measurements of the capacity of 4 of the displayed capacitors and one additional. 
Nominal values are 47, 47, 100, 150, 150. All with stated tolerance of 1%.

```{r }
load(url("https://asta.math.aau.dk/datasets?file=cap_1pct.RData"))
head(capDat, 4)
```

Here we see the first 4 capacity measurements of the first capacitor with nominal value 47.

- Remark: The measured values are consistently below the nominal value minus the 1% tolerance: $47 - 0.47 = 46.53$. 

```{r}
table(capDat$sample)
```


## Transformation

Linearisation:

\begin{align}
  f(x) &\approx f(x_0) + f'(x_0)(x - x_0) \\
  \\
  x_0 &= 1 \\
  f(x) &= \log x \\
  f'(x) &= 1/x \\
  \\
  x    &= m/n \\
  \log \left ( \frac{m}{n} \right ) &\approx \log 1 + \frac{1}{1} \left ( \frac{m}{n} - 1 \right ) \\
  &= \frac{m - n}{n}
\end{align}

\vspace{1em}

$$
\log \left ( \frac{m}{n} \right ) \approx \frac{m - n}{n}
$$

```{r}
n <- 47
m <- seq(47-5*0.01*47, 47+5*0.01*47, length.out = 100)
plot(m, log(m/n), col = "red", type = "l")
lines(m, (m - n)/n, col = "blue", type = "l")
legend("topleft", legend = c("log(m/n)", "(m-n)/n"), lty = 1, col = c("red", "blue"))
```



## Transformation

$$
\log \left ( \frac{m}{n} \right ) \approx \frac{m - n}{n}
$$


Instead of the raw measurement we will consider:

$\mbox{lnError = ln(measuredValue/nominalValue)}$

Remark that by linear approximation:

$\mbox{lnError}\approx\mbox{measuredValue/nominalValue - 1 = }$ $\mbox{(measuredValue-nominalValue)/nominalValue}$

which is the error relative to the nominal value.

**I.e.: `lnError` can be interpreted as relative error.**

## Transformed data

```{r fig.dim=c(4,3)}
capDat = within(capDat, lnError <- log(capacity/nomval))
head(capDat, 2)
tail(capDat, 2)
```

- The resolution on Peters capacitance meter is with 2/1 decimal(s) in the 47/150 nF range, which means that only a limited number of different values(3-8) are observed. Meaning that box- or histogram-plots are noninformative.

## Model considerations

- The measurements are more than 2.7% below the nominal value.
This must be due to a systematic error on the meter.

In this case we have as earlier mentioned two further sources of error:

- $\mbox{ln(measuredValue / nominalValue) = systematicError + productionError + measurementError}$

## Statistical model

$$\mbox{ln(measuredValue / nominalValue) = systematicError + productionError + measurementError}$$

We formulate the model:

- $Y_{ij}=\mu+A_i+\varepsilon_{ij}$

where

- $Y_{ij}$ is the log error measurement
- $\mu$ is the systematic error on the meter
- $A_i$ is the random production error
- $\varepsilon_{ij}$ is the random measurement error
- $i=1,2,3,4,k=5$ is the number of the 5 samples
- $j=1,\ldots,n=100$ is the number of the observation in each sample

## Assumptions

This is the model treated in WMM chapter 13.11, where it is assumed that

- $A_i$ is normally distributed with mean 0 and variance $\sigma_\alpha^2$, which represents the production error
- $\varepsilon_{ij}$ is normally distributed with mean 0 and variance $\sigma^2$, which represents the measurement error

## Estimation of systematic error

The systematic error is simply estimated by the mean

- $\hat\mu=\bar y_{..}$

```{r fig.dim=c(4,3)}
muhat <- mean(capDat$lnError)
muhat
```

The meter systematically reports a value, which is estimated to be `r round(-100*muhat, 2)`% too low.

## Estimation of random error

Notation from WMM chapter 13.3:

- $SSA=n\sum_i (\bar y_{i.}-\bar y_{..})^2$ (related to production error)
- $SSE=\sum_{ij} ( y_{ij}-\bar y_{i.})^2$ (related to measurement error)

Theorem 13.4 states:

- $E(SSA)=(k-1)\sigma^2+n(k-1)\sigma_\alpha^2$
- $E(SSE)=k(n-1)\sigma^2$

## Fit

```{r}
fit <- lm(lnError ~ sample, data = capDat)
anova(fit)
```

where we read

```{r}
SS <- anova(fit)$`Sum Sq`
SSA <- SS[1]
SSE <- SS[2]
```

- SSA = `r format(SSA, digits = 3)` and SSE = `r format(SSE, digits = 3)`

## Solution

Solving the equations 

- `SSA = E(SSA)` and `SSE = E(SSE)`

yields

```{r, include=FALSE}
k <- length(unique(capDat$sample))
n <- nrow(capDat)/k
sigma2 <- SSE/(k*(n-1))
sigma2a <- (1/n)*(SSA/(k-1) - sigma2)

mult <- ceiling(min(log10(sigma2a), log10(sigma2)))

fmt <- function(x) {
  paste0(round(x*10^(-mult), 2), 
         "e", mult, "")
}
fmt_tex <- function(x) {
  paste0(round(x*10^(-mult), 2), 
         " \\times 10^{", mult, "}")
}
fmt(sigma2)
fmt_tex(sigma2)

fmt(sigma2a)
fmt_tex(sigma2a)
```

- $\hat\sigma_\alpha^2=\frac{1}{n}(\frac{SSA}{k-1}-\hat\sigma^2)$ = $`r fmt_tex(sigma2a)`$
- $\hat\sigma^2=\frac{SSE}{k(n-1)}$ = $`r fmt_tex(sigma2)`$


## Summing up

- the meter has an estimated systematic error of $`r round(100*muhat, 2)`$%
- the estimated standard error of the meter is $\sqrt{`r fmt_tex(sigma2)`}$ = `r round(100*sqrt(sigma2), 3)`%
- the estimated standard error of the production is $\sqrt{`r fmt_tex(sigma2a)`}$ = `r round(100*sqrt(sigma2a), 2)`%. 
So the 3-sigma limit is `r round(3*100*sqrt(sigma2a), 2)`%, which is in accordance with the tolerance of 1%. It should be noted that the estimate is insecure, as it is based on 4 degrees of freedom only.

The estimated variance on log error

- $`r fmt_tex(sigma2)` + `r fmt_tex(sigma2a)` = `r fmt_tex(sigma2 + sigma2a)`$

is clearly dominated by the production error.

## Test of no random effect

We have the possibility of testing the hypothesis

- $H_0$: $\sigma_\alpha=0$

This is equivalent to 

- $E(SSA/(k-1))=E(SSE/k/(n-1))=\sigma^2$

Under $H_0$ the statistic 

- $F=\frac{\frac{SSA}{k-1}}{\frac{SSE}{k(n-1)}}$

has an F-distribution with degrees of freedom $(k-1,k(n-1))$

In the actual case $f_{obs}=4067.4$, which is highly significant (p-value=0).

## Lognormal variation

In the preceeding we assumed normal errors after a log  transformation.

Let $X$ be a random variable and $Y=ln(X)$.


We say that $X$ has a lognormal distribution if $Y$ has a normal distribution with - say - mean $\mu$ and standard deviation $\sigma$.

Density plots:

```{r echo=FALSE, fig.dim = c(8, 6)}
curve(dlnorm(x,0,1),0,5,ylab="density",lwd=2,col=2)
curve(dlnorm(x,1,1),0,5,ylab="density",lwd=2,col=2,add=T,lty=2)
legend(x = "topright",
  legend = c(expression(paste("(",
    paste(mu),",",paste(sigma),") = (0, 1)")),expression(paste("(",
    paste(mu),",",paste(sigma),") = (1, 1)"))),
  lty = c(1, 2),col = c(2, 2),lwd = 2)     
```

## Moments of lognormal

If $Y=ln(X)$ has a normal distribution with mean $\mu$ and standard deviation $\sigma$, then Theorem 6.7 of WMM states:

- $E(X)=\exp(\mu+\sigma^2/2)$
- $Var(X)=\exp(2\mu+\sigma^2)(\exp(\sigma^2)-1)$

If we are interested in relative variation, it is common to look at the coefficient of variation

- $CV(X)=\frac{\sigma}{\mu}$

if e.g. CV=0.05 then 95% of our measurements are within

- $\mu\pm 2\sigma=\mu\pm 2*0.05\mu=\mu(1\pm 0.1)$

i.e. most observations are within 10% of the mean.

## CV of Lognormal

If $Y=ln(X)$ has a normal distribution with mean $\mu$ and standard deviation $\sigma$, we calculate CV for $X$ as

- $CV(X)=\frac{E(X)}{\sqrt{Var(X)}}=\sqrt{\exp(\sigma^2)-1}$

In Peter's data we estimated the variance of the log error to $`r fmt_tex(sigma2a)`$, which means that the estimated CV of the capacity measurement is


- CV = $\sqrt{\exp\left (`r fmt_tex(sigma2a)` \right ) -1} = `r round(100*(sqrt(exp(sigma2a) - 1)), 2)`$%.

i.e., **if** we correct for the systematic error of the meter, then our measurements are extremely precise.

## Linear calibration

In our previous analysis, we assumed, that the systematic error on the meter did not depend on nominal value.

To check this assumption consider the model

- $Y=\mbox{ln(measuredValue)}$ is a linear model of $x= \mbox{ln(nominalValue)}$
- $Y=\alpha+\beta x+\varepsilon$

where we have previously assumed slope($\beta$) equal to 1.

## Linear calibration fit

```{r}
fit <- lm(log(capacity) ~ log(nomval), data = capDat)
summary(fit)
```
The slope is more than close to 1. But is actually extremely significantly different from 1 (tvalue=3776.74 >>>> 3). 

Clearly, it is a bit dubious to assume a linear relationship, as we only have 3 nominal values. 

## Calibrated values

If we stick to the linear calibration model, it is sensible to correct our measured errors according to the calibration of the meter:

- $$\mbox{measuredError}=\alpha+\beta *\mbox{correctError}$$
- $$\mbox{correctError}=(\mbox{measuredError}-\alpha)/\beta$$
  
```{r}
ab = coef(fit)
ab
capDat$lnError_c = (capDat$lnError - ab[1])/ab[2]
```
  
```{r echo=FALSE}
save(ab, file = "ab.RData")
```

## Calibrated data

```{r}
head(capDat)
```

The calibrated data now shows that the production error on component s_1_nF47 is in the vicinity of 0.2%. Well below the tolerance 1%.
