# Theoretical Probability Distributions {#prob}

In the previous chapter we focused our considerations entirely on the sample that we have in front of us. However, as we mentioned earlier, the purpose of collecting a sample usually is to infer characterics of the underlying study population. In order to do so, we often use theoretical probability distributions that can be viewed as models that characterize study populations and that will help us when it comes to statistical inference later on. 

## About Distributions {#aboutDist}

Let's go back to Figure \@ref(fig:3-part-symm) for a second. It depicts the distribution of Martian height within our sample. To be more specific, it depits the _empirical frequency distribution_ of the variable height for the 10 Martians that we collected at site I. Another empirical frequency distribution of the same variable is shown in Figure \@ref(fig:MarHeightHisto). Here the Martians that we collected at site II are also included.  We see that it looks slightly different but has the same basic features (it's roughly symmetrical about the mean).^["The term empirical was originally used to refer to certain ancient Greek practitioners of medicine who rejected adherence to the dogmatic doctrines of the day, preferring instead to rely on the observation of phenomena as perceived in experience." Wikipedia, The Free Encyclopedia, s.v. "Empirical Research," (accessed January 07, 2018), `https://en.wikipedia.org/wiki/
Empirical\_Research`.]

Unfortunately, every sample is likely to result in a different empirical frequency distribution and probably none of them will exactly resemble the distribution of height in the entire Martian population. This is a shame because -- in general -- knowing the distribution of your variable of interest in the population you're studying would provide answers to all questions that you may want to investigate!

To make things more complicated, even if we did have the chance to observe the entire study population (which we usually don't), describing all features of the distribution (all bumps and other quirks) would probably be very cumbersome. So we need a way to simplify things and come up with a model that describes what we imagine the distribution of the variable to look like in the study population. 

Models like this are referred to as _theoretical probability distributions_. Their general shape is governed by a mathematical formula that includes parameters that will determine features like location and spread. There are different prototype theoretical probability distributions, and which one you choose depends on the nature of the variable that you're studying.

In this chapter, we'll cover the Binomial distribution that is a model for binary variables and the Normal distribution ^[Later on we'll also cover the $t$--, the $F$-- and the $\chi^2$--distributions in this workshop.] that is the most common model for symmetrical continuous variables.

Before we start, let's recapitulate what we've learned so far:

**Empirical frequency distributions** are distributions of a variable that we observe in a given sample.

**Rheoretical probability distributions** are mathematical models for the distribution of a variable in the study population. By picking a type of theoretical probability distribution for a variable we fix the general shape of distribution. This is like having a model of how we think the world works.

```{exercise, name="Absolute versus relative"}
Empirical frequency distributions can either be reported in terms of absolute or relative frequencies. Which of the two would be more compatible with theoretical probability distributions?
```


```{exercise, name="Empiricial versus Theoretical"}
Can you think of a property that empiricial (relative) frequency distributions and theoretical probability distributions have in common?
```

Let's explore some common distibutions:

## The Binomial Distribution {#sec:the_binomial_distribution}

As mentioned above, the Binomial distribution is a model for the distribution of a binary variable. To explain things in this context, we'll start with the most used and maybe also most boring example in statistics, i.e. heads or tails when tossing a coin a certain number of times -- sorry about this example, bear with us for a minute!

\begin{equation} 
  f\left(k\right) = \binom{n}{k} p^k\left(1-p\right)^{n-k}
  (\#eq:binom)
\end{equation}

The equation is presented in \@ref(eq:binom).

We said that theoretical probability distributions can be seen as models of how we imagine the world to be. When we are tossing a coin, a natural view of the world's status would be that the coin is fair, i.e. for every toss the probability for head is 0.5 -- same as for tails. With this assumption, we can come up with the distribution for the number of heads that we expect when we toss the coin more than once, let's say we toss it twice: the probability of getting two heads, $p(\mbox{2 heads})$, is $0.5\times0.5 = 0.25$, which is also the case for two tails, $p(\mbox{0 heads})$. There are two ways of getting $p(\mbox{1 head})$, either the first toss is heads or the second one is, thus $p(\mbox{1 head}) = (0.5\times0.5) + (0.5\times0.5) = 0.25 + 0.25 = 0.5$. These probabilities are plotted on the left--hand side Figure \@ref(fig:cointoss05).

Once you understand the above procedure, it's easy to extend the consideration to more than two tosses, _cf._ the middle and right--hand side of Figure \@ref(fig:cointoss05) that show all possible outcomes and the associated probabilities for 5 and 10 tosses of a fair coin, respectively.

```{r cointoss05, fig.cap="Changing n in binomial distributuions, from left to right, B(p = 0.5, n = 2), B(0.5, 5), B(0.5, 10).", fig.width = large_width, fig.height = 2.5, echo = FALSE}

source("./data/PlotFuns.R")
p1 <- PlotBinom(0.5, 2) +
  ggtitle("2 coin tosses") +
  scale_x_discrete("Number of heads", expand = c(0,0)) +
  # theme_classic(24) +
  coord_fixed(2) +
  theme(rect = element_blank(),
        plot.background = element_blank(),
        panel.background = element_blank(),
        axis.text = element_text(color = "black"))
# ggsave("binom_0-5_2.pdf",p1, height = 6, width = 8, units = "in")

p2 <- PlotBinom(0.5, 5) +
  ggtitle("5 coin tosses") +
  scale_x_discrete("Number of heads", expand = c(0,0)) +
  # theme_classic(24) +
  coord_fixed(2) +
  theme(rect = element_blank(),
        plot.background = element_blank(),
        panel.background = element_blank(),
        axis.text = element_text(color = "black"))
# ggsave("binom_0-5_5.pdf",p2, height = 6, width = 8, units = "in")

p3 <- PlotBinom(0.5, 10) +
  ggtitle("10 coin tosses") +
  scale_x_discrete("Number of heads", expand = c(0,0)) +
  # theme_classic(24) +
  coord_fixed(2) +
  theme(rect = element_blank(),
        plot.background = element_blank(),
        panel.background = element_blank(),
        axis.text = element_text(color = "black"))
# ggsave("binom_0-5_10.pdf",p3, height = 6, width = 8, units = "in")

ggarrange(p1, p2, p3, 
          ncol = 3, nrow = 1, align = "h")

```

One very cool feature of theoretical probability functions is that we can use them to probe what we image the world to look like with reality! In the coin toss example, this would mean we want to know if the assumption of a fair coin is justified. And we can do this by actually conducting the experiment and tossing the coin!

```{exercise, name="Exploring the Binomal distribution"}
Imagine you throw the coin twice and get one head -- is that fishy? What about getting only one head when throwing the coin 5 or 10 times, respectively?
```


At some point, we would get very suspicious about our fair coin assumption! Maybe someone tweaked it! Maybe the probability to get a head everytime you toss the coin that you are examining is actually only 0.1! Figure \@ref(fig:cointoss01) shows the Binomial distribution for the number of heads that you get with this assumption for 2, 5, and 10 tosses. They seem to be much more in line with the outcome of our experiments!

```{r cointoss01, fig.cap="Changing p and n in binomial distributuions, from left to right, Bin(p = 0.1, n = 2), Bin(0.1, 5), Bin(0.1, 10).", fig.width = large_width, fig.height = 2.5, echo = FALSE}

p1 <- PlotBinom(0.1, 2) +
  ggtitle("2 coin tosses") +
  scale_x_discrete("Number of heads", expand = c(0,0)) +
  # theme_classic(24) +
  coord_fixed(2) +
  theme(rect = element_blank(),
        plot.background = element_blank(),
        panel.background = element_blank(),
        axis.text = element_text(color = "black"))
# ggsave("binom_0-1_2.pdf",p1, height = 6, width = 8, units = "in")

p2 <- PlotBinom(0.1, 5) +
  ggtitle("5 coin tosses") +
  scale_x_discrete("Number of heads", expand = c(0,0)) +
  # theme_classic(24) +
  coord_fixed(2) +
  theme(rect = element_blank(),
        plot.background = element_blank(),
        panel.background = element_blank(),
        axis.text = element_text(color = "black"))
# ggsave("binom_0-1_5.pdf",p2, height = 6, width = 8, units = "in")

p3 <- PlotBinom(0.1, 10) +
  ggtitle("10 coin tosses") +
  scale_x_discrete("Number of heads", expand = c(0,0)) +
  # theme_classic(24) +
  coord_fixed(2) +
  theme(rect = element_blank(),
        plot.background = element_blank(),
        panel.background = element_blank(),
        axis.text = element_text(color = "black"))
# ggsave("binom_0-1_10.pdf",p3, height = 6, width = 8, units = "in")


ggarrange(p1, p2, p3, 
          ncol = 3, nrow = 1, align = "h")

```


```{exercise, name="Exploring the Binomal distribution"}
If the probability of heads in a single toss is 0.1, what's the probability of getting tails?
```


Now here comes the relief! Coin tosses are not the only application for Binomial probability distributions! They can be used for any binary variable that is sampled a certain number of times, e.g. a phenotype that can be present or absent in some brain cells or a disease that can be early--stage or advanced in a cohort of patients. 

The generic way to deal with the Binomial distribution is to call one of the outcomes _the success_^[Don't let the term _success_ confuse you. Which of the two outcomes you choose as the _success_ is arbitrary! As you have only two outcomes, we automatically know that the probability of no success is $1 - p(success)$ in a single trial.] and to assign to the event "success" a certain probability, $p(success)$, that you expect in a single _trial_, where "trial" can mean coin toss, individual cell or patient -- depending on what you're looking at. 

With these definitions, we can make the general statement that the Binomial probability distribution is goverened by two parameters:  

**The number of trials $n$** that determines the potential number of successes which can be any integer from 0 to $n$ (look at the $x$--axes in Figure \@ref(fig:cointoss05) and Figure \@ref(fig:cointoss01)).

**The probability of success $p$** that governs the shape of the distribution.

To specify a particular Binomial distribution with these two parameters, we write Bin($n$, $p$). 

Use the Binomial Distribution App, explore how changing these parameters affects the probability distribution and answer the following questions:

```{exercise, name="Exploring the Binomal distribution"}
What is the probability of finding exactly 9 heads when tossing a fair coin 25 times? What about 16 heads?
```


```{exercise, name="Exploring the Binomal distribution"}
What is the probability of finding 9 or less heads when tossing a fair coin 25 times? What about 16 or more heads?
```

In the Martian dataset, nose colour is a binary variable that can take the two values blue and red. Let's have a look how we can apply the Binomial distribution here!

In our sample of Martians from site I, we observed that `r (table(martian$Nose[martian$Site == "Site I"])/sum(martian$Site == "Site I")*100)[1]`\% have blue noses. If this observation is a good model for the entire Martian population, we would expect the number of blue--nosed Martians to follow a Binomial distribution with $p(blue \, \, nose) = `r (table(martian$Nose[martian$Site == "Site I"])/sum(martian$Site == "Site I"))[1]`$.

As our dataset also contains the $n=10$ Martians that we (independently) collected at site II, we can directly probe this assumption about the state of the world: at site II, we observed `r (table(martian$Nose[martian$Site == "Site II"])/sum(martian$Site == "Site II")*100)[1]` \% blue-nosed Martians. Does this fit with the model?! Not really -- we don't even need to look at the Bin($10$, $0.9$) distribution for this!

Now, consider a situation that is a bit more ambiguous. Imagine we didn't have the site II Martians, so our hypothesis $p(blue \, \, nose) = `r (table(martian$Nose[martian$Site == "Site I"])/sum(martian$Site == "Site I"))[1]`$ still holds next time we go to Mars. This time we collect a sample of 100 Martians. The Binomial distribution for the number of blue--nosed Martians in this sample that is in line with our hypothesis is Bin $(100,  `r (table(martian$Nose[martian$Site == "Site I"])/sum(martian$Site == "Site I"))[1]`)$. It is depicted below. According to this distribution, we most likely expect to observe `r (table(martian$Nose[martian$Site == "Site I"])/sum(martian$Site == "Site I"))[1]*100` blue-nosed Martians, i.e. `r (table(martian$Nose[martian$Site == "Site I"])/sum(martian$Site == "Site I")*100)[1]` \% of the sample. It turns out, however, that instead we observe 85. How likely is this?


```{r echo = FALSE}

prob <- (table(martian$Nose[martian$Site == "Site I"])/sum(martian$Site == "Site I"))[1]
size <- 100
obs <- 85

val <- dbinom(0:size, size, prob)
df.raw <- data.frame(variable = seq(0, size, 1), values = val)
df.martian <- as.data.frame(df.raw)
names(df.martian) <- c("variable", "values")

point.prob <- df.martian$values[df.martian$variable == obs]

lower.bound <- max(df.martian$variable[df.martian$values <= point.prob & df.martian$variable < obs])

df.martian$col <- NULL
df.martian$col[df.martian$values <= point.prob] <- MidRed # Red
df.martian$col[df.martian$values > point.prob] <- MidYellow # Yellow

# Draw basic yellow probability plot:
martian.prob <-
  ggplot(df.martian, aes(x = variable, y = values, fill = col)) +
  scale_x_continuous(paste("Binomial Probability Distribution of ", size," trials\nwith a ", prob, " probability of success", sep = ""),
                     #labels = categories,
                     limit = c(79, 100),
                     labels = seq(80,100,5),
                     breaks = seq(80,100,5),
                     expand = c(0,0)) +
  scale_y_continuous("Probability",
                     limit = c(0, 0.14),
                     labels = seq(0, 0.14, 0.01),
                     breaks = seq(0, 0.14, 0.01),
                     expand = c(0,0))  +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )

# martian.prob +
#   geom_bar(stat = "identity", fill = MidYellow)

outcome <- dbinom(obs, size, prob)

# Now draw lines showing the value we are interested in:
martian.prob +
  geom_bar(stat = "identity", fill = MidYellow) +
  geom_vline(xintercept = obs) +
  geom_hline(yintercept = df.martian$values[df.martian$variable == obs])


```


Use the Binomial Distribution App to answer the following questions: 

```{exercise, name="Exploring the Binomal distribution"}
Imagine that we knew that the true $p(blue \,\, nose)$ of getting a blue-nosed Martian is $0.6$, what is the probability of observing $25$ blue-nosed Martians in a sample of 30 Martians?
```

```{exercise, name="Exploring the Binomal distribution"}
What is the probability of observing $25$ _or more_ blue-nosed Martians if the $p(blue \,\, nose)$ is $0.6$?
```

In the next section we'll look at probably the most famous distrbution in all of statistics, the Normal distribution.

## The Normal Distribution {#normalDist}

The Binomial distribution serves as a model for binary variables. When we are dealing with continuous variables instead, a very common model is provided by the _Normal distribution_^[The equation for the Normal distribution was developed by Karl Friedrich Gauss in 1809.]. It has the following very simple and famous equation (which you don't need to memorise):

\begin{equation} 
  f(x) = \frac{1}{\sigma\sqrt{2\pi} } \; \mbox{e}^{ -\frac{(x-\mu)^2}{2\sigma^2} }.
  (\#eq:norm)
\end{equation} 

To understand this formula, let's look at the components of the right--hand side of equation \@ref(eq:norm). Apart from $x$ (which is just the argument of the function), it contains the number 2 (nothing special)^[In Statistics the convention is that Roman letters, e.g. $\bar{x}$ and $s$, are used to refer to the sample, Greek letters, e.g. $\mu$ and $\sigma$, refer to the population.], the natural constants $\mbox{e}$ and $\pi$ (which are just numbers), and two Greek letters $\mu$ and $\sigma$ that don't have an immediately obvious interpretation.

Just like $n$ and $p$ are the parameters of the Binomial distribution, Bin($n$, $p$), that control the appearance of the probability function (cf. Figures \@ref(fig:cointoss05) and \@ref(fig:cointoss01)), $\mu$ and $\sigma$ are the parameters of the Normal distribution that govern the shape of the associated curve. As a short hand, we can just write $N(\mu, \sigma^2)$. A generic example of a Normal curve is shown in Figure \@ref(fig:genericnormal). 

<!-- , fig.cap = "The generic shape of the Normal distribution. The symmetrical curve is centered on $\mu$, the solid red line, and its width is defined by $\sigma$. The two infliction points are situated at $\mu\pm1\sigma$." -->

```{r genericnormal, fig.cap = "The generic normal curve. The shape of the distribution depends on the two parameters.", echo = FALSE, fig.width = large_width}

  muS <- 0
  sigmaS <- 1
  xRangeS <- 4 # Set the upper limit of the x-axis

  xS <- seq(muS-(xRangeS*2), muS+(xRangeS*2), length=5000) # Define an x series based on inputed mu and xRange
  hxS <- dnorm(xS, mean = muS, sd = sigmaS) # Define the density of the normal curve at each x position
  hstdS <- dnorm(xS) # Define the density of the standard normal curve at each x position

test.dfS <- data.frame(value = xS, dist = hstdS)
new.dfS <- data.frame(value = seq(-1,1,2), dist =c(dnorm(-1),dnorm(1)))


# Base plot:




# Base plot:
q0 <- ggplot(test.dfS, aes (x= value, y = dist)) +
  geom_line() +
  scale_y_continuous("f(x)", expand = c(0,0)) +
  coord_cartesian(xlim = c(-xRangeS, xRangeS), ylim = c(0, 0.5)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
#         axis.text.y = element_blank(),
#         axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )

q <- q0 +
    scale_x_continuous("x", breaks = seq(-xRangeS, xRangeS, 1), labels = c(
    expression(paste(mu, "-4", sigma)), #"$-4\\sigma$", # -4
    expression(paste(mu, "-3", sigma)), #"$-3\\sigma$", # -3
    expression(paste(mu, "-2", sigma)), #"$-2\\sigma$", # -2
    expression(paste(mu, "-1", sigma)), #"$-1\\sigma$", # -1
    expression(mu), #"$\\mu$", # 0
    expression(paste(mu, "+1", sigma)), #"$1\\sigma$", # -1
    expression(paste(mu, "+2", sigma)), #"$2\\sigma$", # -2
    expression(paste(mu, "+3", sigma)), #"$3\\sigma$", # -3
    expression(paste(mu, "+4", sigma)) #"$4\\sigma$" # -4
  ), expand = c(0,0))


q + geom_vline(xintercept = muS, colour = MidRed) +
  # geom_vline(xintercept = c(muS+sigmaS, muS-sigmaS), colour = MidRed, linetype = 2) +
  geom_point(data = new.dfS, size = 4, shape = 1, col = MidRed)


```

Let's take note of the following properties of the curve and the role of the two parameters. Refer to Figure \@ref(fig:changeParam).

### The parameter $\mu$

The parameter $\mu$ dictates the location of the curve on the $x$--axis: the curve peaks at $\mu$ and is decreasing in a symmetrical and monotone way as it moves away from $\mu$. This is further demonstrated in Figure \@ref(fig:conSvarU) which shows two Normal distributions with different values of $\mu$ ($\mu=0$ and $\mu=-3$) while keeping $\sigma = 1$ constant. As we will see, $\mu$ can be interpreted as the (usually unknown) population mean of the continuous variable that we model using the Normal distribution.

```{r changeParam, fig.cap="The shape of the normal distribution is determined by its two parameters. Although they look quite different, these are all Normal distributions. ", echo = FALSE, fig.width = large_width}

### This figure draws four normal curves superimposed on one another.
# Six curves:
# mu <- c(0, 0, 0, -4, -4, -4)
# sigma <- c(0.50, 1.00, 2.00, 0.50, 1.00, 2.00)
# myCol <- c(LightYellow, MidYellow, DarkYellow,
#            LightRed, MidRed, DarkRed)

# Four curves:
mu <- c(0, 0, 0, -3)
sigma <- c(0.50, 1.00, 2.00, 1.00)
myCol <- c(LightBlue, MidBlue, DarkBlue, MidRed)

xRange <- 5 # Set the upper limit of the x-axis

makeNorm <- function(mu, sigma, range) {
  x <- seq(range*-2, range*2, length=5000)
  hx <- dnorm(x, mean = mu, sd = sigma)
  df <- data.frame(hx)
  names(df)[1] <- paste0(mu,"_",sigma)
  return(df)
}


test.df <- mapply(makeNorm, mu, sigma, xRange)
test.df <- do.call(cbind.data.frame, test.df)
test.df <- reshape2::melt(test.df, measure.vars = 1:ncol(test.df), value.name = "dist")
test.df$value <- rep(seq(xRange*-2, xRange*2, length=5000), length(mu))
# test.df <- mapply(makeNorm, 0, 1, 4)
# test.df <- do.call("rbind", test.df)
# test.df <- do.call("cbind", test.df)
# test.df <- ldply(test.df, data.frame)

test.df <- cbind(test.df,
      reshape2::colsplit(test.df$variable, "_", c("mu", "sigma")))

# test.df$ribbon <- 0

test.df <- plyr::ddply(test.df, "variable", transform,
      ribbon = ifelse(value > mu-1.96*sigma & value < mu+1.96*sigma, 1, 0))

sigmapoints.df <- data.frame(value = c(mu + sigma, mu - sigma), variable = rep(paste0(mu,"_",sigma),2))
test.df2 <- test.df
test.df2$value <- round(test.df2$value, 2)
sigmapoints.df <- merge(sigmapoints.df, test.df2, by=c("variable","value"))
sigmapoints.df <- plyr::ddply(sigmapoints.df, c("variable", "value"), summarise, dist = mean(dist))


levels(test.df$variable) <- sub("_", "", levels(test.df$variable))
levels(sigmapoints.df$variable) <- sub("_", "", levels(sigmapoints.df$variable))
# myLabels <-

# Base plot:
gg <- ggplot(test.df, aes(x= value, y = dist, col = variable, fill = variable)) +
  geom_line() +
  # geom_vline(aes(col = variable), xintercept = mu) +
  # geom_vline(xintercept = c(mu+sigma, mu-sigma), colour = MidRed, linetype = 2) +
  scale_x_continuous("x", breaks = seq(-xRange-2, xRange+2, 1), expand = c(0,0)) +
  scale_y_continuous(expression(paste("f(x | ", mu, ", ", sigma, ")")),
                     # "$f(x \\; | \\; \\mu, \\sigma)$", 
                     expand = c(0,0)) +
  # scale_color_brewer("$f(x \\; | \\; \\mu, \\sigma) =$") +
#     scale_color_manual("Parameters",
#                        labels = paste("$f(x \\; | \\; \\mu = ", mu, ",  \\sigma = ", sigma, ")$", sep = ""),
#                        values = myCol) +
  
    scale_color_manual(expression(paste("f(x | ", mu, ", ", sigma, ")")), #"$f(x \\; | \\; \\mu, \\sigma) =$",
                       labels = paste0("f(x | ", mu, ",  ", sigma, ")"),
                       values = myCol) +
  
  
    # scale_color_manual("$f(x \\; | \\; \\mu, \\sigma) =$",
    #                    labels = paste("$f(x \\; | \\; ", mu, ",  ", sigma, ")$", sep = ""),
    #                    values = myCol) +
      # scale_fill_manual("$f(x \\; | \\; \\mu, \\sigma) =$",
      #                  labels = paste("$f(x \\; | \\; ", mu, ",  ", sigma, ")$", sep = ""),
      #                  values = myCol) +

#  For Six curves:
#     coord_cartesian(xlim = c(-xRange-2, xRange-2), ylim = c(0, 1.0)) +
#  For Four curves:
  coord_cartesian(xlim = c(-xRange-1, xRange-2), ylim = c(0, 1.0)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
        legend.position = c(0.9, 0.8),
#         axis.text.y = element_blank(),
#         axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )

gg # + geom_point(aes(x= value, y = dist, col = variable, fill = variable), data = sigmapoints.df, size = 4, shape = 1, show.legend = FALSE)


# Make data frame with 1.96 areas:
# xPoints <- c(1.96*sigma, -1.96*sigma)
#
# sigmapoints.df <- data.frame(
#   dist = mapply(dnorm, xPoints, mu, sigma),
#   value = xPoints,
#   variable = rep(paste0(mu,"_",sigma),2)
# )
# levels(sigmapoints.df$variable) <- sub("_", "", levels(sigmapoints.df$variable))


# paste("$f(x \\; | \\; ", mu, ",  ", sigma, ")$", sep = "")

# Tag region for ribbon

gg <- gg + geom_ribbon(data=subset(test.df, ribbon == 1),
                  aes(x=value,ymax=dist, fill = variable), ymin=0, colour=NA) +
  facet_wrap(~ variable, ncol = 1) +
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        legend.position = c(0.2, 0.9))


```

The effects of changing $\mu$ on the Normal curve. As $\mu$ changes, so does the center of the distribution. The width of both curves is identical, since their $\sigma$ values are the same. The infliction points, corresponding to $\mu\pm1\sigma$, are shown as dots on each curve.

```{r conSvarU, fig.cap="Changing mu.", echo = FALSE, fig.width = large_width}

### This figure draws two normal curves according to changes in mean.

# three curves:
mu <- c(0, -3)
sigma <- c(1.00, 1.00)
myCol <- c(MidBlue, MidRed)
# Colour according to sigma:
# myCol <- c(DarkBlue, DarkYellow, DarkRed,
#                     MidYellow)

xRange <- 5 # Set the upper limit of the x-axis

makeNorm <- function(mu, sigma, range) {
  x <- seq(range*-2, range*2, length=5000)
  hx <- dnorm(x, mean = mu, sd = sigma)
  df <- data.frame(hx)
  names(df)[1] <- paste0(mu,"_",sigma)
  return(df)
}

test.df <- mapply(makeNorm, mu, sigma, xRange)
test.df <- do.call(cbind.data.frame, test.df)
test.df <- reshape2::melt(test.df, measure.vars = 1:ncol(test.df), value.name = "dist")
test.df$value <- rep(seq(xRange*-2, xRange*2, length=5000), length(mu))

test.df <- cbind(test.df,
      reshape2::colsplit(test.df$variable, "_", c("mu", "sigma")))

# test.df$ribbon <- 0

test.df <- plyr::ddply(test.df, "variable", transform,
      ribbon = ifelse(value > mu-1.96*sigma & value < mu+1.96*sigma, 1, 0))

sigmapoints.df <- data.frame(value = c(mu + sigma, mu - sigma), variable = rep(paste0(mu,"_",sigma),2))
test.df2 <- test.df
test.df2$value <- round(test.df2$value, 2)
sigmapoints.df <- merge(sigmapoints.df, test.df2, by=c("variable","value"))
sigmapoints.df <- plyr::ddply(sigmapoints.df, c("variable", "value"), summarise, dist = mean(dist))


levels(test.df$variable) <- sub("_", "", levels(test.df$variable))
levels(sigmapoints.df$variable) <- sub("_", "", levels(sigmapoints.df$variable))

# Base plot:
gg <- ggplot(test.df, aes (x= value, y = dist, col = variable, fill = variable)) +
  geom_line() +
  # geom_vline(aes(col = variable), xintercept = mu) +
  # geom_vline(xintercept = c(mu+sigma, mu-sigma), colour = MidRed, linetype = 2) +
  scale_x_continuous("x", breaks = seq(-xRange-2, xRange+2, 1), expand = c(0,0)) +
  scale_y_continuous(expression(paste("f(x | ", mu, ", ", sigma, ")")),
                     # "$f(x \\; | \\; \\mu, \\sigma)$", 
                     expand = c(0,0)) +
  # scale_color_brewer("$f(x \\; | \\; \\mu, \\sigma) =$") +
#     scale_color_manual("Parameters",
#                        labels = paste("$f(x \\; | \\; \\mu = ", mu, ",  \\sigma = ", sigma, ")$", sep = ""),
#                        values = myCol) +
  
      scale_color_manual(expression(paste("f(x | ", mu, ", ", sigma, ")")), #"$f(x \\; | \\; \\mu, \\sigma) =$",
                       labels = paste0("f(x | ", mu, ",  ", sigma, ")"),
                       values = myCol) +

    # scale_color_manual("$f(x \\; | \\; \\mu, \\sigma) =$",
    #                    labels = paste("$f(x \\; | \\; ", mu, ",  ", sigma, ")$", sep = ""),
    #                    values = myCol) +
    #   scale_fill_manual("$f(x \\; | \\; \\mu, \\sigma) =$",
    #                    labels = paste("$f(x \\; | \\; ", mu, ",  ", sigma, ")$", sep = ""),
    #                    values = myCol) +

#  For Six curves:
#     coord_cartesian(xlim = c(-xRange-2, xRange-2), ylim = c(0, 1.0)) +
#  For Four curves:
  coord_cartesian(xlim = c(-xRange-1, xRange-1), ylim = c(0, 1.0)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
        legend.position = c(0.9, 0.8),
#         axis.text.y = element_blank(),
#         axis.text.x = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )

gg # + geom_point(data = sigmapoints.df, size = 4, shape = 1, show.legend = FALSE)

# gg <- gg + geom_ribbon(data=subset(test.df, ribbon == 1),
#                   aes(x=value,ymax=dist, fill = variable), ymin=0, colour=NA) +
#   facet_wrap(~ variable, ncol = 1) +
#   theme(strip.background = element_blank(),
#         strip.text = element_blank(),
#         legend.position = c(0.2, 0.9))

```

### The parameter $\sigma$

The parameter $\sigma$ determines the width of the curve: the two infliction points are situated at $\mu\pm1\sigma$. This is further demonstrated in Figure \@ref(fig:varSconU) which shows three Normal distributions with different levels of $\sigma$ ($\sigma=0.5$, $\sigma=1$, $\sigma=2$) while keeping the center constant at $\mu =0$. We see that the curve becomes higher and narrower as $\sigma$ decreases and, conversly, becomes flatter and broader as $\sigma$ increases. Thus the Normal curve is like an accordion, it can be anything between thin & tall and wide & short. As we will see, $\sigma$ can be interpreted as the (usually unknown) population standard deviation of the continuous variable that we model using the Normal distribution.

fig cap: The effects of changing $\sigma$ on the Normal curve. As $\sigma$ increases (from light to dark blue curves), so does the spread of the distribution. The curves become wider and flatter whereas the height decreases. The infliction points, corresponding to $\mu\pm1\sigma$, are shown as dots on each curve.

```{r fig:varSconU, fig.cap="Changing the variance.", echo = FALSE, fig.width = large_width}


### This figure draws three normal curves according to different sigmas, superimposed on one another.

# Three curves:
mu <- c(0, 0, 0)
sigma <- c(0.5, 1.0, 2.0)
myCol <- c(LightBlue, MidBlue, DarkBlue)
xRange <- 5 # Set the upper limit of the x-axis

test.df <- mapply(makeNorm, mu, sigma, xRange)
test.df <- do.call(cbind.data.frame, test.df)
test.df <- reshape2::melt(test.df, measure.vars = 1:ncol(test.df), value.name = "dist")
test.df$value <- rep(seq(xRange*-2, xRange*2, length=5000), length(mu))

test.df <- cbind(test.df,
      reshape2::colsplit(test.df$variable, "_", c("mu", "sigma")))

test.df <- plyr::ddply(test.df, "variable", transform,
      ribbon = ifelse(value > mu-1.96*sigma & value < mu+1.96*sigma, 1, 0))

test.df$sigma <- factor(test.df$sigma)

# levels(test.df$sigma) <- paste("sigma =", levels(test.df$sigma))
cutoff.df <- data.frame(sigma = levels(test.df$sigma),
                        lower = mu - qnorm(0.5 + 0.95/2, mu, sigma),
                        upper = mu + qnorm(0.5 + 0.95/2, mu, sigma))
#                         ,
#                         value = 2,
#                         dist = 0.5,
#                         text = c("A", "B", "C"))
#  paste("sigma =", levels(test.df$sigma))
# cutoff.df$text <- as.character(cutoff.df$text)

sigmapoints.df <- data.frame(value = c(mu + sigma, mu - sigma), variable = rep(paste0(mu,"_",sigma),2))
test.df2 <- test.df
test.df2$value <- round(test.df2$value, 2)
sigmapoints.df <- merge(sigmapoints.df, test.df2, by=c("variable","value"))
sigmapoints.df <- plyr::ddply(sigmapoints.df, c("variable", "value"), summarise, dist = mean(dist))

levels(test.df$variable) <- sub("_", "", levels(test.df$variable))
levels(sigmapoints.df$variable) <- sub("_", "", levels(sigmapoints.df$variable))

# Base plot:
gg <- ggplot(test.df, aes (x= value, y = dist, col = variable, fill = variable)) +
  geom_line() +
  scale_x_continuous("x", breaks = seq(-xRange-2, xRange+2, 1), expand = c(0,0)) +
  scale_y_continuous(expression(paste("f(x | ", mu, ", ", sigma, ")")),
                     # "$f(x \\; | \\; \\mu, \\sigma)$", 
                     expand = c(0,0)) +

      scale_color_manual(expression(paste("f(x | ", mu, ", ", sigma, ")")), #"$f(x \\; | \\; \\mu, \\sigma) =$",
                       labels = paste0("f(x | ", mu, ",  ", sigma, ")"),
                       values = myCol) +

  
  #   scale_color_manual("$f(x \\; | \\; \\mu, \\sigma) =$",
  #                      labels = paste("$f(x \\; | \\; ", mu, ",  ", sigma, ")$", sep = ""),
  #                      values = myCol) +
  # scale_fill_manual("$f(x \\; | \\; \\mu, \\sigma) =$",
  #                      labels = paste("$f(x \\; | \\; ", mu, ",  ", sigma, ")$", sep = ""),
  #                      values = myCol) +
  coord_cartesian(xlim = c(-xRange-1, xRange-1), ylim = c(0, 1.0)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
        legend.position = c(0.9, 0.8),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )

gg# + geom_point(data = sigmapoints.df, size = 4, shape = 1, show.legend = FALSE) #+
  # geom_text(data = ann_text, aes(label = lab)) +


## Annotating text in each facet:
# p <- ggplot(mtcars, aes(mpg, wt)) + geom_point() +
#   facet_grid(. ~ cyl)
# p
#
# ann_text <- data.frame(mpg = c(15,20),
#                        wt = c(4,5),
#                        lab = c("text1","text2"),
#                        cyl = c(6,8))
# ann_text$cyl
#
# p + geom_text(data = ann_text, aes(label =lab))

# test.df$sigma <- as.numeric(as.character(test.df$sigma))

# str(test.df)


# ann_text <- plyr::ddply(test.df, "sigma", summarise,
#                   variable = as.character(variable),
#                   dist = dist,
#                   value = value,
#                   mu = mu,
#                   sigma = as.character(sigma),
#                   ribbon = 0,
#                   text = as.character(sigma))

# names(test.df)
# variable
# dist
# value
# mu
# sigma
# ribbon

ann_text <- data.frame(variable = levels(test.df$variable),
                  dist = 0.75,
                  value = -5,
                  mu = mu,
                  sigma = as.character(sigma),
                  ribbon = 1,
                  text = paste0("$\\sigma =$", as.character(sigma)))

# str(ann_text)

# ann_text$lab <- as.character(ann_text$lab)
# str(ann_text)
# str(cutoff.df)

# str(test.df)


test.df <- plyr::ddply(test.df, "sigma", transform,
                 text = as.character(sigma))



# names(test.df)

gg <- gg + geom_ribbon(data=subset(test.df, ribbon == 1),
                  aes(x = value,
                      ymax = dist,
                      fill = variable), ymin=0, colour=NA, show.legend = FALSE) +
  geom_vline(data = cutoff.df, aes(xintercept = upper), colour = MidRed, linetype = 2) +
  geom_vline(data = cutoff.df, aes(xintercept = lower), colour = MidRed, linetype = 2) +
  geom_text(data = ann_text, aes(label = text), x = -5, y = 0.75, col = "black", hjust = 0, size = 3.5) +
  # annotate(geom="text", x=0, y=0.5, label="A",
              # color="red", data = cutoff.df) +
  # annotate("text", data = cutoff.df, aes(label = label)) +
  facet_wrap(~ sigma, ncol = 1) +
  # facet_grid(sigma ~ ., labeller = label_bquote(sigma == .(x))) +
  # guides(fill = "none") +
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        legend.position = "none")

```

## Normal probabilites -- the area under the Normal distribution

The Normal distribution can be used as mathematical description of how we imagine a continous variable $X$ (e.g.\,$X$ could be Martian height) to behave in the study population. It therefore has the same function as the Binomial distribution for binary variables, but there is a big difference! From the Binomial distribution, the expected probability of a certain outcome (e.g. finding 85 blue--nosed Martians in a sample of 100 Martians when $p(blue \,\, nose) = 0.9$) could be directly obtained but with a continous distribution this is not possible. In fact, the values on the $y$--axis cannot be interpreted as probabilities!

What _can_ be interpreted as probabilities though is the _area under the curve_. This is because the entire area under the curve of any Normal distribution integrates to one. ^[Integrating to find the area under the curve for continous distributions is equivalent to summing up the individual probabilities of all different outcomes of a discrete variable.] Take for example the red curve in Figure \@ref(fig:conSvarU). While it is not possible to determine the probability for $X=-2.5$, it easy to see that the probability for $X< -3$ is 50\% because the area under the curve to the left of $\mu$ is just as big as the area under the curve on the other side.

For a generic Normal distribution $N(\mu, \sigma)$, there are a number probabilities that are easy to remember. These are the probabilities that $X$ takes a value within certain _central intervals_ around $\mu$, cf. Table \@ref(tab:Normalperc). The integrations that need to be carried out to calculate the probabilities are shown in Figure \@ref(fig:DifferentNormalCurves). Notice that less than 5\% of the distribution is to be found beyond $\mu\pm2\sigma$. 

Table: (\#tab:Normalperc) Probabilities that the value of the variable falls within some central intervals around $\mu$. This holds for \emph{any} combination of $\mu$ \& $\sigma$!

| Percentage of area under the curve | Range                              | Figure |
|------------------------------------|------------------------------------|------------------------------------------|
| 68.26 %                            | [$\mu-1\sigma, \mu+1\sigma$]       | \@ref(fig:DifferentNormalCurves), top    |
| 95.44 %                            | [$\mu-2\sigma, \mu+2\sigma$]       | \@ref(fig:DifferentNormalCurves), top    |
| 99.72 %                            | [$\mu-3\sigma, \mu+3\sigma$]       | \@ref(fig:DifferentNormalCurves), top    |
| 95 %                               | [$\mu-1.96\sigma, \mu+1.96\sigma$] | \@ref(fig:DifferentNormalCurves), bottom |

figcap: Central Intervals of the Normal distribution. _top_: For any Normal distribution defined by $\mu$ and $\sigma$, 68.26\% of the area under the curve will lie in the interval $[\mu-1\sigma, \mu+1\sigma]$, 95.44\% in the interval $[\mu-2\sigma, \mu+2\sigma]$, and 99.72\% in the interval $[-3\sigma, 3\sigma]$. _bottom_: 95\% of the area under the curve  is found within the interval $[\mu-1.96\sigma, \mu+1.96\sigma]$.

```{r DifferentNormalCurves, fig.cap="Central Intervals of the Normal distribution. _top_: For any Normal distribution defined by. _bottom_: 95% of the area under the curve  is found within the interval.", echo = FALSE, fig.width = large_width}

mu <- muS
sigma <- sigmaS
xRange <- xRangeS

test.df <- test.dfS
new.df <- new.dfS

aBlend <- 0.5

p1 <- q + geom_ribbon(data=subset(test.df, value > -1 & value < 1),
                aes(x=value,ymax=dist),ymin=0,fill=MidBlue, alpha = aBlend) +
  # ggtitle("$\\mu\\pm1\\sigma=68.2\\%$") +
geom_vline(xintercept = mu + sigma, colour = MidRed, linetype = 2) +
geom_vline(xintercept = mu - sigma, colour = MidRed, linetype = 2) +
annotate("text", 
         x = sigma+0.05, 
         y = 0.45, 
         label = expression(mu~1*sigma), 
         hjust = 0, 
         size = 3.5) +
annotate("text", x = sigma+0.05, y = 0.42, label = "= 68.26 %", hjust = 0, size = 3.5) +

geom_ribbon(data=subset(test.df, value > -2 & value < 2   ),
                  aes(x=value,ymax=dist),ymin=0,fill=MidBlue, alpha = aBlend) +
                  # ggtitle("$\\mu\\pm2\\sigma=95.4\\%$") +
geom_vline(xintercept = mu + (2*sigma), colour = MidRed, linetype = 4) +
geom_vline(xintercept = mu - (2*sigma), colour = MidRed, linetype = 4) +
annotate("text", 
         x = 2*sigma+0.05, 
         y = 0.45, 
         label = expression(mu~2*sigma), 
         hjust = 0, 
         size = 3.5) +
annotate("text", 
         x = 2*sigma+0.05, 
         y = 0.42, 
         label = "= 95.44 %", 
         hjust = 0, 
         size = 3.5) +

geom_ribbon(data=subset(test.df, value > -3 & value < 3   ),
                  aes(x=value,ymax=dist),ymin=0,fill=MidBlue, alpha = aBlend) +
                  # ggtitle("$\\mu\\pm3\\sigma=99.8\\%$") +
geom_vline(xintercept = mu + (3*sigma), colour = MidRed, linetype = 3) +
geom_vline(xintercept = mu - (3*sigma), colour = MidRed, linetype = 3) +
annotate("text", 
         x = 3*sigma+0.05, 
         y = 0.45, 
         label = expression(mu~3*sigma), 
         hjust = 0, 
         size = 3.5) +
annotate("text", x = 3*sigma+0.05, y = 0.42, label = "= 99.72 %", hjust = 0, size = 3.5)


p2 <- q +
geom_ribbon(data=subset(test.df, value > qnorm(0.5 - 0.95/2) & value < qnorm(0.5 + 0.95/2)   ),
                  aes(x=value,ymax=dist),ymin=0,fill=MidBlue) +
geom_vline(xintercept = mu + qnorm(0.5 + 0.95/2), colour = MidRed, linetype = 2) +
geom_vline(xintercept = mu - qnorm(0.5 + 0.95/2), colour = MidRed, linetype = 2) +
annotate("text", 
         x = 1.96+0.05, 
         y = 0.45, 
         label = expression(mu~1.96*sigma), 
         # label = expression(paste0(mu, pm, "1.96", sigma)), 
         hjust = 0, 
         size = 3.5) +
annotate("text", x = 1.96+0.05, y = 0.42, label = "= 95.0 %", hjust = 0, size = 3.5)


ggarrange(p1, p2, 
          nrow = 2, ncol = 1, align = "v")

```

Let's have a second look at the 95\% central interval [$\mu-1.96\sigma, \mu+1.96\sigma$]. If we apply it to the Normal curves in Figure _fig:varSconU_ that are centered around $\mu =0$ we can appreciate the influence of $\sigma$. As $\sigma$ increases, the curves get broader and so do the 95\% central intervals: [$-0.98, 0.98$] for $\sigma = 0.5$, [$-1.96, 1.96$] for $\sigma = 1$, and [$-3.93, 3.92$] for $\sigma = 2$. If you want to remember one of them, pick the second one as it will reappear a lot later!

Use the Normal Distribution App to explore the parameters and central intervals of the Normal distribution. Also answer the exercises below!

In our Martian dataset, the continous variable that we looked at was height. As the associated empirical frequency distributions are roughly symmetrical around the mean (cf. Figures \@ref(fig:3-part-symm) and \@ref(MarHeightHisto)), it seems OK to use the Normal distribution as model for the variable height in the underlying population (we'll discuss a more formal check for this in a bit, cf. Section \@ref(sec:qqplots)). 

```{r, echo = FALSE}

H.Means <- tapply(martian$Height, martian$Site, mean)
H.SDs <- plyr::round_any(tapply(martian$Height, martian$Site, sd),1)
mystery <- 215

```

```{exercise, name="Exploring the Normal distribution"}
Consider the model that the height of Martians at site I can be described by the Normal distribution of $\mu = 200$ cm and $\sigma = 6$ cm. 95 % of Martians at site I fall into what interval?
```

```{exercise, name="Exploring the Normal distribution"}
Now, consider that the model that the height of Martians at site II can be described by the Normal distribution with a $\mu$ = 210$ cm and $\sigma = 1$ cm. 95 % of Martians at site II fall into what interval?
```

<!-- $\sigma = `r H.SDs[2]`$cm. -->

```{exercise, name="Exploring the Normal distribution"}
Make a sketch of the two Normal distributions on a common scale in the space provided:
```

```{r echo = FALSE, fig.width = large_width}

fake.df <- data.frame(x = 1:10, y = 1:10)

ggplot(fake.df, aes (x = x, y = y)) +
  geom_point() +
#   geom_histogram(aes(y = ..density..), binwidth = 2, fill = MidYellow) +
  scale_y_continuous(name = "f(x)", expand = c(0,0)) +
  scale_x_continuous("Martian Height (cm)", limits = c(180,220), expand = c(0,0)) +
#   stat_function(fun = dnorm,
#                 args = list(mean = mean(martian$Height[martian$Site == "Site I"]),
#                             sd = sd(martian$Height[martian$Site == "Site I"])),
#                 col = DarkRed) +
  theme(plot.background = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.ticks.y = element_blank(),
        axis.line = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_blank()
  )

```

```{exercise, name="Identifying extreme values using the Normal distribution"}
Imagine that we are presented with a 215 cm tall Martian. Draw a vertical line on your plot from the previous exercise.
```

```{exercise, name="Exploring the Notmal distribution"}
Given our assumptions about the distributions of  Martian height at site I and II, which site do you think this Martian came from, and why?
```

As you can imagine this is pretty tedious. For every new population we have to draw a curve or calculate the inverval whitin which we think 95% of the population will exisit. It would be much easier to have a standardized way of doing this... enter Z--scores.

## Z-scores and the Standard Normal Curve {#z_scores}

Notice that in the exercise above, we had to draw two Normal curves and then position our mystery $`r mystery`$ cm tall Martian in this landscape to compare the likelihood of the two scenarios (sites). Luckily there is an easier way to make this comparison by standardising the value of unknown origin ($`r mystery`$ cm in our case) against our two models under consideration. That is the role of Z--scores:

> Z--scores allow us to normalise disparate data onto a common scale!

The Z--score for any given value $X_i$ with respect to a model that assumes a mean value of $\mu$ and a standard deviation of $\sigma$ can be obtained using the following normalisation: ^[In reality we don't know $\mu$ and $\sigma$ so we just use $\bar{x}$ and $s$. We'll come back to this point when we encounter estimation.]

$$Z_i = \frac{X_i - \mu}{\sigma}.$$

We'll encounter this surprisingly simple yet useful equation many times throughout the workshop! Here we'll demonstrate its value by re--addressing the exercise above.

^[We can also understand Z--scores as $$\frac{X_i - \mu}{\sigma} = \frac{\mbox{signal}}{\mbox{noise}}.$$ The signal is the systematic variation we are interested in and the noise is the random error which we cannot avoid. Keep this in mind and we'll get back to it in the context of hypothesis testing.]

```{exercise, name="Calculating Z--scores"}
Consider the population described by $\mu = 200$ cm and $\sigma = 6$ cm, complete table \@ref(tab:Z1).
```

Table: (\#tab:Z1) Calculating Z--scores for the Martian height at site I.

| $x_i (cm)$  | `r H.Means[1] - 3*H.SDs[1]` | `r H.Means[1] - 2*H.SDs[1]` | `r H.Means[1] - H.SDs[1]` | `r H.Means[1]` | `r H.Means[1] + H.SDs[1]` | `r H.Means[1] + 2*H.SDs[1]` | `r H.Means[1] + 3*H.SDs[1]` |
|-------------|-----------------------------|-----------------------------|---------------------------|----------------|---------------------------|-----------------------------|-----------------------------| 
| $x_i - \mu$ |                             |                             |                           |                |                           |                             |                             |  
| $z_i$       |                             |                             |                           |                |                           |                             |                             |

```{exercise, name="Calculating Z--scores"}
Continuing with the population described by $\mu = 210$ cm and $\sigma = 1$ cm, complete table \@ref(tab:Z2):
```


Table: (\#tab:Z2) Calculating Z--scores for the Martian height at a separate site.

| $x_i (cm)$  | `r H.Means[2] - 3*1` | `r H.Means[2] - 2*1` | `r H.Means[2] - 1` | `r H.Means[2]` | `r H.Means[2] + 1` | `r H.Means[2] + 2*1` | `r H.Means[2] + 3*1` |
|-------------|----------------------|----------------------|--------------------|----------------|--------------------|----------------------|----------------------| 
| $x_i - \mu$ |                      |                      |                    |                |                    |                      |                      |  
| $z_i$       |                      |                      |                    |                |                    |                      |                      |  

From the results of these exercises we can understand that Z--scores simply tell us how many standard deviations a value is away from the mean with respect to a given model. If the model is described by the Normal distribution, Z--scores have the very cool feature that they themselves are Normally distributed!! As they are normalised, the corresponding Normal distribution has a mean of 0 and a standard deviation of 1. This all follows logically from the equation for Z--scores and can be written in a shorthand form as:

\begin{equation} 
  Z_i = \frac{X_i - \mu}{\sigma} \sim N(0, 1).
  (\#eq:Z)
\end{equation} 

The distribution $N(0, 1)$ is called the _standard Normal distribution_ and is shown is Figure _StdNorm_.

figcap: The standard Normal distribution is a Normal distribution with $\mu=0$ and $\sigma=1$, $N(0,1)$ for short.

```{r StdNorm, fig.cap = "The standard Normal distribution is a Normal distribution with $\\mu=0$ and $\\sigma=1$, $N(0,1)$ for short.", echo = FALSE, fig.width = large_width}

q0 + #geom_ribbon(aes(x=value,ymax=dist),ymin=0, fill = NA) +
    scale_x_continuous("x",
                       breaks = seq(-xRangeS, xRangeS, 1),
                       labels = seq(-xRangeS, xRangeS, 1),
                       expand = c(0,0)) +
# ggtitle(  expression(paste0("The standard Normal distribution: (", mu, "= 0, ", sigma, "=1)"))) +
geom_vline(xintercept = mu, colour = MidRed) 

```

Why is this of any use to us? It enables us to immediately judge how unusual a value is with respect to a given model. Recall that 95\% of values of any Normal curve fall in the interval $[\mu-1.96\sigma, \mu+1.96\sigma]$. In the context of Z--scores where $\mu=0$ and $\sigma=1$ this interval becomes $[-1.96, 1.96]$ and we immediately know that any Z--score outside of this interval is quite an unusal value! In general, the larger the absolute value of a Z--score (the further away it is from 0), the more extreme the underlying $X$--value is with repect to the model that is used for the normalisation.^[Also note that Z--scores don't require a unit! When you calculate them, the unit appears in both the numerator and the denominator and cancels out.]

```{exercise, name="Calculating Z--scores"}
Returning to the mystery 215 cm tall Martian from the previous exercise. Calcuate it's Z--scores with respect to our two models from site I and site II.  Which site do you think the Martian came from, and why?
```

Note that you could even plot both Z--scores (e.g. as vertical lines) on the standard Normal curve. By definition they are on the same scale!
