# Testing and Estimation

### A simple example

Consider the computer-generated picture below which is supposed to mimic a photograph taken by a low resolution telescope. Is it reasonable to say that there are two distinct stars in the picture?

The next image is obtained by reducing the resolution of the telescope. Each star is now more blurred. Indeed, their distance in the picture is almost overwhelmed in the blur, and they appear to be an elongated blur caused by a single star.

The last picture is is of the same two stars taken through a very low resolution telescope. Based on this image alone there is just no reason to believe that we are seeing two stars.

What we did just now is an informal version of a statistical test of hypothesis. Each time we tried to answer the question ``Are there two stars or one?" based on a blurry picture. The two possibilities
``There is just one star'' as opposed to ``There are two stars''
are called two hypotheses. The simpler of the two gets the name null hypothesis while the other is called the alternative hypothesis. Here we shall take the first one as the null, and the second one as the alternative hypothesis.

The blurry picture we use to make a decision is called the statistical data (which always contains random errors in the form of blur/noise). Notice how our final verdict changes with the amount of noise present in the data.

Let us follow the thought process that led us to a conclusion based on the first picture. We first mentally made a note of the amount of blur. Next we imagined the centers of the bright blobs. If there are two stars then the are most likely to be here. Now we compare the distance between these centers and the amount of blur present. If the distance seems too small compared to the blur then we pass off the entire bundle as a single star.

This is precisely the idea behind most statistical tests. We shall see this for the case of the two sample t-test.

### Do Hyades stars differ in color from the rest?

Recall in the Hipparcos data set we had 92 Hyades stars and 2586 non-Hyades stars. We want to test the null hypothesis
``The Hyades stars have the same color as the non-Hyades stars''
versus the alternative hypothesis
``They have different colors.''
First let us get hold of the data for both the groups. For this we shall use the little script we had saved earlier. This script is going to load the entire Hipparcos data set, and extract the Hyades stars. So you must make sure that both HIP.dat and hyad.r are saved on your desktop. Right-click on the links in the last line, choose "Save link as..." or "Save target as..." (and be careful that the files do not get saved a .txt files). Then you must navigate R to the desktop. The simplest way to achieve this is to use the File > Change dir... menu item. This will open a pop up like this:

 Changing directory in R (Windows version)

Click on Desktop, and then OK.
```source("hyad.r")
color = B.V
H = color[HyadFilter]
nH = color[!HyadFilter & !is.na(color)]
m = length(H)
n = length(nH)
```
In the definition of `nH` above, we needed to exclude the `NA` values. `H` is a list of m numbers and ```nH ``` is a list of n numbers.

First we shall make an estimate of the ``blur'' present in the data. For this we shall compute the pooled estimate of standard deviation.
```blur.H = var(H)
blur.nH = var(nH)
blur.pool = ((m-1)*var(H) + (n-1)*var(nH))/(m+n-2)
```
Next we shall find the difference of the two means:
```meanDiff = mean(H)-mean(nH)
```
Finally we have to compare the difference with the blur. One way is to form their ratio.
```(meanDiff/sqrt(blur.pool))/sqrt(1/m + 1/n)
```
This last factor (which is a constant) is there only for technical reasons (you may think of it as a special constant to make the ``units match'').

The important question now is ``Is this ratio small or large?'' For the image example we provided a subjective answer. But in statistics we have an objective way to proceed. Before we see that let us quickly learn a one-line shortcut to compute the above ratio using the `t.test` function.
```t.test(H,nH,var.eq=T) #we shall explain the "var.eq" soon
```
Do you see the ratio in the output? Also this output tells us whether the ratio is to be considered small or large. It does so in a somewhat diplomatic way using a number called the p-value. Here the p-value is near 0, meaning
if the colors were really the same then then chance of observing a ratio this large (or larger) is almost 0.
Typically, if the p-value is smaller than 0.05 then we reject the null hypothesis. So we conclude that the mean of the color of the Hyades stars is indeed different from that of the rest.
 A rule of thumb: For any statistical test (not just t-test) accept the null hypothesis if and only if the p-value is above 0.05. Such a test fails to recognize a true null hypothesis at most 5% of the time.
The `var.eq=T` option means we are assuming that the colors of the Hyades and non-Hyades stars have more or less the same variance. If we do not want to make this assumption, we should simply write
```t.test(H,nH)
```
 Remember: `t.test` is for comparing means.

### Chi-squared tests for categorical data

Suppose that you are to summarize the result of a public examination. It is not a reasonable to report the grades obtained by each and every student in a summary report. Instead, we break the range of grades into categories like A,B,C etc and then report the numbers of students in each category. This gives an overall idea about the distribution of grades.

The `cut` function in R does precisely this.
```bvcat = cut(color, breaks=c(-Inf,0.5,0.75,1,Inf))
```
Here we have broken the range of values of the `B.V` variable into 4 categories:
(-Inf, 0.5], (0.5, 0.75], (0.75, 1] and (1,Inf).
The result (stored in `bvcat`) is a vector that records the category in which each star falls.
```bvcat
table(bvcat)
plot(bvcat)
```
It is possible to tabulate this information for Hyades and non-Hyades stars in the same table.
```table(bvcat,HyadFilter)
```
To perform a chi-squared test of the null hypothesis that the true population proportions falling in the four categories are the same for both the Hyades and non-Hyades stars, use the `chisq.test` function:
```chisq.test(bvcat,HyadFilter)
```
Since we already know these two groups differ with respect to the `B.V` variable, the result of this test is not too surprising. But it does give a qualitatively different way to compare these two distributions than simply comparing their means.

The test above is usually called a chi-squared test of homogeneity. If we observe only one sample, but we wish to test whether the categories occur in some pre-specified proportions, a similar test (and the same R function) may be applied. In this case, the test is usually called the chi-squared test of goodness-of-fit. We shall see an example of this next.

Consider once again the Hipparcos data. We want to know if the stars in the Hipparcos survey come equally from all corners of the sky. In fact, we shall focus our attention only on the `RA` values. First we shall break the range of RA into 20 equal intervals (each of width 18 degrees), and find how many stars fall in each bin.
```count = table(cut(RA,breaks=seq(0,360,len=20)))
chisq.test(count)
```

### Kolmogorov-Smirnov Test

There is yet another way (a better way in many situations) to perform the same test. This is called the Kolmogorov-Smirnov test.
```ks.test(RA,"punif",0,360)
```
Here `punif` is the name of the distribution with which we are comparing the data. `punif` denoted the uniform distribution, the range being from 0 to 360. Thus, here we are testing if `RA` is taking all values from 0 to 360 with equal likelihood, or are some values being taken more or less frequently.

The Kolmogorov-Smirnov test has the advantage that we do not need to group the data into categories as for the chi-squared test.
 Remember: `chisq.test` and `ks.test` are for comparing distributions.

### Estimation

Finding (or, rather, guessing about) the unknown based on approximate information (data) is the aim of statistics. Testing hypotheses is one aspect of it where we seek to answer yes-no questions about the unknown. The problem of estimation is about guessing the values of unknown quantities.

There are many methods of estimation, but most start off with a statistical model of the data. This is a statement of how the observed data set is (probabilistically) linked with the unknown quantity of interest.

For example, if I am asked to estimate p based on the data
Head, Head, Head, Tail, Tail, Head, Head, Tail, Head, Tail, Tail, Head
then I cannot make head-or-tail of the question. I need to link this data set with p through a statement like
A coin with probability p was tossed 12 times and the data set was the result.
This is a statistical model for the data. Now the problem of estimating p from the data looks like a meaningful one.
 Statistical models provide the link between the observed data and the unknown reality. They are indispensable in any statistical analysis. Misspecification or over-simplification of the statistical model is the most frequent cause behind misuse of statistics.

#### Estimation using R

You might be thinking that R has some in-built tool that can solve all estimation problems. Well, there isn't. In fact, due to the tremendous diversity among statistical models and estimation methods no statistical software can have tools to cope with all estimation problems. R tackles estimation in three major ways.
1. Many books/articles give formulas to estimate various quantities. You may use R as a calculator to implement them. With enough theoretical background you may be able to come up with your own formulas that R will happily compute for you.
2. Sometimes estimation problems lead to complicated equations. R can solve such equations for you numerically.
3. For some frequently used statistical methods (like regression or time series analysis) R has the estimation methods built into it.
To see estimation in action let us load the Hipparcos data set.
```hip = read.table("HIP.dat", header = TRUE, fill = TRUE)
attach(hip)
Vmag
```
If we assume the statistical model that the `Vmag` variable has a normal distribution with unknown mean μ and unknown variance σ2 then it is known from the literature that a good estimator of μ is and a 95% confidence interval is
where
This means that the true value of μ will lie in this interval with 95% chance. You may think of μ as a peg on the wall and the confidence interval as a hoop thrown at it. Then the hoop will miss the peg in only about 5% of the cases.

 Exercise: Find the estimate and confidence interval for μ based on the observed values of `Vmag` using R as a calculator.

Next we shall see a less trivial example. The data set comes from NASA's Swift satellite. The statistical problem at hand is modeling the X-ray afterglow of gamma ray bursts. First, read in the dataset GRB.dat (right-click, save on your desktop, without any .txt extension).
```dat = read.table("GRB.dat",head=T)
flux = dat[,2]
```
Suppose that it is known that the `flux` variable has an Exponential distribution. This means that its density function is of the form
Here λ is a parameter, which must be positive. To get a feel of the density function let us plot it for different values of λ
```x = seq(0,200,.1)
y = dexp(x,1)
plot(x,y,ty="l")
```
Now let us look at the histogram of the observed `flux` values.
```hist(flux)
```
The problem is estimating λ based on the data may be considered as finding a value of λ such that the the density is as close as possible to the histogram.

We shall first try to achieve this interactively. For this you need to download interact.r on your desktop first.
```source("interact.r")
```
This should open a tiny window as shown below with a slider and a button in it.

 Screenshot of the tiny window

Move the slider to see how the density curve moves over the histogram. Choose a position of the slider for which the density curve appears to be a good approximation.

 Exercise: It is known from the theory of Exponential distribution that the Maximum Likelihood Estimate (MLE) of λ is the reciprocal of the sample mean Compute this estimate using R and store it in a variable called `lambdaHat`. Draw the density curve on top of the histogram using the following commands. ```y = dexp(x,lambdaHat) lines(x,y,col="blue") ```

### Some nonparametric tests

The remaining part of this tutorial deals with a class of inference procedures called nonparametric inference. These may be skipped without loss of continuity. Also the theory part for this will be covered in a later theory class. So you may like to save the rest of this lab until that time. I have tried to make the lab largely self-explanatory though, with a little theoretical discussion to explain terms and concepts.

#### Distributions

In statistics we often come across questions like
``Is a given data from such-n-such distribution?''
Before we can answer this question we need to know what is meant by a distribution. We shall talk about this first.

In statistics, we work with random variables. Chance controls their values.

The distribution of a random variable is the rule by which chance governs the values. If you know the distribution then you know the chance of this random variable taking value in any given range. This rule may be expressed in various ways. One popular way is via the Cumulative Distribution Function (CDF), that we illustrate now. If a random variable (X) has distribution with CDF F(x) then for any given number a the chance that X will be ≤ a is F(a).

Example: Suppose that I choose a random star from the Hipparcos data set. If I tell you that its `Vmag` is a random variable with the following CDF then what is the chance that it takes values ≤ 10? Also find out the value x such that the chance of `Vmag` being ≤ x is 0.4.

For the first question the required probability is F(10) which, according to the graph, is slightly above 0.8.

In the second part we need F(x) = 0.4. From the graph it seems to be about 7.5.

Just to make sure that these are indeed meaningful, let us load the Hipparcos data set and check:
```hip = read.table("HIP.dat", header = TRUE, fill = TRUE)
attach(hip)
n = length(Vmag) #total number of cases
count = sum(Vmag<=10) #how many <= 10
count/n #should be slightly above 0.8
```
This checks the answer of the first problem. For the second
```count = sum(Vmag<=7.5) #how many <= 7.5
count/n #should be around 0.4
```

So you see how powerful a CDF is: in a sense it stores as much information as the entire `Vmag` data. It is a common practice in statistics to regard the distribution as the ultimate truth behind the data. We want to infer about the underlying distribution based on the data. When we look at a data set we actually try to look at the underlying distribution through the data set!

R has many standard distributions already built into it. This basically means that R has functions to compute their CDFs. These functions all start with the letter p.

 Example: A random variable has standard Gaussian distribution. We know that R computes its CDF using the function `pnorm`. How to find the probability that the random variable takes values ≤ 1? The answer may be found as follows: ```pnorm(1) ```

For every p-function there is a q-function that is basically its inverse.

 Example: A random variable has standard Gaussian distribution. Find x such that the random variable is ≤ x with chance 0.3. Now we shall use the function `qnorm`: ```qnorm(0.3) ```

OK, now that we are through our little theory session, we are ready for the nonparametric lab.

### One sample nonparametric tests

#### Sign-test

In a sense this is the simplest possible of all tests. Here we shall consider the data set LMC.dat that stores the measured distances to the Large Magellanic Cloud. (As always, you'll need to save the file on your desktop.)
```LMC = read.table("LMC.dat",head=T)
data = LMC[,2] #These are the measurements
data
```
We want to test if the measurements exceed 18.41 on an average. Now, this does not mean whether the average of the data exceeds 18.41, which is trivial to find out. The question here is actually about the underlying distribution. We want to know if the median of the underlying distribution exceeds 18.41, which is a less trivial question, since we are not given that distribution.

We shall use the sign test to see if the median is 18.41 or larger. First it finds how many of the observations are above this value:
```abv = sum(data>18.41)
abv
```
Clearly, if this number is large, then we should think that the median (of the underlying distribution) exceeds 18.41. The question is how large is ``large enough''? For this we consult the binomial distribution to get the p-value. The rationale behind this should come from the theory class.
```n = nrow(LMC)
pValue = 1-pbinom(abv-1,n,0.5)
pValue
```
We shall learn more about p-values in a later tutorial. For now, we shall use the following rule of thumb:
If this p-value is below 0.05, say, we shall conclude that the median is indeed larger than 18.41.

#### Wilcoxon's Signed Rank Test

As you should already know from the theoretical class, there is a test called Wilcoxon's Signed Rank test that is better (if a little more complicated) than the sign test. R provides the function `wilcox.test` for this purpose. We shall work with the Hipparcos data set this time, which we have already loaded. We want to see if the median of the distribution of the `pmRA` variable is 0 or not.
```wilcox.test(pmRA, mu=0)
```
Since the p-value is pretty large (above 0.05, say) we shal conclude that the median is indeed equal to 0. R itself has come to the same conclusion.

Incidentally, there is a little caveat. For Wilcoxon's Rank Sum Test to be applicable we need the underlying distribution to be symmetric around the median. To get an idea about this we should draw the histogram.
```hist(pmRA)
```
Well, it looks pretty symmetric (around a center line). But would you apply Wilcoxon's Rank Sum Test to the variable `DE`?
```hist(DE)
```
No, this is not at all symmetric!

### Two-sample nonparametric tests

Here we shall be taking about the shape of distributions. Let us first make sure of this concept. Suppose that I make a list of all the people in a capitalist country and make a histogram of their incomes. I should get a histogram like this

 Income histogram of a capitalist country

But if the same thing is done for a socialist country we shall see something like this (where almost everybody is in the middle income group).

 Income histogram of a socialist country

Clearly, the shapes of the histograms differ for the two populations. We must be careful about the term ``shape'' here. For example, the following two histograms are for two capitalist countries, one rich the other poor.

 Income histograms two capitalist countries

Here the histograms have the same shape but differ in location. You can get one from the other by applying a shift in the location.

This is a very common phenomenon in statistics. If we have two populations with comparable structure they often have similar shapes, and differ in just a location shift. In such a situation we are interested in knowing about the amount of shift. if the amount of shift is zero, then the two populations behave identically.

Wilcoxon's rank sum test is one method to learn about the amount of shift based on samples from the two populations.

We shall start with the example done in the theory class, where we had two samples
```x = c(37,49,55,57)
y = c(23,31,46)
m = length(x)
n = length(y)
```
We are told that both the samples come from populations with the same shape, though one may be a shifted version of the other.

Our aim is to check if indeed there is any shift or not.

We shall consider the pooled data (i.e., all the 7 numbers taken together) and rank them
```pool = c(x,y)
pool
r = rank(pool)
r
```
The first m=4 of these ranks are for the first sample. Let us sum them:
```H = sum(r[1:m])
H
```
If the the two distributions are really identical (i.e., if there is no shift) then we should have (according to the theory class) H close to the value
```m*(m+n+1)/2   #Remember: * means multiplication
```
Can the H that we computed from our data be considered ``close enough" to this value? We shall let us R determine that for us.
```wilcox.test(x,y)
```
So R clearly tells us that there is a shift in location. The output also mentions a W and a p-value. tutorial. But W is essentially what we had called H. Well, H was 21, while W is 11. This is because R has the habit of subtracting
m(m+1)/2
from H and calling it W. Indeed for our data set
m(m+1)/2 = 10,
which perfectly accounts for the difference between our H and R's W.

You may recall from the theory class that H is called Wilcoxon's rank sum statistic, while W is the Mann-Whitney statistic (denoted by U in class).

Now that R has told us that there is a shift, we should demand to estimate the amount of shift.
```wilcox.test(x,y,conf.int=T)
```
The output gives two different forms of estimates. It gives a single value (a point estimate) which is 16. Before it gives a 95% confidence interval: -9 to 34. This basically means that
we can say with 95% confidence that the true value of the shift is between -9 and 34.
We shall talk more about confidence intervals later. For now let us see how R got the poit estimate 16. It used the so-called Hodges-Lehman formula, that we have seen in the theoretical class. To see how this method works we shall first form the following "difference table"
```
| 23 31 46
------------
37| 14  6 -9
49| 26 18  3
55| 32 24  9
57| 34 26 11

```
Here the rows are headed by the first sample, and columns by the second sample. Each entry is obtained by subtracting the column heading from the row heading (e.g., -9 = 37-46). This table, by the way, can be created very easily with the function `outer`.
```outer(x,y,"-")
```
 The `outer` function is a very useful (and somewhat tricky!) tool in R. It takes two vectors x, and y, say, and some function f(x,y). Then it computes a matrix where the (i,j)-th entry is f(x[i],y[j]).
Now we take the median of all the entries in the table:
```median(outer(x,y,"-"))
```
This gives us the Hodges-Lehman estimate of the location shift.