R(本の計算を再現)

2023年1月28日 (土)

[R]分散と標準偏差(「まずはこの一冊から 意味が分かる統計学」(ペレ出版)、pp.27-28)

p.27の観測値をファイルtable_p027.csvに保存しておく。

> dtf <- read.csv("table_p027.csv")
> x <- dtf$x
> n <- length(x)
> m <- mean(x)
> v <- sum((x - m) ^ 2) / n
> s <- sqrt(v)
> cat(sprintf("分散 = %f\n", v))
分散 = 163.700000
> cat(sprintf("標準偏差 = %f\n", s))
標準偏差 = 12.794530

2023年1月27日 (金)

[R]共分散(「まずはこの一冊から 意味が分かる統計解析」(ペレ出版)、p.152)

以下のp.152の表の数値をファイルtable_p214.csvに保存して、カレントディレクトリに置いておく。

 i,     h,    w
1, 147.9, 41.7
2, 163.5, 60.2
3, 159.8, 47.0
4, 155.1, 53.2
5, 163.3, 48.3
6, 158.7, 55.2
7, 172.0, 58.5
8, 161.2, 49.0
9, 153.9, 46.7
10, 161.6, 52.5

計算する。

> dtf <- read.csv("table_p214.csv")
> n <- nrow(dtf)
> h <- dtf$h
> w <- dtf$w
> mh <- mean(h)
> mw <- mean(w)
> sxy <- sum((h - mh) * (w - mw)) / n
> spxy <- sum((h - mh) * (w - mw)) / (n - 1)
> cat(sprintf("身長の平均 = %f\n", mh))
身長の平均 = 159.700000
> cat(sprintf("体重の平均 = %f\n", mw))
体重の平均 = 51.230000
> cat(sprintf("sxy = %f\n", sxy))
sxy = 23.730000
> cat(sprintf("cov関数の結果 = %f\n", cov(h, w)))
cov関数の結果 = 26.366667
> cat(sprintf("s'xy = %f\n", spxy))
s'xy = 26.366667

2023年1月26日 (木)

[R]偏差値(「まずはこの一冊から 意味が分かる統計解析」(ペレ出版)、p.47)

以下のp.47の表の数値をファイルtable_p047.csvに保存して、カレントディレクトリに置いておく。

i,  x
1, 61
2, 59
3, 60
4, 67
5, 53

計算する。

> dtf <- read.csv("table_p047.csv")
> n <- nrow(dtf)
> x <- dtf$x
> mx <- mean(x)
> s2 <- sum((x - mx) ^ 2) / n
> s <- sqrt(s2)
> z <- 50 + 10 * (x - mx) / s
> cat(sprintf("平均 %f\n", mx))
平均 60.000000
> cat(sprintf("標準偏差 %f\n", s))
標準偏差 4.472136
> cat("番号 点数 偏差値\n")
番号 点数 偏差値
> for (i in 1:n) {
+ cat(sprintf("%d %d %.1f\n", i, x[i], z[i]))
+ }
1 61 52.2
2 59 47.8
3 60 50.0
4 67 65.7
5 53 34.3

2023年1月22日 (日)

[R]平均(「まずはこの一冊から 意味が分かる統計学」(ペレ出版)、p.27)

以下の観測値をファイルtable_p027.csvに保存しておく。

 i,  x
1, 43
2, 47
3, 52
4, 52
5, 54
6, 61
7, 67
8, 67
9, 68
10, 69
11, 70
12, 71
13, 71
14, 73
15, 76
16, 78
17, 82
18, 84
19, 84
20, 91

計算する。

> dtf <- read.csv("table_p027.csv")
> x <- dtf$x
> m <- mean(x)
> cat(sprintf("%f\n", m))
68.000000

2023年1月 7日 (土)

[R]コレスキー分解を行う(「最小二乗法による実験データ解析」、p.81)

以下は、式(4.115)に基に作成したコレスキー分解を行う関数。行列 B を上三角行列 L の転置行列と上三角行列の積に分解している。細かなチェックは行っていないので注意。

B = R^T R
choldecomp <- function(mxbb) {
m <- ncol(mxbb)
mxrr <- matrix(0., m, m)
mxrr[1, 1] <- sqrt(mxbb[1, 1])
for (j in 2:m) {
mxrr[1, j] <- mxbb[1, j] / mxrr[1, 1]
}
for (i in 2:m) {
k <- 1:(i - 1)
mxrr[i, i] <- sqrt(mxbb[i, i] - sum(mxrr[1:(i - 1), i] ^ 2))
for (j in i:m) {
mxrr[i, j] <- (mxbb[i, j]
- sum(mxrr[k, i] * mxrr[k, j])) / mxrr[i, i]
}
}
return(mxrr)
}

以下、動作確認。

> d <- c(1, 1, 1, 1, 1, 2, 3, 4, 1, 3, 6, 10, 1, 4, 10, 20)
> mxbb <- matrix(d, 4, 4, byrow = TRUE)
> mxrr <- choldecomp(mxbb)
> cat("B\n")
B
> print(mxbb)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 2 3 4
[3,] 1 3 6 10
[4,] 1 4 10 20
> cat("R^T\n")
R^T
> print(t(mxrr))
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 1 1 0 0
[3,] 1 2 1 0
[4,] 1 3 3 1
> cat("R\n")
R
> print(mxrr)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 0 1 2 3
[3,] 0 0 1 3
[4,] 0 0 0 1
> cat("R^T T\n")
R^T T
> print(t(mxrr) %*% mxrr)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 2 3 4
[3,] 1 3 6 10
[4,] 1 4 10 20

Rへの組み込み関数cholでも計算してみる。上の結果(mxrr, R)と一致していることがわかる。

> base::chol(mxbb)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 0 1 2 3
[3,] 0 0 1 3
[4,] 0 0 0 1

2023年1月 6日 (金)

[R]例3.3 配達時間 (3)(「線形回帰分析」、pp.113-114)

簡略化のため、lm関数を使用して計算している。

> dtf <- read.csv("table3_1.csv", header = TRUE)
> xx2i <- dtf$CASE
> xx3i <- dtf$DIS ^ 2 / 1.e5
> yyi <- dtf$DVT
> # (3.64)式
> r <- lm(yyi[-c(11)] ~ xx2i[-c(11)] + xx3i[-c(11)])
> print(summary(r))
Call:
lm(formula = yyi[-c(11)] ~ xx2i[-c(11)] + xx3i[-c(11)])
Residuals:
Min 1Q Median 3Q Max
-3.8002 -1.2991 0.3215 1.3196 3.0112
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.4226 0.7920 8.110 6.62e-08 ***
xx2i[-c(11)] 1.3376 0.1196 11.182 2.65e-10 ***
xx3i[-c(11)] 1.4969 0.1831 8.178 5.78e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.125 on 21 degrees of freedom
Multiple R-squared: 0.9826, Adjusted R-squared: 0.9809
F-statistic: 592.7 on 2 and 21 DF, p-value: < 2.2e-16
> # (3.66)式
> ddi <- rep(0., nrow(dtf))
> ddi[11] <- 1
> r <- lm(yyi ~ xx2i + xx3i + ddi)
> print(summary(r))
Call:
lm(formula = yyi ~ xx2i + xx3i + ddi)
Residuals:
Min 1Q Median 3Q Max
-3.800 -1.279 0.303 1.298 3.011
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.4226 0.7920 8.110 6.62e-08 ***
xx2i 1.3376 0.1196 11.182 2.65e-10 ***
xx3i 1.4969 0.1831 8.178 5.78e-08 ***
ddi 5.4200 2.2536 2.405 0.0255 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.125 on 21 degrees of freedom
Multiple R-squared: 0.9836, Adjusted R-squared: 0.9813
F-statistic: 419.9 on 3 and 21 DF, p-value: < 2.2e-16

2023年1月 5日 (木)

[R]例3.3 配達時間 (1),(2)(「線形回帰分析」、pp.106-113)※lm関数使用版

> dtf <- read.csv("data/table3_1.csv", header = TRUE)
> n <- nrow(dtf)
> xx2i <- as.double(dtf$CASE)
> xx3i <- as.double(dtf$DIS) ^ 2 / 1.e5
> yyi <- as.double(dtf$DVT)
> #
> r <- lm(yyi ~ xx2i + xx3i)
> sr <- summary(r)
> # 求めた回帰係数β,回帰係数の標準誤差s,仮説を検定するためのt値,p値
> print(sr$coef)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.191398 0.8673945 7.137927 3.711853e-07
xx2i 1.410930 0.1276345 11.054459 1.890317e-10
xx3i 1.424706 0.1992415 7.150648 3.609898e-07
> # R^2
> print(sr$r.squared)
[1] 0.9790869
> # R~^2
> print(sr$adj.r.squared)
[1] 0.9771857
> # log L
> print(logLik(r))
'log Lik.' -55.18207 (df=4)
> # AIC
> print(AIC(r))
[1] 118.3641
> # SBIC
> print(BIC(r))
[1] 123.2396
 

2023年1月 4日 (水)

[R]例3.3 配達時間 (1),(2)(「線形回帰分析」、pp.106-113)

式(3.62)のモデルで計算している。表3.2(p.112)の対数尤度などは完全に一致しない。書籍では、途中で値を切り捨てて計算をしているからだと思われる。

> dtf <- read.csv("table3_1.csv", header = TRUE)
> n <- nrow(dtf)
> xx2i <- as.double(dtf$CASE)
> xx3i <- as.double(dtf$DIS) ^ 2 / 1.e5
> yyi <- as.double(dtf$DVT)
> k <- 3
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1., n), xx2i, xx3i), ncol = 3)
> mxbh <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> mxyh <- mxxx %*% mxbh
> mxe <- mxy - mxyh
> ssy2 <- sum(t(mxy) %*% mxy) - sum(yyi) ^ 2 / n
> ssyh2 <- sum(t(mxyh) %*% mxyh) - sum(yyi) ^ 2 / n
> sse2 <- as.vector(t(mxe) %*% mxe)
> s2 <- sse2 / (n - k)
> s <- sqrt(s2)
> mxvv <- s2 * solve(t(mxxx) %*% mxxx)
> mxs <- sqrt(diag(mxvv))
> rr2 <- ssyh2 / ssy2
> rrb2 <- 1 - (n - 1) / (n - k) * (1 - rr2)
> tbh <- as.vector(mxbh / mxs)
> sh2 <- sse2 / n
> logll <- -n / 2 * log(2 * pi) - n / 2 * log(sh2) - sse2 / (2 * sh2)
> p <- k + 1
> aic <- -2 * logll + 2 * p
> sbic <- -2 * logll + p * log(n)
> gcv <- -2 * logll - 2 * n * log(1 - p / n)
> hq <- -2 * logll + 2 * p * log(log(n))
> yyhi <- as.vector(mxyh)
> ei <- yyi - yyhi
> ri <- ei / yyi * 100
> ai2 <- ei ^ 2 / sum(ei ^ 2) * 100
> #
> cat("求めた回帰係数 β\n")
求めた回帰係数 β
> print(mxbh)
[,1]
[1,] 6.191398
[2,] 1.410930
[3,] 1.424706
> cat("回帰係数の標準誤差 s\n")
回帰係数の標準誤差 s
> print(mxs)
[1] 0.8673945 0.1276345 0.1992415
> cat("仮説を検定するためのt値\n")
仮説を検定するためのt値
> print(tbh)
[1] 7.137927 11.054459 7.150648
> cat(sprintf("R^2 = %8.4f\n", rr2))
R^2 = 0.9791
> cat(sprintf("R~^2 = %8.4f\n", rrb2))
R~^2 = 0.9772
> cat(sprintf("log L = %8.4f\n", logll))
log L = -55.1821
> cat(sprintf("AIC = %8.4f\n", aic))
AIC = 118.3641
> cat(sprintf("SBIC = %8.4f\n", sbic))
SBIC = 123.2396
> cat(sprintf("GCV = %8.4f\n", gcv))
GCV = 119.0818
> cat(sprintf("HQ = %8.4f\n", hq))
HQ = 119.7164
> st <- sprintf("%5.2f %5.2f %5.2f %6.2f %5.2f", yyi, yyhi, ei, ri, ai2)
> for (i in 1:n) cat(sprintf("%2d %s\n", i, st[i]))
1 16.68 20.54 -3.86 -23.12 12.29
2 11.50 11.11 0.39 3.36 0.12
3 12.03 12.07 -0.04 -0.34 0.00
4 14.88 11.93 2.95 19.85 7.21
5 13.75 14.98 -1.23 -8.93 1.25
6 18.11 17.62 0.49 2.71 0.20
7 8.00 9.19 -1.19 -14.82 1.16
8 17.83 16.70 1.13 6.36 1.06
9 79.24 78.89 0.35 0.44 0.10
10 21.50 18.46 3.04 14.14 7.64
11 40.33 35.51 4.82 11.95 19.20
12 21.00 20.96 0.04 0.19 0.00
13 13.50 12.76 0.74 5.47 0.45
14 19.75 17.70 2.05 10.39 3.48
15 24.00 21.75 2.25 9.38 4.19
16 29.00 28.88 0.12 0.41 0.01
17 15.35 15.23 0.12 0.80 0.01
18 19.00 16.32 2.68 14.13 5.95
19 9.50 10.44 -0.94 -9.92 0.73
20 35.10 38.62 -3.52 -10.04 10.27
21 17.90 20.58 -2.68 -14.97 5.94
22 52.32 52.22 0.10 0.19 0.01
23 18.75 21.77 -3.02 -16.13 7.56
24 19.83 23.22 -3.39 -17.11 9.52
25 10.75 12.16 -1.41 -13.08 1.63

2023年1月 3日 (火)

[R]推定した回帰モデルのパラメーターの分散と標準誤差(「44の例題で学ぶ計量経済学」、p.171)

> dtf <- read.csv("table3_2.csv", header = TRUE)
> n <- nrow(dtf)
> kk <- 2
> xxm <- mean(dtf$xxi)
> yym <- mean(dtf$yyi)
> ssxy <- sum((dtf$xxi - xxm) * (dtf$yyi - yym))
> ssxx <- sum((dtf$xxi - xxm) ^ 2)
> bh <- ssxy / ssxx
> ah <- yym - bh * xxm
> yyhi <- ah + bh * dtf$xxi
> ui <- dtf$yyi - yyhi
> sh2 <- sum(ui ^ 2) / (n - kk)
> sbh2 <- sh2 / ssxx
> sbh <- sqrt(sbh2)
> sah2 <- sh2 * sum(dtf$xxi ^ 2) / (n * ssxx)
> sah <- sqrt(sah2)
> # β^の分散
> print(sbh2)
[1] 0.045
> # β^の標準誤差
> print(sbh)
[1] 0.212132
> # α^の分散
> print(sah2)
[1] 1.35
> # α^の標準誤差
> print(sah)
[1] 1.161895

2022年12月 5日 (月)

[R]女子大生200人の体重とバストの歪度と尖度(「統計解析のはなし」(東京図書)、p.43)

> dtf <- read.csv("table1_1_1.csv", header = TRUE)
> rdtf <- data.frame(c(0., 0.), c(0., 0.))
> colnames(rdtf) <- c("体重", "バスト")
> rownames(rdtf) <- c("歪度 a3", "尖度 a4")
> for (i in 3:4) {
+ xi <- dtf[, i]
+ nn <- length(xi)
+ xm <- mean(xi)
+ ss <- sqrt((nn - 1) / nn * var(xi))
+ a3 <- sum((xi - xm) ^ 3) / (nn * ss ^ 3)
+ a4 <- sum((xi - xm) ^ 4) / (nn * ss ^ 4) - 3
+ rdtf[1, i - 2] <- a3
+ rdtf[2, i - 2] <- a4
+ }
> print(rdtf)
体重 バスト
歪度 a3 0.6519655 1.190016
尖度 a4 0.8305713 3.170911

より以前の記事一覧

無料ブログはココログ