［R］母分散σ^2が未知の場合の、女子大生200名の身長の母平均の区間推定（「統計解析のはなし」（東京図書）pp.148-150）

> d1 <- c(159, 158, 151, 167, 151)
> d1 <- c(d1, 160, 160, 158, 160, 158)
> xm <- mean(d1)
> s <- sd(d1)
> ah <- 0.05
> nn <- length(d1)
> t9 <- qt(ah / 2, nn - 1, lower.tail = FALSE)
> print(xm)
[1] 158.2
> print(t9)
[1] 2.262157
> print(s)
[1] 4.613988
> print(nn)
[1] 10
> a <- xm - t9 * s / sqrt(nn)
> b <- xm + t9 * s / sqrt(nn)
> # 母平均μの95%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
154.9≦μ≦161.5
> d2 <- c(156, 160, 153, 153, 166)
> d2 <- c(d2, 157, 158, 169, 165, 159)
> xm <- mean(d2)
> s <- sd(d2)
> ah <- 0.05
> nn <- length(d2)
> t9 <- qt(ah / 2, nn - 1, lower.tail = FALSE)
> a <- xm - t9 * s / sqrt(nn)
> b <- xm + t9 * s / sqrt(nn)
> # 母平均μの95%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
155.7≦μ≦163.5
> xm <- mean(c(d1, d2))
> s <- sd(c(d1, d2))
> ah <- 0.05
> nn <- length(c(d1, d2))
> t19 <- qt(ah / 2, nn - 1, lower.tail = FALSE)
> a <- xm - t19 * s / sqrt(nn)
> b <- xm + t19 * s / sqrt(nn)
> # 母平均μの95%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
156.6≦μ≦161.2

［R］母分散σ^2が既知の場合の、女子大生200名の身長の母平均の区間推定（「統計解析のはなし」（東京図書）pp.146-148）

> d <- c(159, 158, 151, 167, 151)
> d <- c(d, 160, 160, 158, 160, 158)
> s <- 5.04
> xm <- mean(d)
> ah <- 0.05
> z <- qnorm(ah / 2, lower.tail = FALSE)
> nn <- length(d)
> print(xm)
[1] 158.2
> print(z)
[1] 1.959964
> print(s)
[1] 5.04
> print(nn)
[1] 10
> a <- xm - z * s / sqrt(nn)
> b <- xm + z * s / sqrt(nn)
> # 母平均μの95%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
155.1≦μ≦161.3
> ah <- 0.01
> z <- qnorm(ah / 2, lower.tail = FALSE)
> a <- xm - z * s / sqrt(nn)
> b <- xm + z * s / sqrt(nn)
> # 母平均μの99%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
154.1≦μ≦162.3

［R］F分布の各αパーセント点（「統計解析のはなし」（東京図書）p.128）

自由度（4,6）の95パーセント点と99パーセント点を求める。

> qf(0.05, 4, 6, lower.tail = FALSE)
[1] 4.533677
> qf(0.01, 4, 6, lower.tail = FALSE)
[1] 9.148301

［Python］連立一次方程式を解く

solveメソッドを使う。以下の連立一次方程式を解く（解はx=3,y=2）。

 x + 2y = 7
3x - 4y = 1

解いてみる。

>>> import numpy as np
>>> mxaa = np.array([[1, 2], [3, -4]])
>>> mxy = np.array([7, 1])
>>> mxaa
array([[ 1,  2],
       [ 3, -4]])
>>> mxy
array([7, 1])
>>> np.linalg.solve(mxaa, mxy)
array([3., 2.])

同じく以下の連立一次方程式を解く（解はx=2,y=-3,z=-14）。

 3x - 3y + z =  1
 3x + 2y     =  0
-1x - 5y + z = -1

解いてみる。

>>> mxaa = np.array([[3, -3, 1], [3, 2, 0], [-1, -5, 1]])
>>> mxy = np.array([1, 0, -1])
>>> mxaa
array([[ 3, -3,  1],
       [ 3,  2,  0],
       [-1, -5,  1]])
>>> mxy
array([ 1,  0, -1])
>>> np.linalg.solve(mxaa, mxy)
array([  2.,  -3., -14.])

［R］データフレームから特定の列を削除する

[ ]演算子を使用して、削除したい列の列番号を、列番号指定箇所に負数で与える。

> no <- 1:3
> name <- c("A", "B", "C")
> age <- c(20, 30, 40)
> blood <- c("AB", "B", "AB")
> dtf1 <- data.frame(no, name, age, blood)
> print(dtf1)
  no name age blood
1  1    A  20    AB
2  2    B  30     B
3  3    C  40    AB
> dtf2 <- dtf1[, -c(2, 3)]
> print(dtf2)
  no blood
1  1    AB
2  2     B
3  3    AB

［R］t分布の各αパーセント点（「統計解析のはなし」（東京図書）p.125）

> n <- c(5, 10)
> a <- c(0.05, 0.01)
> for (i in 1:length(n)) {
+ for (j in 1:length(a)) {
+ cat(sprintf("両側 a = %.2f t_%d(a/2) = %.3f\n", a[j], n[i], qt(a[j] / 2, n[i], lower.tail = FALSE)))
+ cat(sprintf("片側 a = %.2f t_%d(a)   = %.3f\n", a[j], n[i], qt(a[j], n[i], lower.tail = FALSE)))
+ }
+ }
両側 a = 0.05 t_5(a/2) = 2.571
片側 a = 0.05 t_5(a)   = 2.015
両側 a = 0.01 t_5(a/2) = 4.032
片側 a = 0.01 t_5(a)   = 3.365
両側 a = 0.05 t_10(a/2) = 2.228
片側 a = 0.05 t_10(a)   = 1.812
両側 a = 0.01 t_10(a/2) = 3.169
片側 a = 0.01 t_10(a)   = 2.764

［R］χ^2分布の各αパーセント点（「統計解析のはなし」（東京図書）p.121）

> n <- c(5, 6)
> a <- 0.05
> for (i in 1:length(n)) {
+ cat(sprintf("n = %d, a = %.2f\n", n[i], a))
+ cat("両側かつ下側 ")
+ cat(sprintf("χ_%d^2 = %.6f\n", n[i], qchisq(a / 2, n[i], lower.tail = TRUE)))
+ cat("両側かつ上側 ")
+ cat(sprintf("χ_%d^2 = %.6f\n", n[i], qchisq(a / 2, n[i], lower.tail = FALSE)))
+ cat("片側かつ下側 ")
+ cat(sprintf("χ_%d^2 = %.6f\n", n[i], qchisq(a, n[i], lower.tail = TRUE)))
+ cat("片側かつ上側 ")
+ cat(sprintf("χ_%d^2 = %.6f\n", n[i], qchisq(a, n[i], lower.tail = FALSE)))
+ }
n = 5, a = 0.05
両側かつ下側 χ_5^2 = 0.831212
両側かつ上側 χ_5^2 = 12.832502
片側かつ下側 χ_5^2 = 1.145476
片側かつ上側 χ_5^2 = 11.070498
n = 6, a = 0.05
両側かつ下側 χ_6^2 = 1.237344
両側かつ上側 χ_6^2 = 14.449375
片側かつ下側 χ_6^2 = 1.635383
片側かつ上側 χ_6^2 = 12.591587

［R］WebサイトのWebページをダウンロードする

rvestパッケージのread_html関数を使う。

> library(rvest)
> s <- read_html("https://mainichi.jp")
> print(s)
{html_document}
（以下、表示省略）

プロキシ等の設定で、以下のようなエラーが発生して実行できない場合がある。

> s <- read_html("https://mainichi.jp")
 open.connection(x, "rb") でエラー: 
  Timeout was reached: [○○] Connection timeout after ○○ ms

このような場合は、download.file関数を利用した以下のような関数を自作する。read_html関数による戻り値と全く同じ戻り値が得られる。

read_html2 <- function(urlname) {
    filename <- tempfile()
    download.file(urlname, destfile = filename)
    return(read_html(filename))
}

実行してみる。

> library(rvest)
> read_html2 <- function(urlname) {
+     filename <- tempfile()
+     download.file(urlname, destfile = filename)
+     return(read_html(filename))
+ }
> s <- read_html2("https://mainichi.jp")
 URL 'https://mainichi.jp' を試しています 
downloaded 163 KB
> print(s)
{html_document}
（以下、表示省略）

［Octave］正規分布におけるp値

>> pkg load statistics
>> % 平均10、分散3^2（標準偏差3）の正規分布における確率変数が15のときの下側p値
>> normcdf(15, 10, 3)
ans = 0.9522
>> % 平均10、分散3^2（標準偏差3）の正規分布における確率変数が15のときの上側p値
>> 1 - normcdf(15, 10, 3)
ans = 0.047790
>> % 平均10、分散3^2（標準偏差3）の正規分布における確率変数が15のときの両側p値
>> (1 - normcdf(15, 10, 3)) * 2
ans = 0.095581

［Octave］正規分布におけるパーセント点

>> pkg load statistics
>> % 平均10、分散3^2（標準偏差3）の正規分布における下側5パーセント点
>> norminv(0.05, 10, 3)
ans = 5.0654
>> % 平均10、分散3^2（標準偏差3）の正規分布における上側5パーセント点
>> norminv(1 - 0.05, 10, 3)
ans = 14.935
>> % 平均10、分散3^2（標準偏差3）の正規分布における両側5パーセント点
>> norminv(1 - 0.05 / 2, 10, 3)
ans = 15.880

［Octave］t分布におけるp値

>> % 自由度5のt分布における、確率変数が2のときの上側p値
>> 1 - tcdf(2, 5)
ans = 0.050970
>> # 自由度5のt分布における、確率変数が2のときの両側p値
>> (1 - tcdf(2, 5)) * 2
ans = 0.1019

tcdfは下側p値を返すため、上側p値を得たい場合は1からそれをひく。両側p値は上側p値を2倍する。

［R］正規分布の各αパーセント点（「統計解析のはなし」（東京図書）p.114）

> m <- 0
> s <- 1
> a <- c(0.01, 0.05, 0.1)
> cat("両側パーセント点\n")
両側パーセント点
> for (i in 1:length(a)) {
+ cat(sprintf("a = %.2f, z(a/2) = %.2f\n", a[i], qnorm(p = a[i] / 2, mean = m, sd = s, lower.tail = FALSE)))
+ }
a = 0.01, z(a/2) = 2.58
a = 0.05, z(a/2) = 1.96
a = 0.10, z(a/2) = 1.64
> cat("上側パーセント点\n")
上側パーセント点
> for (i in 1:length(a)) {
+ cat(sprintf("a = %.2f, z(a) = %.2f\n", a[i], qnorm(p = a[i], mean = m, sd = s, lower.tail = FALSE)))
+ }
a = 0.01, z(a) = 2.33
a = 0.05, z(a) = 1.64
a = 0.10, z(a) = 1.28

［Python］ファイルのハッシュ値を得る

hashlibライブラリを使う。以下はPHPの圧縮ファイルのハッシュ値（SHA256）を求めた例。

>>> import hashlib
>>> f = open('php-8.2.4-Win32-vs16-x64.zip', 'br')
>>> b = f.read()
>>> hashlib.sha256(b).hexdigest()
'a3601fe23adfb4985be52eb3a7365716350e4b857c47583673a9aa53162885a3'
>>> f.close()

［R］MLBにおける2001～2010年の得点差分と勝率の散布図（「Rによるセイバーメトリクス入門」（技術評論社）pp.95-99）

> library(tidyverse)
> library(Lahman)
> my_teams <- Teams %>% filter(yearID > 2000 & yearID < 2011) %>%
+ select(teamID, yearID, lgID, G, W, L, R, RA)
> my_teams %>% tail()
    teamID yearID lgID   G  W  L   R  RA
295    SFN   2010   NL 162 92 70 697 583
296    SLN   2010   NL 162 86 76 736 641
297    TBA   2010   AL 162 96 66 802 649
298    TEX   2010   AL 162 90 72 787 687
299    TOR   2010   AL 162 85 77 755 728
300    WAS   2010   NL 162 69 93 655 742
> my_teams <- my_teams %>% mutate(RD = R - RA, Wpct = W / (W + L))
> run_diff <- ggplot(my_teams, aes(x = RD, y = Wpct)) + geom_point() +
+ scale_x_continuous("Run differential") +
+ scale_y_continuous("Winning percentage")
> print(run_diff)
> crcblue <- "#2905A1"
> linfit <- lm(Wpct ~ RD, data = my_teams)
> print(linfit)
Call:
lm(formula = Wpct ~ RD, data = my_teams)
Coefficients:
(Intercept)           RD  
  0.4999909    0.0006216  
> run_diff + geom_smooth(method = "lm", se = FALSE, color = crcblue)
`geom_smooth()` using formula = 'y ~ x'

［Python］一般逆行列を求める

linalg.pinvメソッドを使う。

>>> import numpy as np
>>> mx = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
>>> print(mx)
[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]
 [13 14 15 16]]
>>> np.linalg.det(mx)
-1.820448242817726e-31
>>> np.linalg.inv(mx)
array([[ 1.50119988e+15, -3.75299969e+14, -3.75299969e+15,
         2.62709978e+15],
       [-1.95155984e+16,  1.95155984e+16,  1.95155984e+16,
        -1.95155984e+16],
       [ 3.45275971e+16, -3.79052969e+16, -2.77721977e+16,
         3.11498974e+16],
       [-1.65131986e+16,  1.87649984e+16,  1.20095990e+16,
        -1.42613988e+16]])
>>> np.linalg.pinv(mx)
array([[-0.285 , -0.145 , -0.005 ,  0.135 ],
       [-0.1075, -0.0525,  0.0025,  0.0575],
       [ 0.07  ,  0.04  ,  0.01  , -0.02  ],
       [ 0.2475,  0.1325,  0.0175, -0.0975]])

一般逆行列の定義を満たしているか否か確認。

>>> mxi = np.linalg.pinv(mx)
>>> np.dot(np.dot(mx, mxi), mx)
array([[ 1.,  2.,  3.,  4.],
       [ 5.,  6.,  7.,  8.],
       [ 9., 10., 11., 12.],
       [13., 14., 15., 16.]])

当然、正方行列でなくとも求めることができる。

>>> mx = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
>>> print(mx)
[[1 2 3 4]
 [5 6 7 8]]
>>> np.linalg.det(mx)
（表示省略）
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square
>>> np.linalg.inv(mx)
（表示省略）
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square
>>> np.linalg.pinv(mx)
array([[-5.50000000e-01,  2.50000000e-01],
       [-2.25000000e-01,  1.25000000e-01],
       [ 1.00000000e-01, -1.38777878e-17],
       [ 4.25000000e-01, -1.25000000e-01]])
>>> mxi = np.linalg.pinv(mx)
>>> np.dot(np.dot(mx, mxi), mx)
array([[1., 2., 3., 4.],
       [5., 6., 7., 8.]])

linalg.pinvメソッドはムーア・ペンローズ一般逆行列を数値的に求めるメソッドのため、求まった一般逆行列は近似値であることに注意。ムーア・ペンローズ一般逆行列の定義の4つの式について、比較をした例は以下のとおり。すべての要素がTrueにならないことがわかる。

>>> np.dot(np.dot(mx, mxi), mx) == mx
array([[False, False, False,  True],
       [False, False, False, False]])
>>> np.dot(np.dot(mxi, mx), mxi) == mxi
array([[False, False],
       [False,  True],
       [False, False],
       [False, False]])
>>> np.dot(mx, mxi).T == np.dot(mx, mxi)
array([[ True, False],
       [False,  True]])
>>> np.dot(mxi, mx).T == np.dot(mxi, mx)
array([[ True, False, False, False],
       [False,  True, False,  True],
       [False, False,  True,  True],
       [False,  True,  True,  True]])

［Python］逆行列を求める

linalg.invメソッドを使う。

>>> mx = np.array([[3]])
>>> mx
array([[3]])
>>> np.linalg.inv(mx)
array([[0.33333333]])
>>> mx = np.array([[1, 2], [3, 4]])
>>> mx
array([[1, 2],
       [3, 4]])
>>> np.linalg.inv(mx)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])
>>> mx = np.array([[3, -3, 1], [3, 2, 0], [-1, -5, 1]])
>>> mx
array([[ 3, -3,  1],
       [ 3,  2,  0],
       [-1, -5,  1]])
>>> np.linalg.inv(mx)
array([[ 1. , -1. , -1. ],
       [-1.5,  2. ,  1.5],
       [-6.5,  9. ,  7.5]])

正則ではない行列（行列式の値が0）の場合は、linalg.invメソッドはエラーを返す。

>>> mx = np.array([[1, 4, 7], [2, 5, 8], [3, 6, 9]])
>>> mx
array([[1, 4, 7],
       [2, 5, 8],
       [3, 6, 9]])
>>> np.linalg.det(mx)
0.0
>>> np.linalg.inv(mx)
Traceback (most recent call last):
（表示省略）
numpy.linalg.LinAlgError: Singular matrix

［C#］ディレクトリ（フォルダー）を作成／削除する

作成にはDirectoryクラス（System.IO名前空間）のCreateDirectoryメソッドを、削除にはDeleteメソッドを使う。以下、実行例。DドライブのルートディレクトリにTempディレクトリを一時的に作成して削除している。

> Directory.Exists(@"D:\Temp")
false
> Directory.CreateDirectory(@"D:\Temp")
[Temp]
> Directory.Exists(@"D:\Temp")
true
> Directory.Delete(@"D:\Temp")
> Directory.Exists(@"D:\Temp")
false

［C#］ディレクトリの存在を確認する

Directoryクラス（System.IO名前空間）のExistsメソッドを使用する。第一引数に指定したディレクトリが存在すればtrueを、存在しなければfalseを返す。以下の例のとおりファイルを指定すると、そのファイルが存在したとしてもfalseを返す。

> Directory.Exists(@"C:\Windows")
true
> Directory.Exists(@"C:\Windowss")
false
> Directory.Exists(@"C:\Windows\explorer.exe")
false

［R］データフレームに列を追加する

transform関数を使用する。追加する際には列の名前（等号記号の左側）を指定する必要があることに注意。

> no <- 1:3
> name <- c("A", "B", "C")
> dtf <- data.frame(no, name)
> print(dtf)
  no name
1  1    A
2  2    B
3  3    C
> age <- c(20, 30, 40)
> dtf <- transform(dtf, age = age)
> print(dtf)
  no name age
1  1    A  20
2  2    B  30
3  3    C  40

［R］例外処理を行う

式の評価や関数の実行時に例外（エラー）が発生した際、処理を分岐させる方法。try関数をsilentオプションをTRUEにして式の評価や関数の実行を行えばよい。

以下の例では、「2 * 3」（構文的に正しい）と「2 * "3"」（数値と文字列をかけているため構文的に正しくない）をそれぞれ評価した例。前者であれば、try関数の戻り値は、本来戻される戻り値を返すが、後者の例外が発生すると、try関数は出力を抑制してtry-errorオブジェクトを返す。try-errorオブジェクトか否かは、最後に処理された実行の戻り値を返す.Last.value変数を使って判定をしている。try関数の戻り値で判定しても同じ。

> 2 * 3
[1] 6
> try(2 * 3, silent = TRUE)
[1] 6
> print(.Last.value)
[1] 6
> class(.Last.value)
[1] "numeric"
> .Last.value == "try-error"
[1] FALSE
> try(2 * "3", silent = TRUE)
> print(.Last.value)
[1] "Error in 2 * \"3\" :  二項演算子の引数が数値ではありません \n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in 2 * "3":  二項演算子の引数が数値ではありません >
> class(.Last.value)
[1] "try-error"
> .Last.value == "try-error"
[1] TRUE
> r <- try(2 * "3", silent = TRUE)
> class(r) == "try-error"
[1] TRUE

［R］Webサイトの証明書を調べる

opensslパッケージのdownload_ssl_cert関数を使う。以下は首相官邸のWebサイトの証明書を調べた例。

> library(openssl)
> r <- download_ssl_cert("www.kantei.go.jp")
> print(r)
[[1]]
[x509 certificate] *.kantei.go.jp
md5: d3d1602d86f62630c270db1a641313c1
sha1: 71b0f7c5b22eed8f29e5d29b916363a6e98bbe93
[[2]]
[x509 certificate] GlobalSign RSA OV SSL CA 2018
md5: a75d8b210911ee173efd25e1e00ed6fd
sha1: dfe83023062b997682708b4eab8e819aff5d9775

［R］画像ファイルを読み込む

imagerパッケージのload.image関数を使う。以下の例ではコマンドプロンプトでImageMagickのロゴの画像をあらかじめ作成し、それをカレントディレクトリに置いてある。
コマンドプロンプトで、画像ファイルを作成する

>magick logo: magicklogo.png

Rで実行する。

> library(imager)
> img <- load.image("magicklogo.png")
> print(img)
Image. Width: 640 pix Height: 480 pix Depth: 1 Colour channels: 3 
> plot(img)

［R］ファイルのハッシュ値を得る

opensslパッケージを使う。以下はPHPの圧縮ファイルのハッシュ値（SHA256）を得た例。sha256関数の戻り値はロウ型のため、as.character関数で文字列に変換して表示している。

> library(openssl)
> sha256(file("php-8.2.4-Win32-vs16-x64.zip"))
sha256 a3:60:1f:e2:3a:df:b4:98:5b:e5:2e:b3:a7:36:57:16:35:0e:4b:85:7c:47:58:36:73:a9:aa:53:16:28:85:a3 
> as.character(sha256(file("php-8.2.4-Win32-vs16-x64.zip")))
[1] "a3601fe23adfb4985be52eb3a7365716350e4b857c47583673a9aa53162885a3"

［C#］ファイルの存在を確認する

Fileクラス（System.IO名前空間）のExistsメソッドを使用する。第一引数に指定したファイルが存在すればtrueを、存在しなければfalseを返す。以下の例のとおりディレクトリを指定すると、そのディレクトリが存在したとしてもfalseを返す。

> File.Exists(@"C:\Windows\explorer.exe")
true
> File.Exists(@"C:\Windows\explorerr.exe")
false
> File.Exists(@"C:\Windows")
false
> Directory.Exists(@"C:\Windows")
true

［R］ggplot2で作成した図を画像ファイルに出力する

ggsave関数を使う。使用すると、最後に表示させた図を画像ファイルに保存する。ファイル名で拡張子を指定すれば、それを自動で判別してその保存形式で保存する。widthオプションなどを指定することで画像のサイズを指定することができる。

> library(ggplot2)
> dtf <- data.frame(xp = 0:5, yp = (0:5) ^ 2)
> ggplot(dtf, aes(x = xp, y = yp)) + geom_point()
> ggsave("fig_64.jpg", width = 6, height = 4, units = "cm")
> ggsave("fig_44.png", width = 4, height = 4, units = "cm")

plotオプションを使うことで、過去に作成した任意の図を画像ファイルに出力することができる。

> g1 <- ggplot(dtf, aes(x = xp * 4, y = yp)) + geom_point()
> g2 <- ggplot(dtf, aes(x = xp * 40, y = yp)) + geom_point()
> ggsave("fig_46_1.png", plot= g1, width = 4, height = 6, units = "cm")
> ggsave("fig_46_2.png", plot= g1, width = 4, height = 6, units = "cm")

［R］パッケージをインストールするサーバーを指定する

options関数にreposオプションを付けて使う。reposオプションには、このページ（https://cran.r-project.org/mirrors.html）に掲載されているCRANのミラーサイトのURLを指定する。

日本国内のミラーサイトを指定して、ggplot2パッケージをインストールしてみる。

> options(repos = "https://ftp.yz.yamagata-u.ac.jp/pub/cran/")
> install.packages("ggplot2")

このoptions関数でCRANのミラーサイトを指定すると、install.packages関数でインストールする際は、ミラーサイトを選択させるダイアログは表示されずに、すぐにこのミラーサイトから読み込みを始める。