odd Solutions
2.1
a) case: m is middle value; even case: middle values are m d and m + d for some d, so m is the median. b) Suppose L values are less than m and E values equal to m. Then there are also L values greater than m. Since m d + m + d = 2m, the sum of the values is xi = xi <m xi + xi >m xi + xi =m xi = 2Lm + Em = (2L + E)m, so the mean is m. c) {5, 1, 0, 2, 4} has both mean and median equal to 0 but is not symetric.
2.6 SS (c) =
2(xi c) = 2 ((
xi ) nc)) = 0 c =
1 n
(xi ) = x.
2.7 ve 0s and ve 10s. 2.8 Have all the values be the same. Variance is zero. 2.9
> award <- scan() 1: 37 60 75 115 135 140 149 150 238 17: 1050 1100 1139 1150 1200 1200 1250 28: Read 27 items > sum(award^2) [1] 24657511 > mean(award) + c(-1,1) * 2 * sd(award) [1] -466.4174 1961.1581 290 340 410 600 750 1576 1700 1825 2000 750 750
1015
Index
3.3 Let X be any continuous random variable with range [1, 2]. Let Y be determined as follows: ip a coin. If it is heads, Y = 0. Otherwise Y = X. 3.5 This is equivalent to looking for non-negative integer solutions to x + y + z = 12. By analogy to Example 3.2.5, there are 12 + 1 (122 ) such solutions. 2 3.6
> (choose(13,2)*choose(4,2)*choose(4,2)*choose(11,1)*choose(4,1)) / choose(52,5) [1] 0.04753902
3.7
> (choose(13,1)*choose(4,3)*choose(12,1)*choose(4,2)) / choose(52,5) [1] 0.001440576
3.8
> (choose(13,1)*choose(4,3)*choose(12,2)*choose(4,1)*choose(4,1)) / choose(52,5) [1] 0.02112845
3.9
> > > > > > p<-rep(NA,4) p[1] <- choose(4,1) * choose(13,5) / choose(52,5) p[2] <- choose(4,2) * ( choose(26,5) - (choose(13,5)+choose(13,5))) / choose(52,5) p[4] <- choose(4,1) * choose(13,2) * 13 * 13 * 13 / choose(52,5) p[3] <- 1 - sum(p[-3]) # sum of all probabilities must be 1 rbind(1:4,p) 1.000000000 2.0000000 3.0000000 4.0000000 p 0.001980792 0.1459184 0.5883553 0.2637455
3.10
1016
c 2007 Randall Pruim ([e-mail])
Index
> birthdayprob <- function(n) { # calculates prob of shared birthday + last <- 366 - n; + 1 - ( prod(seq(last,365)) / 365^n ); + } > > birthdayprob(10); [1] 0.1169482 > cbind(20:25,sapply(20:25,birthdayprob)) [,1] [,2] [1,] 20 0.4114384 [2,] 21 0.4436883 [3,] 22 0.4756953 [4,] 23 0.5072972 [5,] 24 0.5383443 [6,] 25 0.5686997
3.13 P(AB C) = P(A)+P(B)+P(C)P(AB)P(AC)P(B C)+P(AB C) 3.14
P(Bad) = 3/18; P(AL 1) = 8/18
2/18 8/18 ;
P(Bad | AL 1) = 2/8 =
P(AL 1 | Bad) = 2/3 =
2/18 3/18
3.17 If A and B are independent, then
P(A | B c ) = = P(A B c ) P(A) P(A B) = c) P(B 1 P(B) P(A) P(A) P(B) P(A)(1 P(B)) = , 1 P(B) 1 P(B)
so A and B c are independent. Part (b) follows from this (or can be proven analogously).
3.18 P(S) =
11 = 0.12088 91
P(A) =
50 91
15:21 -- December 4, 2007
1017
Index 6 = 0.0659 91 11 50 P(A) P(S) = = 0.0664 91 91 P(S and A) = P(S|A) =
6 91 50 91
= 0.12
3.19
P(D) = P(D and AA) + P(D and Aa) + P(D and aa) = P(AA) P(D|AA) + P(Aa) P(D|Aa) + P(aa) P(D|aa) = (.25)(.01) + (.5)(.05) + (.25)(.5) = .1525
P(C+ and T +) P(T +)
3.20 Goal: Compute P(C + |T +) =
Data:
P(C+) = 2/3 [based on prior information] P(T + |C+) = .7 [sensitivity] Results: P(T |C) = .9 [specicity] P(T + and C+) = P(C+) P(T + |C+) = (2/3)(0.7) = 14/30 P(T + and C) = P(C) P(T + |C) = (1/3)(0.1) = 1/30 So P(C + |T +) =
14/30 15/30
= 14/15 = 0.9333
For part b, compute P(C + |T ) instead. Same methods work. (Answer is 6/15)
3.21 P(All 5 numbers the same) = 1 1 1 1 6 6 6 3.23 Let p = P(D), then P(D | +) =
1 6
Now we can let R do the work.
P(D +) P(D) P(+ | D) p .98 = = P(+) P(D +) + P(H +) p .98 + (1 p) .01
1018
c 2007 Randall Pruim ([e-mail])
Index
> p <- c(.01,.10) > ( p * .98 ) / ( p * .98 + (1-p)*0.01 ); [1] 0.4974619 0.9158879
3.24
> factorial(10) / ( factorial(3) * factorial(3) * factorial(2) ) [1] 50400
3.25
> choose(90,4)/choose(100,4) # prob only good ones selected [1] 0.6516305 > 1 - choose(90,4)/choose(100,4) # lot is rejected [1] 0.3483695 > f <- function(x) { 1 - choose(100-x,4)/choose(100,4) } > myplot <- xyplot(sapply(10:100,f) ~ 10:100, type=l, + xlab="number of defective parts", + ylab="probability of rejecting", + lwd=2, + col="navy");
1.0
probability of rejecting
0.8
0.6
0.4
20
40
60
80
100
number of defective parts
3.26
15:21 -- December 4, 2007
1019
Index
> 1/8 + # HHH-+ 1/16 + # or THHH+ 1/16 # or -THHH [1] 0.25
3.29 The wines can be presented in 24 = 4 3 2 1 dierent orders, so the probability of guessing correctly is 1/24. 3.31
> 8*5*4 / choose(17,3) # 1 sock of each kinds means no pairs [1] 0.2352941 > 1 - (8*5*4 / choose(17,3)) # so this is prob of getting a pair [1] 0.7647059 # or do it this way > ( choose(8,2)*9 + choose(5,2) * 12 + choose(4,2) * 13 + + choose(8,3) + choose(5,3) + choose(4,3) ) / choose(17,3) [1] 0.7647059
3.32 Let F be the event that a double-headed (fake) coin is selected, then
a) P(F | H) =
P(F
and
P(H)
H)
=
1 1 3 1 1 0+ 3 2 + 1 1 3
=
2 3
1 0+ 1 1 + 1 1 3 3 4 3 1 0+ 1 1 + 1 1 3 3 2 3
b) P(H | H) = c) P(F | HH) =
P(H and H) P(H) P(F
=
P(HH) P(H)
=
=
5 6
and HH) P(HH)
=
1 1 3 0+ 1 1 + 1 1 3 4 3
=
4 5
= 0.80
3.33
> 1 - dbinom(0,4,1/6) # P(at least one 6 in 4 tries) [1] 0.5177469 > pgeom(3,1/6) # P(fail at most 3 times before getting a 6) [1] 0.5177469 > 1 - dbinom(0,24,1/36) # P(at least one double 6 in 24 tries)
1020
c 2007 Randall Pruim ([e-mail])
Index
[1] 0.4914039 > pgeom(23,1/36) [1] 0.4914039
# P(fail at most 23 times before getting double 6)
3.34 3.38
value of X probability
1 .4
2 .3
3 .2
4 .1
> dbinom(10,10,.8) [1] 0.1073742 > 1 - pnbinom(4,10,.80) [1] 0.1298396 > pbinom(9,14,.8) [1] 0.1298396 > > 1 - pnbinom(4,10,.70) [1] 0.4157988 > pbinom(9,14,.7) [1] 0.4157988
# make 10 straight # at least 5 misses = at least 15 shots # 9 or fewer makes in 14 tries -> at least 15 shots
# at least 5 misses = at least 15 shots # 9 or fewer makes in 14 tries -> at least 15 shots
3.39
> pp <- dbinom(5,5,.80); pp # a) prob make 5 straight (= success) [1] 0.32768 > dgeom(1,pp) # b) succeed with 1 miss (correct answer) [1] 0.2203058 > 0.20 * pp # miss first then make 5 straight (_not_ the answer) [1] 0.065536 > > 1 - pgeom(1,pp) # c) miss more than one shot before success [1] 0.4520142 > probs <- dgeom(0:15,pp) # d) > myplot <- xyplot(probs~0:15,main="Freddie",xlab="misses",ylab="probability"); ################################################################# > pp <- dbinom(5,5,.70) # a) prob make 5 straight (= success) > pp
15:21 -- December 4, 2007
1021
Index
[1] 0.16807 > > dgeom(1,pp) # b) succeed with 1 miss (correct answer) [1] 0.1398225 > 0.20 * pp # miss first then make 5 straight (_not_ the answer) [1] 0.033614 > > 1 - pgeom(1,pp) # c) miss more than one shot before success [1] 0.6921075 > probs <- dgeom(0:15,pp) # d) > myplot <- xyplot(probs~0:15,main="Frank",xlab="misses",ylab="probability")
Freddie
Frank
0.3
0.15
probability
probability
0.2
0.10
0.1
0.05
0.0 0 5
10
15
0
5
10
15
misses
misses
3.40
> 1 [1] > 1 [1] > 1 [1] - pbinom(11,20,.25) # 11 or fewer correct fails 0.0009353916 - pbinom(11,20,1/3) # 11 or fewer correct fails 0.01297330 - pbinom(11,20, .5 + .4 * 1/3 + .1 * 1/4) 0.7865717
3.42
a) P(X k) =
x=k (1
)x =
(1)k 1(1)
= (1 )k .
P(X=x and Xk) P(Xk)
b) Assuming x k, P(X = x | X k) = P(X = x k)
=
(1)x (1)k
= (1 )x k =
1022
c 2007 Randall Pruim ([e-mail])
Index
3.43
> binom.test(8,50,.25) Exact binomial test data: 8 and 50 number of successes = 8, number of trials = 50, p-value = 0.1899 alternative hypothesis: true probability of success is not equal to 0.25 95 percent confidence interval: 0.07170077 0.29112631 sample estimates: probability of success 0.16
3.45
> binom.test(428,428+152,.75); Exact binomial test data: 428 and 428 + 152 number of successes = 428, number of trials = 580, p-value = 0.5022 alternative hypothesis: true probability of success is not equal to 0.75 95 percent confidence interval: 0.7001255 0.7732932 sample estimates: probability of success 0.737931
3.46
> binom.test(10,10+17); Exact binomial test
15:21 -- December 4, 2007
1023
Index
data: 10 and 10 + 17 number of successes = 10, number of trials = 27, p-value = 0.2478 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.1940072 0.5763204 sample estimates: probability of success 0.3703704 > binom.test(36,7+36); Exact binomial test data: 36 and 7 + 36 number of successes = 36, number of trials = 43, p-value = 8.963e-06 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.6929891 0.9319479 sample estimates: probability of success 0.8372093
In the rst case, the p-value is large, so there is not enough evidence to reject the null hypothesis (that the wasps have no preference between the two options). It could be the that the dierence in wasp behavior is simply due to random chance. In the second case, the p-value is very small, so we reject the null hypothesis (that wasps have no preference between the two options).
3.47
> qbinom(.975,200,.5) [1] 114 > qbinom(.025,200,.5) [1] 86 > pbinom(85:86,200,.5) [1] 0.02001860 0.02798287 > 1 - pbinom(114:115,200,.5) [1] 0.02001860 0.01406270
1024
c 2007 Randall Pruim ([e-mail])
Index
> # define a function to calculate power for given sample size. > power = function(size,null=.50,alt=.55){ + leftCritical <- -1 + qbinom(.025, round(size),null) + rightCritical <- 1 + qbinom(.975, round(size),null) + leftPower <- pbinom( leftCritical, round(size),alt); + rightPower <- 1 - pbinom( rightCritical - 1, round(size),alt); + + result <- leftPower + rightPower; + + # the rest of this allows for more verbose output + attr(result,"power") = leftPower + rightPower; + attr(result,"null") = null; + attr(result,"alt") = alt; + attr(result,"leftCritical") = leftCritical; + attr(result,"rightCritical") = rightCritical; + attr(result,"leftPower") = leftPower; + attr(result,"rightPower") = rightPower; + + return (result); + } > # just the power values: > as.numeric(power(200)); [1] 0.2619829 > as.numeric(power(400)); [1] 0.4806564 > # more verbose output: > power(200); [1] 0.2619829 attr(,"power") [1] 0.2619829 attr(,"null") [1] 0.5 attr(,"alt") [1] 0.55 attr(,"leftCritical") [1] 85
15:21 -- December 4, 2007
1025
Index
attr(,"rightCritical") [1] 115 attr(,"leftPower") [1] 0.0002581707 attr(,"rightPower") [1] 0.2617247 > > power(400); [1] 0.4806564 attr(,"power") [1] 0.4806564 attr(,"null") [1] 0.5 attr(,"alt") [1] 0.55 attr(,"leftCritical") [1] 179 attr(,"rightCritical") [1] 221 attr(,"leftPower") [1] 2.478852e-05 attr(,"rightPower") [1] 0.4806316 > # find sample size with 90% power > uniroot(function(size){ as.numeric(power(size)) -.90} ,c(400,5000)) $root [1] 1075.5 $f.root [1] -0.001766926 $iter [1] 26 $estim.prec [1] 6.103516e-05
1026
c 2007 Randall Pruim ([e-mail])
Index
3.48 With a sample of size 200, the power of a test of the null hypothesis that the coin is fair against an alternative where the coin has a probability of .55 of coming up heads is approximately 26%. With a sample size of 400, the power increases to approximately 48%. To achieve 90% power, we need a sample size of approximately 1075. 3.49 It turns out that some of these tests dont actually have = .05. This is because the binomial distribution is discrete, and for some combinations of n and , we cant approximate = 0.05 as accurately as for others. So some of these tests might have an actual value that is less than 0.05, which makes the test less powerful. 3.50 Many examples exist. t(x) = x2 works for any random variable that is not constant. 3.51
The rst equality is the denition of expected value. 1 2: we omit the rst term from the sum since it is 0. 2 3: factor out n x
n x
from binomial coecient. Note that = n! n (n 1)! n n1 = = x!(n x)! x (x 1)!(n x)! x x1 .
3 4: change of variable, y = x 1. 4 5: sum is now all the probabilities of a binomial random variable with parameters n 1 and , so the sum is 1.
3.52
(3.5) denition. (3.6) expand the square. (3.7) distribute over sum and split into multiple summations. (3.8) denition of expected value and factoring of constants.
15:21 -- December 4, 2007
1027
Index (3.9) denition of expected value and sum of probabilities must be 1. (3.10) simple algebra and renaming.
a2 E(X 2 )
3.53 Var(aX + b) = E((aX + b)2 ) E(aX + b)2 = E(a2 X 2 + 2abX + b2 ) (a E(X) + b)2 =
+ 2ab E(X) + b2 a2 E(X)2 2ab E(X) b2 = a2 Var(X).
3.55
> vals <- 0:50; > probs.x <- dnbinom(vals,3,.5); > probs.y <- dnbinom(vals,3,.2); > var.x <- sum(vals^2 * probs.x) - sum(vals [1] 6 > var.y <- sum(vals^2 * probs.y) - sum(vals [1] 58.67976 > # better approximation using more terms: > vals <- 0:500; > probs.x <- dnbinom(vals,3,.5); > probs.y <- dnbinom(vals,3,.2); > var.x <- sum(vals^2 * probs.x) - sum(vals [1] 6 > var.y <- sum(vals^2 * probs.y) - sum(vals [1] 60
* probs.x)^2; var.x; * probs.y)^2; var.y;
* probs.x)^2; var.x; * probs.y)^2; var.y;
3.56 E(X) =
E(X 2 ) =
n 1 1 n+1 k=1 k n = n (1 + 2 + + n) = 2 . n 1 1 n(n+1)(2n+1) 21 2 2 k=1 k n = n (1 + 2 + + n ) = n 6
=
(n+1)(2n+1) . 6
3.57
Part (a): This is just a uniform random variable, so the mean is
n+1 2
= 128.
1 2 4 8 16 32 64 Part (b): E(X) = 1 255 +2 255 +3 255 +4 255 +5 255 +6 255 +7 255 +8 128 7.03. 255
1028
c 2007 Randall Pruim ([e-mail])
Index
# expected value of binary search > vals <- 1:8; > probs <- 2^(vals-1) / 255; > sum(vals*probs); [1] 7.031373 > # expected value of linear search > vals <- 1:255; > probs <- rep(1/255,255) > sum(vals*probs); [1] 128
3.63 For every x, y, and z,
P(X = x and Y = y | Z = z) = P(X = x | Z = z) P(Y = y | Z = z) So once we know the value of Z, knowing the value of X or Y gives no additional information about the other.
3.65 We have already handled the case where r = 1. If r > 1, then we can express X NBinom(r, ) as X = X1 + X2 + Xr where Xi NBinom(1, ) = Geom(), and the Xi s are independent. We already know that E(Xi ) = 1/ 1, so
r
E(X) =
i=1
E(Xi ) = r(
1 r 1) = r
3.70
> f <- function(prob) { prob + prob * (1-prob)^3/(1-(1-prob)^2) } > g <- function(prob) {f(prob) - .50} > uniroot(g,c(.20,.5))$root when # g=0, f=.5 [1] 0.2929236
3.72
15:21 -- December 4, 2007
1029
Index
> d <- c(9,13,14,9); dim(d) = c(2,2); d [,1] [,2] [1,] 9 14 [2,] 13 9 > fisher.test(d) Fishers Exact Test for Count Data data: d p-value = 0.238 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.1141432 1.7066640 sample estimates: odds ratio 0.4533593 > phyper(9,23,22,22) [1] 0.1490436
3.74
> d <- c(61,103,69,44); dim(d) = c(2,2); d [,1] [,2] [1,] 61 69 [2,] 103 44 > fisher.test(d) Fishers Exact Test for Count Data data: d p-value = 0.0001362 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.2236475 0.6366559 sample estimates: odds ratio 0.3790412
1030
c 2007 Randall Pruim ([e-mail])
Index
> phyper(61,61+103,44+69,61+69) [1] 7.214999e-05
3.75 The sum is telescoping and simplies to log10 (10) log10 (1) = 1 0 = 1. 4.4 The cdf of an exponential distribution is F (x) = 1ex , so the p-quantile q satises
1 p = eq p = 1 eq
ln(1 p) = q ln(1 p) = q ,.
So the median is rst, second, and third quartiles are tively.
ln(.75) ln(.5) , ,
and
ln(.25) ,
respec-
4.5 limx F (x) = 0, limx F (x) = 1, and
so F is a legitimate cdf.
> qcauchy(.9) [1] 3.077684
d dx F (x)
= f (x) =
1 (x2 +1)
0 for all x,
4.6
a) P(X 1) = F (1) = 1/4 b) P(.5 X 1) = F (1) F (.5) = 1/4 1/16 = 3/16 c) P(X > 1.5) = 1 F (1.5) = 1 9/16 = 7/16 d) The median of X: m2 /4 = 1/2 = m = 2 e) The pdf of X: f (x) = x/2
15:21 -- December 4, 2007
1031
Index
> f <- function(x) {x/2;} # define pdf > integrate(f,lower=0,upper=2); # check it is a pdf 1 with absolute error < 1.1e-14 > xf <- function(x) x^2/2; > integrate(xf,lower=0,upper=2); # expected value 1.333333 with absolute error < 1.5e-14 > xxf <- function(x) x^3/2; > integrate(xxf,lower=0,upper=2); # E(X^2) 2 with absolute error < 2.2e-14 > # compute the variance using E(X^2) - E(X)^2 > integrate(xxf,lower=0,upper=2)$value + (integrate(xf,lower=0,upper=2)$value)^2 [1] 0.2222222
4.7 This can all be done by hand to give exact values, but here is code to show how to get R to do the numerical calculations.
> f <- function(x) {1/x^4}; > k <- 1 / integrate(f,1,Inf)$value; k; [1] 3 > f <- function(x) {k / x^4}; > integrate(f,1,Inf); 1 with absolute error < 1.1e-14 > integrate(f,2,3); 0.08796296 with absolute error < 9.8e-16 > xf <- function(x) { x * f(x) }; > # find the median > g <- function(x) { integrate(f,1,x)$value - .5 } > uniroot(g,c(1,10))$root; [1] 1.259922 > > Ex <- integrate(xf,1,Inf)$value; Ex; # E(X) [1] 1.5 > xxf <- function(x) { x^2 * f(x) }; > Exx <- integrate(xxf,1,Inf)$value; Exx; # E(X^2) [1] 3
1032
c 2007 Randall Pruim ([e-mail])
Index
> Exx - (Ex)^2 [1] 0.75 > sqrt(Exx - (Ex)^2) [1] 0.8660254 # variance # st dev
4.8
> dpois(0,6/3); [1] 0.1353353 > 1-pexp(20/60,rate=6); [1] 0.1353353 > dpois(2,6/3); [1] 0.2706706 > 1-ppois(6,6); [1] 0.3936972 > dpois(6,6); [1] 0.1606231 > ppois(5,6); [1] 0.4456796 > ppois(30,24) - ppois(19,24); [1] 0.7238911
# 0 customers in 1/3 hour # first customer after 1/3 hour # 2 customers in 1/3 hour # more than 6 in 1 hour # exactly 6 in 1 hour # less than 6 in 1 hour # 20 to 30 customers in 4 hours
The Poisson model might not be good because (1) arrivals may not be independent (people come in groups or in response to something), and (2) the rate may not be constant over the time period.
15:21 -- December 4, 2007
1033
Index
4.10
MX (t) = E(etX )
n
=
x=1
etx
n
1 n
=
1 n
(et )x
x=1
1 et (et )n+1 = n (1 et ) et 1 ent = n (1 et )
4.11 4.12 Here is Mathematica code for this problem.
In[4]:= Integrate[E^(t*y)*y, {y, 0, 1}] + Integrate[E^(t*y)*(2 - y), {y, 1, 2}] Out[4]= (1 + e^t (-1 + t))/t^2 + (e^t (-1 + e^t - t))/t^2 In[5]:= Integrate[E^(t*y)*y, {y, 0, 1}] Out[5]= (1 + e^t (-1 + t))/t^2 In[6]:= Integrate[E^(t*y)*(2 - y), {y, 1, 2}] Out[6]= (e^t (-1 + e^t - t))/t^2 In[7]:= Together[ Integrate[E^(t*y)*y, {y, 0, 1}] + Integrate[E^(t*y)*(2 - y), {y, 1, 2}]]
1034
c 2007 Randall Pruim ([e-mail])
Index
Out[7]= (-1 + e^t)^2/t^2
4.13 MX (t) = e+e . 4.14 Uses integration by parts. 4.19 From Mathematica:
In[16]:= D[(1 - alpha*t)^(-k), t] Out[16]= alpha k (1 - alpha t)^(-1 - k) In[17]:= D[D[(1 - alpha*t)^(-k), t], t]
t
Out[17]= -alpha^2 (-1 - k) k (1 - alpha t)^(-2 - k) So the mean is k and the variance is (k + 1)k2 k 2 2 = k2 .
4.23
a) Binom(10, 1/2). b) Norm(1, 1). c) Exp( = 1/2). d) Gamma( = 3, = 1/2) = Gamma( = 3, = 2).
4.25
> pexp(1+1) - pexp(1-1); [1] 0.8646647 > pexp(1/2 + 1/2,rate=2) - pexp(1/2 - 1/2,rate=2) [1] 0.8646647 > punif(.5+ sqrt(1/12),0,1) - punif(.5 - sqrt(1/12),0,1)
15:21 -- December 4, 2007
1035
Index
[1] 0.5773503 > s <- sqrt(8/(6^2*7)) > pbeta(2/6 + s,shape1=2,shape2=4) - pbeta(2/6 - s,shape1=2,shape2=4) [1] 0.652183
Things to notice: Both exponential distributions give the same answer. This is because the exponential distributions are all constant multiples of each other. Compared to 68% for the normal distribution, the beta and uniform distribution are lower (not as much piled in the middle), and the exponential distribution is higher.
4.26
# a = shape; b = scale (acts like 1/lambda from exponential dist) > > a <- 2; b <- 3; > m <- b * gamma(1 + 1/a); m # mean [1] 2.658681 > v <- b^2 * ( gamma(1 + 2/a) - gamma(1+1/a)^2 ); v; # var [1] 1.931417 > qweibull(.5,a,b) # median [1] 2.497664 > pweibull(m,a,b) # less than mean [1] 0.5440619 > pweibull(6,a,b) - pweibull(1.5,a,b) # in a range [1] 0.7604851 > pweibull(m+sqrt(v),a,b) - pweibull(m - sqrt(v) ,a,b) # near mean [1] 0.6743336 > # with roles of parameters reversed > > a <- 3; b <- 2; > m <- b * gamma(1 + 1/a); m # mean [1] 1.785959 > v <- b^2 * ( gamma(1 + 2/a) - gamma(1+1/a)^2 ); v; # var [1] 0.4213315 > qweibull(.5,a,b) # median
1036
c 2007 Randall Pruim ([e-mail])
Index
[1] 1.769994 > pweibull(m,a,b) [1] 0.5093739 > pweibull(6,a,b) - pweibull(1.5,a,b) [1] 0.655816 > pweibull(m+sqrt(v),a,b) - pweibull(m - sqrt(v) ,a,b) [1] 0.6677128 >
# less than mean # in a range # near mean
4.27
> m <- 5/(2+5); m; [1] 0.7142857 > v <- 5 * 2 / ( (5+2)^2 * (5 + 2 + 1) ) > qbeta(.5,5,2); [1] 0.73555 > pbeta(m,5,2); [1] 0.4515550 > pbeta(.4,5,2) - pbeta(.2,5,2) [1] 0.03936 > pbeta(m+sqrt(v), 5, 2) - pbeta(m-sqrt(v), 5, 2) [1] 0.662012 # mean # var # median # less than mean # in a range # near mean
5.1 E(X w ) = E( (wi Xi )) =
wi = 1.
E(wi Xi ) =
wi E(Xi ) =
wi = precisely when
5.2 SE = 10/ 16 = 2.5. X Norm(0, SE). P(|X | < 2) =
> pnorm(2,sd=2.5) - pnorm(-2,2.5) [1] 0.7881412
5.3
> x <- c(1,2,4,4,9) > mu <- sum(x * .2); mu [1] 4 # population mean
15:21 -- December 4, 2007
1037
Index
> v <- sum(x^2 *.2) - mu^2; v # population variance [1] 7.6 > pairsums <- outer(x,x,"+") # compute 25 sums > pairmeans <- pairsums/2; > # sampling distribution with SRS > srs.means <- as.vector(pairmeans[lower.tri(pairmeans)]); srs.means; [1] 1.5 2.5 2.5 5.0 3.0 3.0 5.5 4.0 6.5 6.5 > iid.means <- as.vector(pairmeans); iid.means; [1] 1.0 1.5 2.5 2.5 5.0 1.5 2.0 3.0 3.0 5.5 2.5 3.0 4.0 4.0 6.5 2.5 3.0 4.0 4.0 [20] 6.5 5.0 5.5 6.5 6.5 9.0 > > srs.mean <- sum(srs.means * .1); srs.mean [1] 4 > srs.var <- sum(srs.means^2 * .1) - srs.mean^2; srs.var; [1] 2.85 > v/2 * (5-2) / (5-1) [1] 2.85 > sqrt(v/2 * (5-2) / (5-1)) [1] 1.688194 > > var(srs.means) # INCORRECT variance [1] 3.166667 > # sampling distribution with iid sample > iid.mean <- sum(iid.means * .04); iid.mean; [1] 4 > iid.var <- sum(iid.means^2 * .04) - iid.mean^2; iid.var; [1] 3.8 > v/2 [1] 3.8 > sqrt(v/2) [1] 1.949359 > > var(iid.means) # INCORRECT variance [1] 3.958333
5.4
1038
c 2007 Randall Pruim ([e-mail])
Index a) Let T be the event that a subject would answer true to question A. Let A be the even that a subject is asked question A. Then = P(subject answers true) = P(T | A) P(A) + P(T c | Ac ) P(Ac ) = + (1 )(1 ) = + 1 + b) =
+1 (21) .
= (2 1) + (1 ) .
c) Given an iid random sample, X = number of true responsees Binom(n, ) which has mean n and variance n(1). So = X has mean (and hence is unbiased) and n (1) variance n . If we have a simple random sample, then the expected value remains unchanged because X is still a sum of (dependent) Bernoulli random variables. d) A natural estimator of is = unbiased. e) see above.
n f ) Var() = Var( +1 ) = (21)2 . We can use the identity in (a) to express this in (21) terms of and (but not ) if we like.
(1)
+1 (21) .
E() =
E()+1 (21)
=
+1 (21)
= , so is
g) The variance of each estimator has n in the denominator, ...