Bayesian Computation with R Practice #1-6

Published on by anesth

練習問題1-6
1.
a)
> library(LearnBayes)
> data(studentdata)
> attach(studentdata)

> hist(Dvds)

b)
> summary(Dvds)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s
0.00 10.00 20.00 30.93 30.00 1000.00 16

c)
> table(Dvds)
Dvds
0 1 2 2.5 3 4 5 6 7 8 9 10 11 12 13 14
26 10 13 1 18 9 27 14 12 12 7 78 3 20 7 4
15 16 17 17.5 18 20 21 22 22.5 23 24 25 27.5 28 29 30
46 1 3 1 4 83 3 3 1 3 2 31 3 1 1 45
31 33 35 36 37 40 41 42 45 46 48 50 52 53 55 60
1 1 12 4 1 26 1 1 5 1 2 26 1 2 1 7
62 65 67 70 73 75 80 83 85 90 97 100 120 122 130 137
1 2 1 4 1 3 4 1 1 1 1 10 2 1 2 1
150 152 157 175 200 250 500 900 1000
6 1 1 1 8 1 1 1 1

> barplot(table(Dvds),xlab=”Dvds”,ylab=”Count”)

2.
a)
> boxplot(Height~Gender, ylab=”Height”)

b)
> output <- boxplot(Height~Gender) > output
$stats
[,1] [,2]
[1,] 57.75 65
[2,] 63.00 69
[3,] 64.50 71
[4,] 67.00 72
[5,] 73.00 76

$n
[1] 428 219

$conf
[,1] [,2]
[1,] 64.19451 70.6797
[2,] 64.80549 71.3203

$out
[1] 56 76 55 56 76 54 54 84 78 77 56 63 77 79 62 62 61 79 59 61 78 62

$group
[1] 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2

$names
[1] “female” “male”

c)
> male.Height <- Height[Gender == "male"] > female.Height <- Height[Gender == "female"] > male.Height
[1] 70.000 65.000 72.000 72.000 66.000 71.000 71.000 71.000 68.000 70.000 71.000
[12] 73.000 72.000 67.000 72.000 73.000 71.000 73.000 72.000 68.000 66.000 76.000
[23] 71.000 68.000 73.000 70.000 71.000 74.000 76.000 70.000 72.000 70.000 75.000
[34] 71.000 71.000 72.000 74.000 70.000 74.000 70.000 71.000 72.000 68.000 70.500
[45] 70.000 73.000 68.000 74.000 69.000 63.000 72.000 71.000 75.000 70.000 77.000
[56] 69.000 68.000 69.000 75.000 70.000 69.000 70.000 75.000 71.000 68.000 70.000
[67] 74.000 79.000 NA NA 69.000 69.000 72.000 65.000 62.000 66.000 74.000
[78] 68.000 67.000 71.000 71.500 70.000 70.000 73.000 70.000 68.000 71.000 72.000
[89] 72.000 70.000 72.000 70.000 73.000 72.000 70.000 67.000 70.000 71.000 68.000
[100] 71.000 72.000 70.000 68.000 71.000 62.000 68.000 70.000 72.000 73.000 61.000
[111] 67.000 74.000 66.929 72.000 79.000 66.000 68.000 70.000 69.000 74.000 68.000
[122] 68.000 67.000 72.000 67.000 72.000 72.000 66.000 NA 73.000 71.000 71.000
[133] 72.000 71.000 70.000 70.000 66.000 72.000 70.000 67.000 75.000 74.000 67.000
[144] 72.000 72.000 68.000 73.000 71.000 70.000 71.000 74.000 66.000 72.000 69.000
[155] 67.000 69.000 73.000 70.500 70.000 74.000 68.000 70.000 73.000 59.000 71.000
[166] 72.000 61.000 72.000 70.000 72.000 71.000 71.000 74.000 70.000 70.000 67.000
[177] 75.000 72.000 73.000 71.000 67.000 75.000 69.000 78.000 70.000 70.000 74.000
[188] 72.000 76.000 68.000 73.000 66.000 70.000 72.000 70.000 74.000 75.500 74.000
[199] 62.000 66.000 73.000 70.000 71.000 70.000 71.750 71.000 72.000 73.000 72.500
[210] 66.000 72.000 70.000 72.000 65.000 71.000 71.000 70.000 71.000 74.000 72.000
[221] 68.000 69.000

> female.Height
[1] 67.00 64.00 61.00 61.00 63.00 61.00 64.00 66.00 63.00 67.00 63.50 65.00
[13] 62.00 64.00 61.00 56.00 63.00 65.00 64.00 59.00 61.50 64.00 63.00 66.00
[25] 66.00 63.00 67.00 76.00 61.00 63.00 62.00 65.00 68.00 63.00 62.00 66.00
[37] 62.00 67.00 64.50 66.00 65.00 NA 65.00 65.00 67.00 64.00 67.00 64.00
[49] 66.00 65.00 64.00 68.00 62.00 55.00 69.00 68.00 62.00 62.00 61.00 68.00
[61] 60.00 68.00 56.00 67.00 69.00 72.00 76.00 64.00 72.00 66.00 63.00 63.00
[73] 64.00 69.00 54.00 64.00 62.00 64.00 NA 64.00 66.00 65.00 64.50 65.00
[85] 61.00 69.00 64.00 62.00 68.00 62.00 64.50 67.00 62.00 65.00 66.00 62.00
[97] 67.00 66.00 61.00 63.00 66.00 67.00 62.50 59.75 64.00 67.00 54.00 60.00
[109] 63.00 63.00 64.00 65.00 67.00 62.00 66.00 67.00 64.00 68.00 64.00 61.00
[121] 63.00 64.00 71.00 61.00 61.00 71.00 60.50 58.50 66.50 66.00 64.00 60.00
[133] 65.00 63.00 69.00 66.00 67.00 62.00 67.00 65.00 67.00 64.00 67.00 64.00
[145] 64.00 66.00 64.00 66.00 65.50 63.00 63.00 68.00 67.00 69.00 60.00 62.00
[157] 68.00 69.00 61.00 66.00 70.00 71.00 66.00 67.00 62.00 62.00 64.00 70.00
[169] 65.50 62.50 84.00 60.00 62.00 65.00 66.00 60.00 67.50 58.00 64.00 73.00
[181] 64.00 67.00 63.00 68.00 60.00 62.00 63.00 66.00 68.00 63.00 78.00 62.00
[193] 67.00 63.00 62.00 64.00 68.00 67.50 64.00 64.00 64.00 64.00 70.00 65.00
[205] 62.00 67.00 65.00 60.00 66.00 65.00 65.00 67.00 66.00 69.00 72.00 65.00
[217] 67.00 65.00 65.00 62.00 63.00 65.00 63.00 67.00 69.00 72.00 68.00 60.00
[229] 62.00 65.00 66.00 64.00 67.00 66.00 67.00 62.00 64.00 62.00 63.00 64.00
[241] 64.00 61.00 64.00 62.00 67.00 65.00 64.00 63.00 62.00 72.00 65.00 69.00
[253] 70.00 64.00 66.00 60.50 66.00 64.00 65.00 62.00 64.00 65.00 62.00 62.00
[265] 62.00 63.00 64.00 64.00 NA 64.00 60.00 77.00 58.00 61.00 70.00 65.00
[277] 69.50 72.00 63.00 60.00 60.00 68.50 60.00 59.00 64.00 66.00 68.00 66.00
[289] 61.00 62.00 64.00 66.00 62.00 67.00 61.50 66.00 64.00 67.00 66.00 62.00
[301] 62.00 64.00 67.50 65.00 65.00 69.50 65.00 62.00 63.00 65.00 60.00 64.00
[313] 67.00 65.00 63.00 61.00 60.50 63.00 65.00 68.00 64.00 65.00 62.00 68.00
[325] 64.00 62.00 65.00 62.00 56.00 63.00 68.00 60.00 69.50 69.00 65.00 66.00
[337] 64.00 62.00 65.00 65.00 68.00 66.00 66.50 63.00 65.00 62.00 60.00 64.00
[349] 63.00 64.00 65.00 63.00 70.00 65.00 69.00 63.00 62.00 66.00 71.00 64.00
[361] 62.00 67.00 67.00 64.00 63.00 62.00 68.00 60.00 69.00 70.00 67.00 67.00
[373] 67.50 63.00 62.00 64.00 69.00 64.00 68.00 64.00 68.00 67.00 69.00 63.00
[385] 66.00 61.00 NA 65.00 61.00 63.00 64.00 65.00 65.00 69.00 65.00 63.00
[397] 65.00 66.00 68.00 61.00 NA NA NA 66.00 68.00 65.00 60.00 64.00
[409] 69.00 64.00 66.00 67.00 67.00 64.00 64.00 67.00 61.00 67.00 64.00 67.00
[421] 66.50 61.00 63.00 66.00 57.75 66.00 64.00 68.00 67.00 68.00 65.00 68.00
[433] 71.00 66.00 67.00

> mean(male.Height)
[1] NA
> mean(female.Height)
[1] NA
このようにmean()関数を使うと、NAとなる。データ内にNAが含まれているため、除去しないと平均値が算定できない。
そこで、na.rm=TRUEオプションを使う。

> mean(male.Height, na.rm=TRUE)
[1] 70.50767
> mean(female.Height, na.rm=TRUE)
[1] 64.75701
> diff.Height <- mean(male.Height, na.rm=TRUE) - mean(female.Height, na.rm=TRUE) > diff.Height
[1] 5.750657

3.
a)
> plot(ToSleep, WakeUp)

b)
> plot(ToSleep)
> plot(ToSleep, WakeUp)
> fit <- lm(WakeUp ~ ToSleep) > fit

Call:
lm(formula = WakeUp ~ ToSleep)

Coefficients:
(Intercept) ToSleep
7.9628 0.4247

> abline(fit)

c)
> fit

Call:
lm(formula = WakeUp ~ ToSleep)

Coefficients:
(Intercept) ToSleep
7.9628 0.4247

> fit$coefficients[“ToSleep”]*0+fit$coefficients[“(Intercept)”]
ToSleep
7.962764