---
title: "Computing Distribution Quantiles, and fun with Chi-Squared"
output: html_document
---

```{r setup, include=TRUE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Let's say that we have some data distributed with $N(0, 1.2^2)$. Let's generate the data:
```{r generateData}
N = 140
sigma <- 1.2
mu <- 0
dat <- rnorm(n=N, mean=mu, sd= sigma)
```

Suppose the Null Hypothesis is that $\sigma = 1.0$. How can we test it? Let's construct an interval in which $\sum_i \left( \frac{X_i-\bar{X}}{1.0}\right)^2$ will fall 95\% of the time assuming the Null Hypothesis is true. That statistic is distributed according to $\chi^2(N-1)$ (i.e., $\chi^2(139)$) if the Null Hypothesis is true.

Compute the statistic:

```{r}
sum_res <- sum(  ((dat-mean(dat))/1.0)^2 )
sum_res
```

Now, compute the interval:

```{r}
c(qchisq(.025, 139), qchisq(0.975, 139))
```

We can reject the Null Hypothesis.

Note that, *assuming we decide ahead of time which interval to compute*, we could go with a different interval as well. For example, here's another legitimate interval.

```{r}
c(qchisq(.01, 139), qchisq(0.96, 139))
```

Why can't we pick an choose which interval to use after the fact?



How to compute stuff like `qchisq`, `qt`, or `qnorm` anyway? Let's say we need `qnorm(.95)`

Here's how to compute `qt`. First, generate lots of samples from the normal distribution:

```{r}
N = 10000
samples <- rnorm(n=N)
```

Now, we need to find a number $t$ such that samples from the normal are lower than $t$ 95% of the time. We can do it by sorting our samples, and then take sample number $round(.95N)$ in the sorted vector.

```{r}
sorted_samples <- sort(samples)
sorted_samples[round(.95*N)]
```

Did we get it right?

```{r}
qnorm(.95)
```

Looks like it. 

Now, let's try to do the same thing with $\chi^2(10)$. I'll use a loop even though it's a little inefficient

```{r}
chisq10 <- c()
for (i in 1:N){
  chisq10 <- c(chisq10, sum(rnorm(10)^2))
}

sorted_samples <- sort(chisq10)
sorted_samples[round(.95*N)]

```

Did we get it right?

```{r}
qchisq(.95, df=10)
```

Looks like it!
