some statistics notes

Some Statistics Notes PhL 2.7-10.15 These notes are an addendum to already existing noted in Ng Notes and Fourier Appendix G. A. Review of Ng notes and Throwing one Die 2 B. Theorems about Variance and Covariance 3 Theorems 1 through 6 4 Theorem 7 (Standard Deviation of the Mean Theorem) 7 C. The Binomial Distribution and the Bernoulli Trial 8 Theorem 8: (Bernoulli Trial SD of the Mean Theorem) 11 (a) I try but fail to show that Binomial Distribution approaches a Gaussian for large N. 12 (b) Various Binomial Notations and Access from Excel or Maple 13 (c) The Central Limit Theorem and my Second Failed Attempt to show Bin → Gaussian 15 (d) Application of Bernoulli/Binomial to HAS-BLED 19 __________________________ Summary Section __________________________________ Here is a quick statement of the 8 theorems: Theorem 1: If X and Y are independent random variables, with p(x) and p(y) arbitrary, then var(X+Y) = var(X) + var(Y). Theorem 2: cov(X,Y) = E(XY) - E(X)E(Y) = Σx,y xy p(x,y) - μxμy . Theorem 3: var(X+Y) = var(X) + var(Y) + 2 cov(X,Y) . Theorem 4: var(aX+bY) = a2var(X) + b2var(Y) + 2ab cov(X,Y) Theorem 5. var(ΣiaiXi) = Σij aiaj cov(Xi,Xj) Corollary 5a. var(ΣiXi) = Σij cov(Xi,Xj) just setting all ai = 1 Corollary 6a: var(ΣiXi) = Σi var(Xi) if all Xi are independent Theorem 7 (Standard Deviation of the Mean Theorem) For N repeats of the same experiment which have outcomes x (for random variable X), if the experiments are independent (they are) then σ() = σ/ where σ() is the standard deviation of the mean and σ is the standard deviation for one experiment. The pdf shape for the single experiment can be anything. Theorem 8: (Bernoulli Trial SD of the Mean Theorem) . The standard deviation of the mean where the underlying process is a simple Bernoulli trial is given by σ() = . Application: This theorem applies to "population studies" (including medical trials). Each "experiment" involves a different person (like a die throw), and the outcome might be how tall that person is, or whether or not that person had a major bleed over the period of a year. For height, you assume for 7 billion people there exists some μ = and σ. You do an experiment with N people, just measuring height. You then plot all your heights in a histogram chart (they are all between 4' and 7', with mean 6', say). From this chart you come up with μ1 and σ1 -- these are your experiment's estimates for the true values μ and σ. Then your conclusion is that the average height is μ1 = 6' and the standard deviation of mean is σ() = σ/. Obviously large N reduces this standard deviation. The true mean μ for the full population lies within a few σ() values of your μ1 experimental estimate. This conclusion is valid regardless of the nature of the underlying probability function which describes human height distribution. It could be a binomial distribution, a normal distribution, or any other kind of distribution, symmetric or asymmetric about the mean. The symbol "σ" means , and every distribution, symmetric or not about the mean, has some variance var = Σx p(x) (x-μ)2. __________________________ End of Summary Section __________________________________ A. Review of Ng notes and Throwing one Die First, from the wiki page http://en.wikipedia.org/wiki/Standard_error SEM = standard error of the mean = σ/ σ = standard deviation (they call it s) "of the population" = From my Ng notes we have, var(X) = !Syntax Error, I (x-μx)2 p(x) dx. μx = !Syntax Error, I x p(x) dx . var(X1) = (1/n) Σi=1n (x1(i)– μX1)2 μX1 = (1/n) Σi=1n x(i)1 p = 1/n σx2 ≡ varx = E(X2) - μx2 In my Fourier Transform doc I have a huge appendix G on probability theory, hurray, I forgot all about it! I have just read it. I understand how to compute things for throwing 1 or 2 or N dice. Dice Review. If you roll a fair die one time, you get this pmf, For this distribution, you can compute whatever statistical object you want. For example, μ = Σx x p(x) = Σx x (1/6) = (1/6) Σxx = (1/6) [ 1 + 2 + 3 +4 + 5 + 6] = (1/6) [21] = 7/2 = 3.5 which seems pretty reasonable, see picture above. Another object: E(X2) = Σx x2 p(x) = (1/6) Σxx2 = (1/6) [ 12 + 22 + 32 +42 + 52 + 62] = (1/6) [ 1 + 4 + 9 +16 + 25+36] = (1/6) [ 91] = 91/6 = [7*13/6] = 15.16666 var = E(X2) - μx2 = 91/6 - (7/2)2 = 15.16666 - (3.5)*(3.5) = 2.91666 Well, in Appendix G I compute this by brute force enumerating all the terms and find that var(Z) = 5.833 Notice that this is exactly twice the variance of the single die experiment! 2* 2.91666 = 5.83332 Why is this the case?? I will now answer that question with Theorem 1 below. B. Theorems about Variance and Covariance Consider rolling two such dice and examine z = x+y. We then have p(z) = Σzx,y p(x,y) // App G E(Z) = Σz z p(z) = Σz z Σzx,y p(x,y) = Σz Σzx,y [zp(x,y)] = Σx,y (x+y)p(x,y) Here I have used my Lemma (G.15a). Next, E(Z2) = Σz z2 p(z) = Σz z2 Σzx,y p(x,y) = Σz Σzx,y [z2p(x,y)] = Σx,y (x+y)2p(x,y) Therefore var(Z) = E(Z2) - [E(z)]2 = Σx,y (x+y)2p(x,y) - [ Σx,y (x+y)p(x,y)]2 Now suppose we have p(x,y) = p(x)p(y). Then var(Z) = Σx,y (x+y)2p(x)p(y) - [ Σx,y (x+y)p(x)p(y)]2 = ΣxΣy (x+y)2p(x)p(y) - [ ΣxΣy (x+y)p(x)p(y)]2 = Σxp(x) Σy p(y)(x+y)2 - [ Σxp(x)Σyp(y) (x+y)]2 = T1 - (T2)2 T1 = Σxp(x) Σy p(y)(x2 + 2xy + y2) = Σxp(x) Σy p(y) x2 + Σxp(x) Σy p(y) 2xy + Σxp(x) Σy p(y) y2 = Σxp(x) x2 + 2 Σxp(x)x Σy p(y) y + Σy p(y) y2 = E(X2) + E(Y2) + 2 E(X)E(Y) T2 = Σxp(x)Σyp(y) (x+y) = Σxp(x)Σyp(y) x + Σxp(x)Σyp(y)y = Σxp(x) x + Σyp(y)y = E(X) + E(Y) Therefore var(Z) = T1 - (T2)2 = E(X2) + E(Y2) + 2 E(X)E(Y) - [E(X) + E(Y)]2 = E(X2) - E(X)2 + E(Y2) - E(Y)2 = var(X) + var(Y) So there is the theorem I have been looking for. Theorems 1 through 6 Theorem 1: If X and Y are independent random variables, with p(x) and p(y) arbitrary, then var(X+Y) = var(X) + var(Y). Back up and take a different pathway: E(Z) = Σx,y (x+y)p(x,y) E(Z2) = Σx,y (x+y)2p(x,y) Now recall that cov(X,Y) = Σx,y (x-μx)(y-μy) p(x,y) = Σx,y xy p(x,y) - μxΣx,yy p(x,y) - μy Σx,y x p(x,y) + μxμyΣx,y p(x,y) = Σx,y xy p(x,y) - μxΣy y Σx p(x,y) - μy Σx x Σy p(x,y) + μxμy = Σx,y xy p(x,y) - μxΣy y p(y) - μy Σx x p(x) + μxμy = Σx,y xy p(x,y) - μxμy - μy μx + μxμy = Σx,y xy p(x,y) - μxμy = E(XY) - E(X)E(Y) This is another new fact to me, so Theorem 2: cov(X,Y) = E(XY) - E(X)E(Y) = Σx,y xy p(x,y) - μxμy . Obviously if X and Y are independent, the result is 0, and if X = Y we get the usual result. But this Theorem 2 is a more general result that both of these claims. Now back up to E(Z) = Σx,y (x+y)p(x,y) E(Z2) = Σx,y (x+y)2p(x,y) Then write E(Z2) = E[(X+Y)2] = E(X2) + E(Y2) + 2E(XY) = E(X2) + E(Y2) + 2 [ cov(X,Y) + E(X)E(Y) ] E(Z) = E(X+Y) = E(X) + E(Y) Then var(Z) = E(Z2) - [E(Z)]2 = E(X2) + E(Y2) + 2 [ cov(X,Y) + E(X)E(Y) ] - { E(X) + E(Y)}2 = var(X) + var(Y) + 2 [ cov(X,Y) + E(X)E(Y) ] - 2 E(X)E(Y) = var(X) + var(Y) + 2 cov(X,Y) So another general theorem: Theorem 3: var(X+Y) = var(X) + var(Y) + 2 cov(X,Y) . Theorem 4: var(aX+bY) = a2var(X) + b2var(Y) + 2ab cov(X,Y) Now generalize for adding more terms. How would that go? Z = Σi aiXi Z2 = Σi aiXi Σj ajXj = Σij aiajXiXj E(Z) = Σi ai E(Xi) E(Z2) = Σij aiaj E(XiXj) Then var(Z) = E(Z2) - [E(Z)]2 = Σij aiaj E(XiXj) - [ Σi ai E(Xi)]2 Now write E(XiXj) = cov(Xi,Xj) + E(Xi)E(Xj) to get var(Z) = Σij aiaj { cov(Xi,Xj) + E(Xi)E(Xj)} - [ Σi ai E(Xi)]2 = Σij aiaj cov(Xi,Xj) + Σij aiajE(Xi)E(Xj) - [ Σi ai E(Xi)]2 = Σij aiaj cov(Xi,Xj) + ΣiaiE(Xi)ΣjajE(Xj) - [ Σi ai E(Xi)]2 = Σij aiaj cov(Xi,Xj) + [ΣiaiE(Xi)][ΣjajE(Xj)] - [ Σi ai E(Xi)]2 = Σij aiaj cov(Xi,Xj) Suppose we had only two Xi. This would be var(a1X1 + a2X2) = a12 cov(X1,X1) + 2 a1a2 cov(X1,X2) + a22 cov(X2,X2) = a12 var(X1) + 2 a1a2 cov(X1,X2) + a22 var(X2) and this reduces to Theorem 4. So we have an intermediate new theorem: Theorem 5. var(ΣiaiXi) = Σij aiaj cov(Xi,Xj) This also appears on the wiki variance page Good. Now break the sum into pieces Σij aiaj cov(Xi,Xj) = Σi=j aiaj cov(Xi,Xj) + Σi≠j aiaj cov(Xi,Xj) = Σi ai2 var(Xi) + Σi≠j aiaj cov(Xi,Xj) = Σi ai2 var(Xi) + 2 Σi>j aiaj cov(Xi,Xj) and this also appears on that same wiki variance page. Corollary 5a. var(ΣiXi) = Σij cov(Xi,Xj) just setting all ai = 1 Now suppose all the Xi are independent. Then we get Theorem 6. var(ΣiaiXi) = Σi ai2 var(Xi) if all Xi are independent. Corollary 6a: var(ΣiXi) = Σi var(Xi) if all Xi are independent This was the theorem I was looking for when I started looking into things. Example: Suppose you roll a die two times. Then var(X+Y) = var(X) + var(Y) = 2 var(X) And if you roll it N times, var( X+ Y + ...) = N var(X) Theorem 7 (Standard Deviation of the Mean Theorem) For N repeats of the same experiment which have outcomes x (for random variable X), if the experiments are independent (they are) then σ() = σ/ where σ() is the standard deviation of the mean and σ is the standard deviation for one experiment. The pdf shape for the single experiment can be anything. Proof: Suppose you roll the dice N times (or any other experiment, X being some outcome). Then T ≡ X1 + X2 + .....XN = the function used in Corollary 6a. so var(T) = Σi var(Xi) = Nσ2 if all the processes are the same, like rolling a die N times. Now this can be restated as T/N = (X1 + X2 + .....XN)/N = (1/N)X1 + .... ai = 1/N. Then apply Theorem 6 var(T/N) = Σi (1/N2) var(Xi) = (1/N2)Σi var(Xi) = (1/N2) N var(X) = var(X)/N = σ2/N Now T/N is the random variable for the mean value of all your samples. Call that . Then σ2() = σ2/N = the variance of the mean value of your repeated X measurements. or σ() = σ/ where σ is the SD of the "underlying experiment". QED. This theorem and its proof are based on http://en.wikipedia.org/wiki/Standard_error . Example: Suppose you roll a die N = 100 times. We know that σ2 = 2.91666 from above, so we then know that σ = 1.708. The mean is μ = 3.5. Then σ() = 1.708/ 10 = .1708 You could then add and subtract this to get your 1 std deviation range for the result. Comment: Nothing has been said so far about any kind of "distribution". We have not specified anything to be a Gaussian distribution, for example. We assume that our random variable has SOME distribution, and therefore it has some variance and some σ. You still have variance and σ if the single event distribution is not symmetrical about its mean. In the special case that you have a symmetric distribution for Xi, then I think T/N will also have a symmetric distribution. If the distribution is Gaussian, then you get 95% CI if you use 1.96 σ() on either side. "observation from a population" = observation of face-up number on a rolled die. Comment: In HAS-BLED, the thing that corresponds to rolling the die is more like flipping a coin. For each of the N people in a HAS-BLED bin, there is a probability of a bleed which we call p. I have now exhausted the wiki page on "standard error". C. The Binomial Distribution and the Bernoulli Trial Now I think the binomial distribution is going to get involved for HAS-BLED. For each person in the group of N = 500, either he bleeds or he doesn't bleed. Suppose p = probability of a person bleeding. Note added: Flipping a weighted coin represents any Bernoulli Trial type process with A or B outcome. It is the repeated Bernoulli Trial which results in the binomial distribution. Now digress to flipping N weighted coins where p = probability of a head. To get n heads, the probability is (N,n) pmqN-m . In other words probN(n heads) = (N,n) pnqN-n If p = 1/2, you get probN(n heads) = (N,n) (1/2)n = [ N!/(n! (N-n)! ] 2-n The curve is symmetric if p = 1/2 because you can see that probN(n) = probN(N-n). Here are some plots N = 5 N = 10 N = 20 But if p = 0.25, things are different, the distribution is not symmetric: N = 5 N = 10, p = 1/4 N = 20, p = 1/4 Now the distribution is NOT symmetrical. For small N it gets smashed into the left wall. However, even if p ≠ q, for large N the distribution becomes close to symmetrical again. N = 100, p=1/4 The mean value is Np = 100*1/4 = 25. Here is some wiki stuff I could derive for this distribution, and we now rename N to be n, and n to be k Now the roll of the die here is the test of one patient for bleed or no bleed. Let's do this no-bleed = 0 bleed = 1 // like a coin toss. The PMF for this underlying experiment looks like this The Bernoulli Trial Picture What is σ for this thing? Well var = E(X2) - [E(X)]2 E(X) = p, probability of a bleed p(x) = p for x = 1 // this is the pmf = q for x = 0 [E(X)] = Σx x p(x) = 0*q + 1*p = p E(X2) = Σx x2 p(x) = 02 *q + 12*p = p var = p - p2 = p(1-p) = pq = σ2 If you repeat this experiment (like coin toss) N time. Now we want to know the standard deviation of the mean for doing this N times. The answer is σ() = σ/ = / = The mean here is = E(x) = p, as just shown above. So I have just proven: Theorem 8: (Bernoulli Trial SD of the Mean Theorem) . The standard deviation of the mean where the underlying process is a simple Bernoulli trial is given by σ() = . Proof. We know from Theorem 7 that in general σ() = σ/ where σ is for the single experiment, For the Bernoulli experment, as shown above σ = . QED. Here is my older statement. Theorem 8: (Bernoulli Trial SD of the Mean Theorem) For an experiment with a binary outcome with X = 1 or X = 0, we assume there is some pmf which can be characterized entirely by prob(1) = p. If this experiment is repeated N times, the mean value of the outcome will be p if N is a very large number. But for finite N, there is some distribution for the outcome (which is the binomial distribution) and we find that σ() = where σ() is the standard deviation of the mean p. In some sense, the mean you measure with finite N will lie in the distribution and the half width of the distribution in some sense is this σ(). Even for very small N, this is all true, but for such small N the binomial distribution is not symmetric about the mean, as my plots above show, so you have to think about how to construct something like a confidence interval for small N. I first say the claim of this theorem on the wiki page for margin of error, where it appears this way: Go back now to probN(n heads) = (N,n) pnqN-n (a) I try but fail to show that Binomial Distribution approaches a Gaussian for large N. In what sense does this approach a Gaussian for large N? I sort of see it in my graphs for p = 1/2. Stirling tells us that So we can certainly claim that N! ≈ (N/e)N I suppose for some reasonable p value, everything of interest is large so then ≈ We end up with one factor of on the bottom. Power of e are these e–N * en * eN-n = 1 so they all go away Then we have NN * n-n * (N-n)-(N-n) I guess we now assume N >> n of interest so ≈ NN * n-n * (N)-(N-n) = n-n * (N)-(-n) = n-n Nn = (N/n)n This seems wrong, but continue. Then ≈ (1/) * (N/n)n and then probN(n heads) ≈ (1/) * (N/n)n pnqN-n The mean value is x0 = Np, and let x = np. Then n = x/p so probN(n heads) ≈ (1/) * (N/n)n pnqN-n = (1/) * (Np/x)(x/p) p(x/p)qN-x/p = (1/) * (x0/x)(x/p) p(x/p)qN-x/p I don't see how this would ever have the form exp( - (x-x0)2/A ). Keep going = (1/) * (x0/x)(x/p) p(x/p)q(Np-x)/p = (1/) * (x0/x)(x/p) p(x/p)q(x0-x)/p = (1/) * e(x/p)ln(x0/x) eln(x/p)eln[(x0-x)/p] = (1/) * e(x/p)[ln(x0)-ln(x)] eln(x/p)eln[(x0-x)/p] = (1/) * e(x/p)[ln(x0)-ln(x)] e[lnx - lnp] e[ln(x0-x)-lnp] = (1/) * e(x/p)[ln(x0)-ln(x)] elnx eln(x0-x) e-2lnp For large x the leading term dominates since it has an exposed x. (b) Various Binomial Notations and Access from Excel or Maple I need to find a source on this issue. What is b(n,p,j) ? Here is one notation from wiki on binomial, so yes, in my notation the variables are N, n and p so probN(n heads) = (N,n) pnqN-n = f(n; N,p) Here is another notation from my 12 chapter PDF book, and I would have (N,n) pnqN-n = b(N,p,n) Excel has this So the last argument is just a flag which I would set false. Then I guess BINOMDIST(n,N,p,FALSE) = b(N,p,n) In maple we have binomial(n r) = n!/r!/(n-r)! which is just the coefficient part which Excel calls COMBIN(n,r). Maple doesn't have the full function, you have to make it yourself like so b := (n,p,j) -> binomial(n j) pj (1-p)n-j So, what does my big PDF book have to say on this? μ = 0 and σ = 1 of the above So we now have a term for what I am doing "Bernoulli Trials". Well OK this thing Sn* is set up to have a mean of 0, AND it is scaled in a certain way. (c) The Central Limit Theorem and my Second Failed Attempt to show Bin → Gaussian Now here is the big connection claim: The author however does not prove this. I would write it out as b(n,p,<x0 + x>) → exp(-x2/2) / as n→∞ Now back up a little and go back to b := (n,p,j) -> binomial(n j) pj (1-p)n-j b(n,p,j) When you start, n and j are integers, but the gamma functions continue this thing to general values. So perhaps of interest is b(n,p,x) where we let x be a continuum Now you know that the SD of this thing is σ = and the mean is μ = np. So the claim is b(n,p,μ+xσ) → exp(-x2/2) / as n→∞ Write this as b(n,p,y) → exp(-x2/2) / as n→∞ y = μ + σx At x = 0, y = μ and the left side has it max. RHS also max at x = 0. Good. At x = 1, y = μ+σ and we are one SD off the peak, more or less, and so is RHS. Let's try one more time on this. LHS = σ b(n,p,y) = σ (n,y) py qn-y = σ py qn-y Now look at our binomial plot for large n = 1000 In our range of non-zero-ness, we are going to have y large as well as n large. And n-y ≈ n/2 so it is also large. Then go ahead on all 3 with Stirling, n! = nn e-n = = = Then we have LHS = σ py qn-y y = μ + σx Now the big question: what do you do next? All I can think to do is expo everything. num = e(n+1/2)ln(n) σ eyln(p) e(n-y)ln(q) den = e(y+1/2)ln(y) e(n-y+1/2)ln(n-y) Then can combine into a single exp: LHS = exp{ (n+1/2)ln(n) + yln(p) + (n-y)ln(q) - (y+1/2)ln(y)- (n-y+1/2)ln(n-y) } If n, y and n-y are all large, as I have assumed, then the 1/2 contributions should not matter, then LHS = exp{ n ln(n) + yln(p) + (n-y)ln(q) - y ln(y) - (n-y) ln(n-y) } As before, I am completely stumped, no idea what to do next. How get this to be exp(-x2/2) / ? It does seem that the terms like y ln(p) and (n-y)ln(q) are smaller than the other terms, so try LHS ≈ exp{ n ln(n) - y ln(y) - (n-y) ln(n-y) } y = μ + σx and at least now things don't depend on p and q. Maybe regroup terms, LHS ≈ exp{ n ln(n) - n ln(n-y) - y ln(y) + y ln(n-y) } LHS ≈ exp{ n [ ln(n) - ln(n-y)] - y [ ln(y) - ln(n-y)] } LHS ≈ exp{ n [ ln(n) - ln(n-y)] - y [ ln(y) - ln(n-y)] } ln(n) - ln(n-y) = ln = - ln = - ln (1 - y/n) ln(y) - ln(n-y) = ln = - ln = - ln(n/y-1) Then we have LHS ≈ exp{ n [ - ln (1 - y/n)] - y [ - ln(n/y-1] } LHS ≈ exp{ - n ln (1 - y/n)] +y ln(n/y-1] } Once again, I hit a dead end. I will need to find a proof. The Central Limit Theorem says that you end up with a normal distribution when you add a bunch of random variables no matter what you had for p(x) for that variable! So it is a very general theorem, and proofs I am looking at are rather involved affairs. Wiki has a "simple proof". (d) Application of Bernoulli/Binomial to HAS-BLED Status: For the HAS-BLED study we have this The numbers are the left are number of people in each bin. These numbers are all > 20 except the last one, but I think it would be just find to do 1.96 SD and just assume normal distribution for each. The key idea is this For example, try applying this to the bin with 22 people in it. We have p = .091 rate of bleeds n = 22 Recall, σ() = σ/ = / = My conclusion is that rate = .091 ± .120 ≈ 0 to .21 Then the measured rate is 9.1%, and the actual rate with 95% confidence is in the range 0% to 21%. Now let's try the next bin up: Then the rate is .089 ± .042 = (.047, ,131) = 4.7% to 13.1% I am now ready to do this in Excel. It took me all day to get ready?