STAT 35000 Simulation and Statistical Inference Project

STAT35000 Project 2Simulation and Statistical Inference
Abstract
There’s an ongoing debate in the academic community about whether
Calculus is a necessary pre-requisite for Statistics. But in age of ubiquitous
computing resources (not to mention open source programming languages
like R), there’s a fair argument to be made that all you really need is
simulation, which is a numerical technique for conducting experiments on
the computer involving random sampling from probability distributions.
• In this project, you will perform simulations to learn statistical inference
procedures. Then create a short write-up describing your findings.
• Most of what you will do in R were explained in R lectures. You can also
search the R Help pages online.
• All calculations and plots should be done with R and the relevant codes
and output (including numerical output and graphs) pasted into your MS
Word report. Use R output to support your answers.
• In the write-up you submit, you should use clear and complete sentences,
your numerical answers should have units attached, and your tables/graphs
should be clearly labeled.
John and Jack found a coin on the sidewalk. They argued about whether
the coin is fair. John claimed a 40% chance for the coin to land on heads based
on his careful observation of the coin. Jack believed the chance would be lower.
To support his claim Jack tossed the coin 50 times and the observed sample is
given below with 1 representing heads and 0 tails:
## [1] 0 0 0 1 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## [36] 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0
Let p∗ be the true probability for the coin to land on heads. Note that p∗ is
a characteristic of the coin, and we want to infer about this unknown parameter
(i.e. Is it lower than 40%?). Let X be the random variable reflecting a single
coin tossing result: it equals 1 if the coin lands on heads, and is 0 otherwise.
1. (10pts) What is the distribution of the random variable X?
1
2
2. (5pts) Calculate the sample proportion p̂∗ (name it phatstar in R) which
is the proportion of “1”s in the given sample above. (Note that this hat
version is your “best” guess for the unknown p∗ based on the sample, hence
is a statistic. These two notations, p∗ and p̂∗ , are indeed different.)
3. How can we help Jack use the data given above to show evidence for/against
John’s guess? We use the proof by contradiction idea. Assume John’s
claim is correct, then we have in the computer a virtual coin with the
claimed 40% chance for heads to show up. And you can simulate to see
what happens based on John’s claim, then compare it to the data given
in above sample to check whether what happens based on John’s claim
is consistent with the observed data. Since the observed sample already
happened, it must be true. And if John’s claim leads to things inconsistent with the observed data, then it must be John’s claim that is wrong.
Here is how to implement:
(i) (10pts) Assume John’s claim of 40% chance is true, use simulations
to flip this virtual coin 50 times. Repeat this process N = 1000
times. The result virtual is a 50 × 1000 matrix where each column
is a repetition.
virtual = replicate(1000, rbinom(50, 1, 0.4))
Why this: If John’s claim of p∗ = 0.4 is correct, then these 1000
samples should be comparable to the above sample Jack obtained as
based on John’s claim they all should be from the same distribution
Ber(0.4). Next we will use one way to check whether the simulated
1000 samples and the observed single sample are indeed comparable.
(ii) (15pts) Calculate the sample proportion for each of these N = 1000
samples, denoted by {p̂k50 , k = 1, 2 · · · , N = 1000}. Now John’s claim
leads us to N = 1000 simulated sample proportions.
phat50s = colSums(virtual)/50
Why this: If John’s claim of p∗ = 0.4 is correct, then the sample proportion of Jack’s observed data, p̂∗ , should be comparable
to most of the {p̂k50 , k = 1, 2 · · · , N = 1000}. (Do not be confused
about the notations: p∗ is the true unknown heads probability of
the coin, p̂∗ is the sample proportion for Jack’s particular sample.
{p̂k50 , k = 1, 2 · · · , N = 1000} are the sample proportions from simulations based on John’s claim.)
(iii) (15pts) Plot the histogram of {p̂k50 , k = 1, 2 · · · , N = 1000} and mark
the observed p̂∗ from Jack’s data in the histogram:
hist(phat50s, prob = TRUE)
points(phatstar, 0, col = “red”, pch = 20)
3
REFERENCES
Why this: Visually check whether observed proportion in Jack’s
sample is comparable to most of the 1000 simulated sample proportions based on John’s claim. If it is in the extreme region of the
histogram, you probably would conclude it is not comparable so that
Jack could use this graph to argue against John’s claim.
(iv) (15pts) To be precise, calculate the percentage of those 1000 p̂k50 ’s
that is smaller than (in the direction of Jack’s hypothesis) the observed p̂∗ in Jack’s data:
mean(phat50s < phatstar) Why this: This percentage is a measure of how extreme Jack’s observed p̂∗ is among the 1000 simulated p̂k50 ’s. This percentage is actually related to an important concept in statistics, p-value! (Do not get confused about the p in the name of p-value and other p’s above. ) (v) (10pts) John will accept Jack’s objection if the p-value is smaller than 5%, otherwise Jack has to accept John’s guess. Based on this rule, what is your conclusion for the coin’s heads proabbility? 4. In the calculation of the above percentage, we actually used simulations +···+X50 to approximate the probability of p̂50 = X1 +X250 < p̂∗ given that {Xi , i = 1, 2, · · · , 50} are independent and identically distributed as X ∼ Ber(0.4). We could use the following theoretical ways to re-calculate it. (Note that p̂50 is the random variable notation for sample proportion of 50 tosses (p̂150 , · · · , p̂1000 obtained earlier are 1000 realizations of this random 50 variable). Again p̂∗ is a real value, it is the sample proportion of Jack’s particular sample.) (a) (10pts) Use Pr(p̂50 < p̂∗ ) = Pr(X1 + X2 + · · · + X50 < 50p̂∗ ) to calculate the probability. Hint: X1 + X2 + · · · + X50 follows Binomial distribution. +···+X50 to approx(b) (10pts) Use central limit theorem for p̂50 = X1 +X250 imate the probability. References