Topics

References

  • Wasserman (2004), Chpter 6
  • Motulsky (2014), Chapters 4, 12, 15-17

Statistical Inference

Given the sample \(X_1,...,X_n\), infer the distribution \(F\) that generated it.

Statistical model: a set of probability distributions

Standard Error

Estimate \(\hat{\theta}\) of a paramter \(\theta\) from data \(X_1,...,X_n\) is stochastic.

  • Bias: \(\mbox{bias}(\hat{\theta}) = E(\hat{\theta}) - \theta\)

  • Standard error: \(\mbox{se}(\hat{\theta}) = \sqrt{V(\hat{\theta})}\)

  • Estimated standard error: \(\hat{\mbox{se}}\)

\(X_1,...,X_n \sim \mbox{Bernoulli}(p)\)

  • Estimate \(\hat{p}=\frac{1}{n}\sum_i X_i\) is unbiased.

  • Standard error: \(\mbox{se}(\hat{p}) = \sqrt{V(\hat{p})} = \sqrt{\frac{p(1-p)}{n}}\)

  • Estimated standard error: \(\hat{\mbox{se}}(\hat{p}) = \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\)

n = 20   # sample size
m = 100  # repeated sampling
p = 0.1  # true p
phat = rbinom(n=m, size=n, prob=p)/n
plot(phat)
lines(c(1,m), c(p,p), type="l", col="red")

mean(phat)-p  # bias
[1] -0.0025
sqrt(var(phat))  # measured standard error
[1] 0.06604827
sqrt(p*(1-p)/n)  # theoretical standard error
[1] 0.06708204

Confidence Interval of Proportion

From sample {0,1,0,0,1,0,…}, infer \(p = \mbox{Prob}(X=1)\)

  • sum \(S = \sum_i X_i\)

  • Estimate \(\hat{p} = \frac{S}{n}\)

What is 95% confidence interval?

  • Range where true \(p\) resides 95% probability??

  • If we compute the intervals from multiple samples, 95% of them should include the true \(p\).

Modified Wald Method

Find a range \([-z, z]\) that covers \(1-\alpha\) of standard Gaussian.

  • e.g., \(z=1.96\) for \(alpha=0.05\)

  • calculate the center \(p’ = \frac{S+0.5z^2}{n+z^2}\)

  • calculate the width \(W = z\sqrt{\frac{p’(1–p’)}{n+z^2}}\)

  • Confidence interval is \([p’–W, p’+W]\)

n = 10
m = 100
p = 0.2
z = 1.96  # for alpha=0.05
plot(c(1,m), c(p,p), type="l", col="red", ylim=c(0,1))
CI = matrix(0, m, 2)
for (j in 1:m) {
    s = rbinom(n=1, size=n, prob=p)
    phat = s/n
    pp = (s+0.5*z^2)/(n+z^2)
    w = z*sqrt(pp*(1-pp)/(n+z^2))
    CI[j,1] = pp - w
    CI[j,2] = pp + w
    lines(c(j,j), CI[j,], col="blue")
}

Confidence Interval of the Mean

  • sample mean \(\hat{\mu} = \frac{1}{n}\sum_{i=1}^n X_i\) varies with samples

  • t distribution presents how a sample mean deviates from the true mean.

\[f(t,\nu) \propto (1+\frac{t^2}{\nu})^{–\frac{\nu+1}{2}}\]

  • \(\nu\): degree of freedom: Gaussian with \(\nu \rightarrow \infty\)
x = seq(-5, 5, 0.1)
f = dnorm(x)  # Gaussian as reference
plot(x, f, type="l")
nu = c(1, 2, 5, 10)
for (i in 1:length(nu)){
    f = dt(x, nu[i])
    lines(x, f, col=i+1)
}

  • Find \([-t^*,t^*]\) that covers \(1-\alpha\) of t distribution with \(\nu=n-1\).

  • Width: \(W = t^* \mbox{se}(\hat{\mu}) = t^* \frac{\hat{\sigma}}{\sqrt{n}}\)

  • Confidence interval: \([\hat{\mu}–W, \hat{\mu}+W]\)

Exercise

1. Confidence Interval of the Mean

Take m samples of n data from a Gaussian distribution and estimate their means.

  1. See how they are distributed for different n.
  1. Compare them with t distribution with different \(\nu\)
  1. Compute a 95% confidence intervals of the means and check how many of them miss the true mean.
LS0tCnRpdGxlOiAiMy4gU3RhdGlzdGljYWwgSW5mZXJlbmNlIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgojIFRvcGljcwoKKiBTdGF0aXN0aWNhbCBpbmZlcmVuY2UKKiBTdGFuZGFyZCBlcnJvcgoqIENvbmZpZGVuY2UgaW50ZXJ2YWwgb2YgcHJvcG9ydGlvbgoqIENvbmZpZGVuY2UgaW50ZXJ2YWwgb2YgdGhlIG1lYW4KCiMjIFJlZmVyZW5jZXMKCiogV2Fzc2VybWFuICgyMDA0KSwgQ2hwdGVyIDYKKiBNb3R1bHNreSAoMjAxNCksIENoYXB0ZXJzIDQsIDEyLCAxNS0xNwoKIyBTdGF0aXN0aWNhbCBJbmZlcmVuY2UKCkdpdmVuIHRoZSBzYW1wbGUgJFhfMSwuLi4sWF9uJCwgaW5mZXIgdGhlIGRpc3RyaWJ1dGlvbiAkRiQgdGhhdCBnZW5lcmF0ZWQgaXQuCgoqIGNoYXJhY3Rlcml6YXRpb24gb2YgZ2l2ZW4gZGF0YQoqIGZpbmRpbmcgdGhlIHN0cnVjdHVyZSBiZWhpbmQgZGF0YQoqIHByZWRpY3Rpb24gb2YgbmV3IGRhdGEKClN0YXRpc3RpY2FsIG1vZGVsOiBhIHNldCBvZiBwcm9iYWJpbGl0eSBkaXN0cmlidXRpb25zCgoqIE5vbnBhcmFtZXRyaWMgbW9kZWw6IG5vIHByZWZpeGVkIGZvcm0gYW5kIHBhcmFtZXRlcml6YXRpb24KICAgICsgZS5nLiwgaGlzdG9ncmFtCgoqIFBhcmFtZXRyaWMgbW9kZWw6IGRpc3RyaWJ1dGlvbnMgd2l0aCBhIHBhcnRpY3VsYXIgZm9ybSAkZih4O1x0aGV0YSkkIGFuZCBwYXJhbWV0ZXJzICQgXHRoZXRhXGluXFRoZXRhJC4KCiogUGFyYW1ldGVyIGVzdGltYXRpb246IFdlIG9mdGVuIGFzc3VtZSBhIHBhcmFtZXRlcmljIG1vZGVsICRGKFg7XHRoZXRhKSAkIGFuZCBpbmZlciBpdHMgcGFyYW10ZXIgJFx0aGV0YSQuCgogICAgKiBleGFtcGxlczogZXN0aW1hdGUgdGhlIG1lYW4gJHAkIG9mIGEgQmVybm91bGxpIGRpc3RyaWJ1dGlvbgogICAgKiBlc3RpbWF0ZSB0aGUgbWVhbiAkXG11JCBhbmQgdGhlIHZhcmlhbmNlICRzaWdtYV4yJCBvZiBhIE5vcm1hbCBkaXN0cmlidXRpb24KICAgIAojIyBTdGFuZGFyZCBFcnJvcgoKRXN0aW1hdGUgJFxoYXR7XHRoZXRhfSQgb2YgYSBwYXJhbXRlciAkXHRoZXRhJCBmcm9tIGRhdGEgJFhfMSwuLi4sWF9uJCBpcyBzdG9jaGFzdGljLgoKKiBCaWFzOiAkXG1ib3h7Ymlhc30oXGhhdHtcdGhldGF9KSA9IEUoXGhhdHtcdGhldGF9KSAtIFx0aGV0YSQKCiogU3RhbmRhcmQgZXJyb3I6ICRcbWJveHtzZX0oXGhhdHtcdGhldGF9KSA9IFxzcXJ0e1YoXGhhdHtcdGhldGF9KX0kCgoqIEVzdGltYXRlZCBzdGFuZGFyZCBlcnJvcjogJFxoYXR7XG1ib3h7c2V9fSQKCiMjIyAkWF8xLC4uLixYX24gXHNpbSBcbWJveHtCZXJub3VsbGl9KHApJAoKKiBFc3RpbWF0ZSAkXGhhdHtwfT1cZnJhY3sxfXtufVxzdW1faSBYX2kkIGlzIHVuYmlhc2VkLgoKKiBTdGFuZGFyZCBlcnJvcjogJFxtYm94e3NlfShcaGF0e3B9KSA9IFxzcXJ0e1YoXGhhdHtwfSl9ID0gXHNxcnR7XGZyYWN7cCgxLXApfXtufX0kCgoqIEVzdGltYXRlZCBzdGFuZGFyZCBlcnJvcjogJFxoYXR7XG1ib3h7c2V9fShcaGF0e3B9KSA9IFxzcXJ0e1xmcmFje1xoYXR7cH0oMS1caGF0e3B9KX17bn19JAoKYGBge3J9Cm4gPSAyMCAgICMgc2FtcGxlIHNpemUKbSA9IDEwMCAgIyByZXBlYXRlZCBzYW1wbGluZwpwID0gMC4xICAjIHRydWUgcApwaGF0ID0gcmJpbm9tKG49bSwgc2l6ZT1uLCBwcm9iPXApL24KcGxvdChwaGF0KQpsaW5lcyhjKDEsbSksIGMocCxwKSwgdHlwZT0ibCIsIGNvbD0icmVkIikKbWVhbihwaGF0KS1wICAjIGJpYXMKc3FydCh2YXIocGhhdCkpICAjIG1lYXN1cmVkIHN0YW5kYXJkIGVycm9yCnNxcnQocCooMS1wKS9uKSAgIyB0aGVvcmV0aWNhbCBzdGFuZGFyZCBlcnJvcgpgYGAKCiMjIENvbmZpZGVuY2UgSW50ZXJ2YWwgb2YgUHJvcG9ydGlvbgoKRnJvbSBzYW1wbGUgezAsMSwwLDAsMSwwLC4uLn0sIGluZmVyICRwID0gXG1ib3h7UHJvYn0oWD0xKSQKCiogc3VtICRTID0gXHN1bV9pIFhfaSQKCiogRXN0aW1hdGUgJFxoYXR7cH0gPSBcZnJhY3tTfXtufSQKCiMjIyBXaGF0IGlzIDk1JSBjb25maWRlbmNlIGludGVydmFsPwoKKiBSYW5nZSB3aGVyZSB0cnVlICRwJCByZXNpZGVzIDk1JSBwcm9iYWJpbGl0eT8/CgoqIElmIHdlIGNvbXB1dGUgdGhlIGludGVydmFscyBmcm9tIG11bHRpcGxlIHNhbXBsZXMsIDk1JSBvZiB0aGVtIHNob3VsZCBpbmNsdWRlIHRoZSB0cnVlICRwJC4KCiMjIyBNb2RpZmllZCBXYWxkIE1ldGhvZAoKRmluZCBhIHJhbmdlICRbLXosIHpdJCB0aGF0IGNvdmVycyAkMS1cYWxwaGEkIG9mIHN0YW5kYXJkIEdhdXNzaWFuLiAgCgoqIGUuZy4sICR6PTEuOTYkIGZvciAkYWxwaGE9MC4wNSQKCiogY2FsY3VsYXRlIHRoZSBjZW50ZXIgJHDigJkgPSBcZnJhY3tTKzAuNXpeMn17bit6XjJ9JAoKKiBjYWxjdWxhdGUgdGhlIHdpZHRoICRXID0gelxzcXJ0e1xmcmFje3DigJkoMeKAk3DigJkpfXtuK3peMn19JAoKKiBDb25maWRlbmNlIGludGVydmFsIGlzICRbcOKAmeKAk1csIHDigJkrV10kCgpgYGB7cn0KbiA9IDEwCm0gPSAxMDAKcCA9IDAuMgp6ID0gMS45NiAgIyBmb3IgYWxwaGE9MC4wNQpwbG90KGMoMSxtKSwgYyhwLHApLCB0eXBlPSJsIiwgY29sPSJyZWQiLCB5bGltPWMoMCwxKSkKQ0kgPSBtYXRyaXgoMCwgbSwgMikKZm9yIChqIGluIDE6bSkgewogICAgcyA9IHJiaW5vbShuPTEsIHNpemU9biwgcHJvYj1wKQogICAgcGhhdCA9IHMvbgogICAgcHAgPSAocyswLjUqel4yKS8obit6XjIpCiAgICB3ID0geipzcXJ0KHBwKigxLXBwKS8obit6XjIpKQogICAgQ0lbaiwxXSA9IHBwIC0gdwogICAgQ0lbaiwyXSA9IHBwICsgdwogICAgbGluZXMoYyhqLGopLCBDSVtqLF0sIGNvbD0iYmx1ZSIpCn0KYGBgCgojIyBDb25maWRlbmNlIEludGVydmFsIG9mIHRoZSBNZWFuIAoKKiBzYW1wbGUgbWVhbiAkXGhhdHtcbXV9ID0gXGZyYWN7MX17bn1cc3VtX3tpPTF9Xm4gWF9pJCB2YXJpZXMgd2l0aCBzYW1wbGVzCgoqIHQgZGlzdHJpYnV0aW9uIHByZXNlbnRzIGhvdyBhIHNhbXBsZSBtZWFuIGRldmlhdGVzIGZyb20gdGhlIHRydWUgbWVhbi4KCiQkZih0LFxudSkgXHByb3B0byAoMStcZnJhY3t0XjJ9e1xudX0pXnvigJNcZnJhY3tcbnUrMX17Mn19JCQKCiogJFxudSQ6IGRlZ3JlZSBvZiBmcmVlZG9tOiBHYXVzc2lhbiB3aXRoICRcbnUgXHJpZ2h0YXJyb3cgXGluZnR5JAoKYGBge3J9CnggPSBzZXEoLTUsIDUsIDAuMSkKZiA9IGRub3JtKHgpICAjIEdhdXNzaWFuIGFzIHJlZmVyZW5jZQpwbG90KHgsIGYsIHR5cGU9ImwiKQpudSA9IGMoMSwgMiwgNSwgMTApCmZvciAoaSBpbiAxOmxlbmd0aChudSkpewogICAgZiA9IGR0KHgsIG51W2ldKQogICAgbGluZXMoeCwgZiwgY29sPWkrMSkKfQpgYGAKCiogRmluZCAkWy10XiosdF4qXSQgdGhhdCBjb3ZlcnMgJDEtXGFscGhhJCBvZiB0IGRpc3RyaWJ1dGlvbiB3aXRoICRcbnU9bi0xJC4KCiogV2lkdGg6ICRXID0gdF4qIFxtYm94e3NlfShcaGF0e1xtdX0pID0gdF4qIFxmcmFje1xoYXR7XHNpZ21hfX17XHNxcnR7bn19JAoKKiBDb25maWRlbmNlIGludGVydmFsOiAkW1xoYXR7XG11feKAk1csIFxoYXR7XG11fStXXSQgCgojIEV4ZXJjaXNlCgojIyAxLiBDb25maWRlbmNlIEludGVydmFsIG9mIHRoZSBNZWFuCgpUYWtlIG0gc2FtcGxlcyBvZiBuIGRhdGEgZnJvbSBhIEdhdXNzaWFuIGRpc3RyaWJ1dGlvbiBhbmQgZXN0aW1hdGUgdGhlaXIgbWVhbnMuCgpgYGB7cn0KCmBgYAoKMSkgU2VlIGhvdyB0aGV5IGFyZSBkaXN0cmlidXRlZCBmb3IgZGlmZmVyZW50IG4uCgpgYGB7cn0KCmBgYAoKMikgQ29tcGFyZSB0aGVtIHdpdGggdCBkaXN0cmlidXRpb24gd2l0aCBkaWZmZXJlbnQgJFxudSQKCmBgYHtyfQoKYGBgCgozKSBDb21wdXRlIGEgOTUlIGNvbmZpZGVuY2UgaW50ZXJ2YWxzIG9mIHRoZSBtZWFucyBhbmQgY2hlY2sgaG93IG1hbnkgb2YgdGhlbSBtaXNzIHRoZSB0cnVlIG1lYW4uCgpgYGB7cn0KCmBgYAo=