［R］複数のパラメーターと定数を持つ関数の値が最大・最小となるパラメーターを求める

optim関数を使う。第一引数にパラメーターの初期値を、第二引数に関数を指定する。関数に与えたい定数は、関数の引数に作成しておく。

optim(パラメーターの初期値, 関数, …)

以下は、自作関数fの値が最大・最小となるパラメーターを求めた例。この例では、パラメーター（未知数）は2個（ベクトルd）、関数を決める定数はc1とc2の2個。optim関数はデフォルトでは関数の値が最小となるパラメーターを推定する。推定するパラメーター（$par）は数値的に解いて求めているため、必ずしも切りのよい値にはならないことに注意。

> f <- function(d, c1, c2) {(d[1] - (c1 + c2)) ^ 2 + (d[2] - c1) ^ 4 - 2}
> optim(c(0, 0), f, method = "BFGS", c1 = 1, c2 = 2)
$par
[1] 2.999486 0.966386
$value
[1] -1.999998
$counts
function gradient 
      19       16 
$convergence
[1] 0
$message
NULL
> f(c(3, 1), 1, 2)
[1] -2
> optim(c(0, 0), f, method = "BFGS", c1 = 3, c2 = 1)
$par
[1] 3.999491 2.969365
$value
[1] -1.999999
$counts
function gradient 
      24       21 
$convergence
[1] 0
$message
NULL
> f(c(4, 3), 3, 1)
[1] -2

関数の値が最大となるパラメーターを求めたい場合は、controlオプションに「list(fnscale = -1)」と指定すればよい。

> f <- function(d, c1, c2) {-((d[1] - (c1 + c2)) ^ 2 + (d[2] - c1) ^ 4 - 2)}
> optim(c(0, 0), f, control = list(fnscale = -1), method = "BFGS", c1 = 1, c2 = 1)
$par
[1] 2.0003449 0.9741651
$value
[1] 1.999999
$counts
function gradient 
      24       21 
$convergence
[1] 0
$message
NULL

［R］数値を二進法で表記する

以下は、数値（1,3,5,7,11,13）を二進法で表記した例。「0」と「1」からなる文字列で表記している。intToBits関数は、与えた数値を二進数に変換して返すが、戻り値はraw型でかつ長さが32のベクトルのため、最初の4ビット分だけを表示するように工夫している。

> paste0(rev(as.character(as.integer(intToBits(1)))[1:4]), collapse = "")
[1] "0001"
> paste0(rev(as.character(as.integer(intToBits(3)))[1:4]), collapse = "")
[1] "0011"
> paste0(rev(as.character(as.integer(intToBits(5)))[1:4]), collapse = "")
[1] "0101"
> paste0(rev(as.character(as.integer(intToBits(7)))[1:4]), collapse = "")
[1] "0111"
> paste0(rev(as.character(as.integer(intToBits(11)))[1:4]), collapse = "")
[1] "1011"
> paste0(rev(as.character(as.integer(intToBits(13)))[1:4]), collapse = "")
[1] "1101"

［R］一行で複数のコマンドを実行する

;（セミコロン）を使う。セミコロンは、コマンドを区切るための記号として機能する。

> n <- 0; n <- n + 2; n <- n ^ 16; print(n)
[1] 65536
> s
エラー: オブジェクト 's' がありません
> s <- "セナディア役の鈴木みのりさん、かわいい"; print(s)
[1] "セナディア役の鈴木みのりさん、かわいい

上記は対話型実行環境で手順を減らすための手法であり、これのスクリプト内での使用は推奨しない。

一般に、スクリプトで一行に複数のコマンドを書くことはスクリプトの可読性を著しく損なうことから、スクリプトでは一行一コマンドを徹底すること。

［R］コマンドの実行に要した時間を計測する

system.timeの引数にコマンドを与えて実行すると、そのコマンドを実行しつつ実行に要した時間を計測して返す。以下は、それぞれ100次、1000次、2000次の正方行列の逆行列を求めた際に要した時間を計測した例。「経過」に表示されている時間が、実際に要した時間（単位：秒）。この関数の戻り値はproc_time型。

> n <- 100; mx <- matrix(rnorm(n * n), n, n)
> dim(mx)
[1] 100 100
> system.time(mxi <- solve(mx))
  ユーザ システム     経過 
       0        0        0 
> n <- 1000; mx <- matrix(rnorm(n * n), n, n)
> dim(mx)
[1] 1000 1000
> system.time(mxi <- solve(mx))
  ユーザ システム     経過 
    0.58     0.00     0.59 
> n <- 2000; mx <- matrix(rnorm(n * n), n, n)
> dim(mx)
[1] 2000 2000
> system.time(mxi <- solve(mx))
  ユーザ システム     経過 
    5.45     0.01     5.50 
> r <- system.time(mxi <- solve(mx))
> length(r)
[1] 5
> class(r)
[1] "proc_time"
> str(r)
 'proc_time' Named num [1:5] 5.79 0.03 5.84 NA NA
 - attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ...
> print(r)
  ユーザ システム     経過 
    5.79     0.03     5.84 

［R］複数の各ベクトルの要素同士によるすべての組合せを作成する

tidyrパッケージ（tidyverseパッケージに含まれている）のexpand_gridを使う。戻り値はtibble。以下では簡潔な表示にするためデータフレームに変換して表示している。

> s1 <- c("ルアン・メェイ", "帰忘の流離人", "キャストリス")
> s2 <- c("セナディア", "ヘリア")
> s3 <- c("カンタレラ", "シャコンヌ")
> library(tidyr)
> expand_grid(h3 = s1, hs = s2) |> as.data.frame()
              h3         hs
1 ルアン・メェイ セナディア
2 ルアン・メェイ     ヘリア
3   帰忘の流離人 セナディア
4   帰忘の流離人     ヘリア
5   キャストリス セナディア
6   キャストリス     ヘリア
> expand_grid(h3 = s1, hs = s2, ww = s3) |> as.data.frame()
               h3         hs         ww
1  ルアン・メェイ セナディア カンタレラ
2  ルアン・メェイ セナディア シャコンヌ
3  ルアン・メェイ     ヘリア カンタレラ
4  ルアン・メェイ     ヘリア シャコンヌ
5    帰忘の流離人 セナディア カンタレラ
6    帰忘の流離人 セナディア シャコンヌ
7    帰忘の流離人     ヘリア カンタレラ
8    帰忘の流離人     ヘリア シャコンヌ
9    キャストリス セナディア カンタレラ
10   キャストリス セナディア シャコンヌ
11   キャストリス     ヘリア カンタレラ
12   キャストリス     ヘリア シャコンヌ
> expand_grid(h3 = s1, hs = s2, ww = s3) |> class()
[1] "tbl_df"     "tbl"        "data.frame"

［R］複数の各ベクトルの要素同士によるすべての組合せを作成する

expand.gridを使う。戻り値はデータフレーム。

> s1 <- c("ルアン・メェイ", "帰忘の流離人", "キャストリス")
> s2 <- c("セナディア", "ヘリア")
> s3 <- c("カンタレラ", "シャコンヌ")
> expand.grid(h3 = s1, hr = s2)
              h3         hr
1 ルアン・メェイ セナディア
2   帰忘の流離人 セナディア
3   キャストリス セナディア
4 ルアン・メェイ     ヘリア
5   帰忘の流離人     ヘリア
6   キャストリス     ヘリア
> expand.grid(h3 = s1, hs = s2)
              h3         hs
1 ルアン・メェイ セナディア
2   帰忘の流離人 セナディア
3   キャストリス セナディア
4 ルアン・メェイ     ヘリア
5   帰忘の流離人     ヘリア
6   キャストリス     ヘリア
> expand.grid(h3 = s1, hs = s2, ww = s3)
               h3         hs         ww
1  ルアン・メェイ セナディア カンタレラ
2    帰忘の流離人 セナディア カンタレラ
3    キャストリス セナディア カンタレラ
4  ルアン・メェイ     ヘリア カンタレラ
5    帰忘の流離人     ヘリア カンタレラ
6    キャストリス     ヘリア カンタレラ
7  ルアン・メェイ セナディア シャコンヌ
8    帰忘の流離人 セナディア シャコンヌ
9    キャストリス セナディア シャコンヌ
10 ルアン・メェイ     ヘリア シャコンヌ
11   帰忘の流離人     ヘリア シャコンヌ
12   キャストリス     ヘリア シャコンヌ
> class(expand.grid(h3 = s1, hs = s2, ww = s3))
[1] "data.frame"

［Octave］重み付き非線形最小二乗法を行う

optimパッケージのnlinfit関数を使う。

以下の例では、与えられた8組の測定値を使用して、非線形最小二乗法により回帰モデルy = b1 + b2 * x * exp(b3 * x)（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。

>> % 計算の準備
>> pkg load optim
>> x = [1 2 3 5 6 7 8 9]';
>> y = [1 4 10 25 32 49 64 75]';
>> w = repmat(1, 1, length(x))';
>> w(3) = 100;
>> w(5) = 100;
>> f = @(b, x) (b(1) + b(2) * x .* exp(b(3) * x));
>> b0 = [1 1 1];
>>
>> % 重みを考慮しない計算(図の青実線)
>> [b, res, jj, covm, rse] = nlinfit(x, y, f, b0);
>> yest1 = f(b, x);
>> disp(b)     % 推定した回帰係数
  -3.8908
   3.3950
   0.1085
>> disp(rse)   % 回帰値の標準誤差
6.2601
>> disp(covm)  % 推定した回帰係数の分散共分散行列
   6.0731e+00  -1.5657e+00   4.4934e-02
  -1.5657e+00   5.5201e-01  -1.7014e-02
   4.4934e-02  -1.7014e-02   5.3738e-04
>>
>> % 重みを考慮した計算(図の赤実線)
>> [b, res, jj, covm, rse] = nlinfit(x, y, f, b0, [], 'weights', w);
>> yest2 = f(b, x);
>> disp(b)     % 推定した回帰係数
  -0.030452
   2.060791
   0.159292
>> disp(rse)   % 回帰値の標準誤差
13.991
>> disp(covm)  % 推定した回帰係数の分散共分散行列
   1.3955e+00  -3.7178e-01   2.1961e-02
  -3.7178e-01   1.1270e-01  -6.9215e-03
   2.1961e-02  -6.9215e-03   4.3259e-04
>>
>> % 図の描画
>> plot(x, yest1, 'b', x, yest2, 'r', x, y, 'ko', 'MarkerFaceColor', 'k')

［R］重み付き線形最小二乗法を行う

lm関数を使う。weightsオプションに重みを指定すると、その重みを考慮した計算を行う。

以下の例では、与えられた8組の測定値を使用して、線形最小二乗法により回帰モデルy = b1 + b2 * x + b3 * x ^ 2（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。

> # 計算の準備
> x <- c(1, 2, 3, 5, 6, 7, 8, 9)
> y <- c(1, 4, 10, 25, 32, 49, 64, 75)  # 測定値（図の黒丸）
> w <- rep(1, length(x))
> w[c(3, 5)] <- 100
> print(data.frame(x, y, w))
  x  y   w
1 1  1   1
2 2  4   1
3 3 10 100
4 5 25   1
5 6 32 100
6 7 49   1
7 8 64   1
8 9 75   1
> 
> # 重みを考慮しない計算（図の青実線）
> r1 <- lm(y ~ 1 + x + I(x ^ 2))
> print(r1)
Call:
lm(formula = y ~ 1 + x + I(x^2))
Coefficients:
(Intercept)            x       I(x^2)  
    -0.6731       0.7042       0.8792  
> sm <- summary(r1)
> print(sm)
Call:
lm(formula = y ~ 1 + x + I(x^2))
Residuals:
       1        2        3        4        5        6        7        8 
 0.08963 -0.25225  0.64741  0.17136 -3.20434  1.66150  2.76888 -1.88219 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.6731     2.8545  -0.236  0.82295   
x             0.7042     1.3677   0.515  0.62859   
I(x^2)        0.8792     0.1351   6.507  0.00128 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.225 on 5 degrees of freedom
Multiple R-squared:  0.9954,    Adjusted R-squared:  0.9936 
F-statistic: 546.6 on 2 and 5 DF,  p-value: 1.399e-06
> print(r1)                              # lm関数による計算結果の概要
Call:
lm(formula = y ~ 1 + x + I(x^2))
Coefficients:
(Intercept)            x       I(x^2)  
    -0.6731       0.7042       0.8792  
> print(sm$coef)                         # 推定した回帰係数,その標準誤差,t値,p値
              Estimate Std. Error    t value    Pr(>|t|)
(Intercept) -0.6730552  2.8545448 -0.2357837 0.822953999
x            0.7041995  1.3677286  0.5148679 0.628589629
I(x^2)       0.8792277  0.1351261  6.5067205 0.001280624
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 24.75789
> print(sum((y - fitted(r1)) ^ 2))
[1] 24.75789
> print(sm$df)                           # モデルのランク,自由度,回帰係数の数
[1] 3 5 3
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 4.951578
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 2.225214
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            (Intercept)          x      I(x^2)
(Intercept)   8.1484260 -3.5141134  0.31168334
x            -3.5141134  1.8706816 -0.18061352
I(x^2)        0.3116833 -0.1806135  0.01825906
> 
> # 重みを考慮した計算（図の赤実線）
> r2 <- lm(y ~ 1 + x + I(x ^ 2), weights = w)
> print(r2)
Call:
lm(formula = y ~ 1 + x + I(x^2), weights = w)
Coefficients:
(Intercept)            x       I(x^2)  
      8.024       -2.828        1.143  
> sm <- summary(r2)
> print(sm)
Call:
lm(formula = y ~ 1 + x + I(x^2), weights = w)
Weighted Residuals:
     1      2      3      4      5      6      7      8 
-5.340 -2.942  1.701  2.535 -2.128  4.754  5.433 -0.173 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   8.0244     3.7465   2.142  0.08512 . 
x            -2.8277     1.8037  -1.568  0.17773   
I(x^2)        1.1432     0.1958   5.839  0.00208 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.542 on 5 degrees of freedom
Multiple R-squared:  0.9966,    Adjusted R-squared:  0.9953 
F-statistic: 734.5 on 2 and 5 DF,  p-value: 6.701e-07
> print(r2)                              # lm関数による計算結果の概要
Call:
lm(formula = y ~ 1 + x + I(x^2), weights = w)
Coefficients:
(Intercept)            x       I(x^2)  
      8.024       -2.828        1.143  
> print(sm$coef)                         # 推定した回帰係数,その標準誤差,t値,p値
             Estimate Std. Error   t value    Pr(>|t|)
(Intercept)  8.024356  3.7464592  2.141851 0.085117824
x           -2.827713  1.8037101 -1.567721 0.177732033
I(x^2)       1.143186  0.1957718  5.839381 0.002083779
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 103.1604
> print(sum(w * (y - fitted(r2)) ^ 2))
[1] 103.1604
> print(sm$df)                           # モデルのランク,自由度,回帰係数の数
[1] 3 5 3
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 20.63208
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 4.542255
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            (Intercept)          x      I(x^2)
(Intercept)  14.0359563 -6.6722356  0.70866822
x            -6.6722356  3.2533702 -0.35090188
I(x^2)        0.7086682 -0.3509019  0.03832661
> 
> # 図の描画
> plot(x, fitted(r1), type = "l", col = "blue", xlab = "x", ylab = "y")
> lines(x, fitted(r2), col = "red")
> points(x, y, pch = 20)

［R］組合せを得る

n個からr個とった組合せの数を得るにはchoose関数を使うが、その組合せ自体を得るにはcombn関数を使う。以下は5個の中から3個を選んだ例。₅C₃=10で異なる10通りの選び方がある。

> choose(5, 3)
[1] 10
> combn(1:5, 3)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    1    1    1    1    1    2    2    2     3
[2,]    2    2    2    3    3    4    3    3    4     4
[3,]    3    4    5    4    5    5    4    5    5     5
> combn(c("a", "b", "c", "d", "e"), 3)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] "a"  "a"  "a"  "a"  "a"  "a"  "b"  "b"  "b"  "c"  
[2,] "b"  "b"  "b"  "c"  "c"  "d"  "c"  "c"  "d"  "d"  
[3,] "c"  "d"  "e"  "d"  "e"  "e"  "d"  "e"  "e"  "e"  

［R］2組のベクトルの各要素同士によるすべての組合せによる演算結果を得る

outer関数もしくは%o%演算子を使う。outer関数は特に演算子を指定しなければ積を計算する。%o%演算子も積を計算する。これを和や差にしたい場合はouter関数の第三引数に、その演算子を指定する。戻り値は行列。行と列を入れ替えたければt関数を使えばよい。第三引数には関数を与えてもよい。

> x <- c(1, 2, 3)
> y <- c(10, 20, 30, 40)
> outer(x, y)
     [,1] [,2] [,3] [,4]
[1,]   10   20   30   40
[2,]   20   40   60   80
[3,]   30   60   90  120
> x %o% y
     [,1] [,2] [,3] [,4]
[1,]   10   20   30   40
[2,]   20   40   60   80
[3,]   30   60   90  120
> outer(x, y, "+")
     [,1] [,2] [,3] [,4]
[1,]   11   21   31   41
[2,]   12   22   32   42
[3,]   13   23   33   43
> outer(x, y, "-")
     [,1] [,2] [,3] [,4]
[1,]   -9  -19  -29  -39
[2,]   -8  -18  -28  -38
[3,]   -7  -17  -27  -37
> outer(x, y, "-") |> t()
     [,1] [,2] [,3]
[1,]   -9   -8   -7
[2,]  -19  -18  -17
[3,]  -29  -28  -27
[4,]  -39  -38  -37
> f <- function(x, y) {x ^ 2 + y ^ 3}
> outer(x, y, f)
     [,1] [,2]  [,3]  [,4]
[1,] 1001 8001 27001 64001
[2,] 1004 8004 27004 64004
[3,] 1009 8009 27009 64009

［Octave］重み付き非線形最小二乗法を行う

optimパッケージのleasqr関数を使う。

以下の例では、与えられた8組の測定値を使用して、非線形最小二乗法により回帰モデルy = b1 + b2 * x * exp(b3 * x)（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。

>> % 計算の準備
>> pkg load optim
>> x = [1 2 3 5 6 7 8 9]';
>> y = [1 4 10 25 32 49 64 75]';
>> w = repmat(1, 1, length(x))';
>> w(3) = 100;
>> w(5) = 100;
>> f = @(x, b) b(1) + b(2) * x .* exp(b(3) * x);
>> b0 = [1 1 1]';
>>
>> % 重みを考慮しない計算(図の青実線)
>> [yest1, b, irlt, iter, corm, covm] = leasqr(x, y, b0, f, 1.0e-8, 100);
>> disp(b)     # 推定した回帰係数
  -3.8908
   3.3949
   0.1085
>> disp(irlt)  # 収束判定(収束したら1)
1
>> disp(covm)  # 推定した回帰係数の分散共分散行列
   6.0731e+00  -1.5657e+00   4.4934e-02
  -1.5657e+00   5.5200e-01  -1.7014e-02
   4.4934e-02  -1.7014e-02   5.3738e-04
>>
>> % 重みを考慮した計算(図の赤実線)
>> [yest2, b, irlt, iter, corm, covm] = leasqr(x, y, b0, f, 1.0e-8, 100, w);
>> disp(b)     # 推定した回帰係数
   0.2421
   1.9985
   0.1623
>> disp(irlt)  # 収束判定(収束したら1)
1
>> disp(covm)  # 推定した回帰係数の分散共分散行列
   6.8549e-01  -2.3642e-01   1.6110e-02
  -2.3642e-01   8.1785e-02  -5.5772e-03
   1.6110e-02  -5.5772e-03   3.8044e-04
>>
>> % 図の描画
>> plot(x, yest1, 'b', x, yest2, 'r', x, y, 'ko', 'MarkerFaceColor', 'k')

［Octave］計算機イプシロンを求める

計算機イプシロンとは、不等式「1 + epsm > 1」が成り立つ最小の正数epsmのこと。これを手計算で求めてみる。

>> epsm = 1.0;
>> while (1.0 + epsm > 1.0)
> epsm /= 2.0;
> end
>> epsm *= 2.0;
>> disp(epsm)
2.2204e-16

この環境下では約2.2×10^-16であることがわかる。

［R］ウェブサイトから読み取ったHTMLファイルをファイルに出力する

xml2パッケージのwrite_html関数を使う。以下は外務省のウェブサイトのトップページをrvestパッケージのread_html関数で読み取り、それをファイルに出力した例。

> library(rvest)
> library(xml2)
> html <- read_html("https://www.mofa.go.jp")
> tmp <- tempfile(fileext = ".html")
> write_html(html, tmp)
> readLines(tmp) |> head(5)
[1] "<!DOCTYPE html>"                                                        
[2] "<html lang=\"en\">"                                                     
[3] "<head>"                                                                 
[4] "<meta http-equiv=\"Content-Type\" content=\"text/html; charset=UTF-8\">"
[5] "<meta charset=\"UTF-8\">"     

［R］前進代入法で連立方程式を解く

forwardsolve関数を使う。以下は、連立方程式を行列で表記し、上三角行列の係数行列をA、定数項行列をy、求める解をxとおき、x = Ayと考えてxを求めた例。関数だけではなく、solve関数を使った例、逆行列を求めて正規方程式を解いた例も示してあるが、得られた解は同じであることがわかる。

> set.seed(7)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxaa[upper.tri(mxaa)] <- 0
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    9    0    0    0
[2,]    4    8    0    0
[3,]    2    4    2    0
[4,]    1    9    3    1
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxy)
     [,1]
[1,]    6
[2,]    1
[3,]    9
[4,]    3
> forwardsolve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> solve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> print(solve(t(mxaa) %*% mxaa) %*% (t(mxaa) %*% mxy))
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667

forwardsolve関数は、与えられた行列の下三角の成分だけを使って解く。以下のとおり、すべての成分を使って解くsolve関数とは得られた解が違うことがわかる。

> set.seed(7)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    9    3    2    7
[2,]    4    8    5    1
[3,]    2    4    2    5
[4,]    1    9    3    1
> print(mxy)
     [,1]
[1,]    6
[2,]    1
[3,]    9
[4,]    3
> forwardsolve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> solve(mxaa, mxy)
            [,1]
[1,] -0.97516556
[2,]  0.06125828
[3,]  0.49337748
[4,]  1.94370861

［R］ベクトルの要素すべてが同じかどうか調べる

setequal関数を使う。指定した2つのベクトルの要素がまったく同じであればTRUEを返す。集合として等しいかどうかを判定しており、最後の例のとおり、要素の順番が異なっても属する要素が同じであればTRUEを返す。

> chara1 <- c("ルアン・メェイ", "帰忘の流離人", "イレイナ")
> chara2 <- c("ルアン・メェイ", "帰忘の流離人", "イレイナ", "アストラ")
> setequal(chara1, chara2)
[1] FALSE
> chara2 <- c("ルアン・メェイ", "帰忘の流離人", "イレイナ")
> setequal(chara1, chara2)
[1] TRUE
> chara2 <- c("イレイナ", "ルアン・メェイ", "帰忘の流離人")
> setequal(chara1, chara2)
[1] TRUE

［R］ベクトルに特定の要素が含まれているかどうか調べる

is.element関数を使う。第一引数に含まれているかどうか調べたい要素を、第二引数には調べる対象の集合をベクトルで指定する。調べたい要素は複数同時に指定できる。%in%演算子を使っても同じ結果が得られる。

> chara <- c("ルアン・メェイ", "帰忘の流離人", "イレイナ", "アストラ")
> is.element("キャストリス", chara)
[1] FALSE
> is.element("ルアン・メェイ", chara)
[1] TRUE
> is.element(c("アストラ", "ホタル", "帰忘の流離人"), chara)
[1]  TRUE FALSE  TRUE
> c("アストラ", "ホタル", "帰忘の流離人") %in% chara
[1]  TRUE FALSE  TRUE

［R］相関行列を求める

cor関数を使う。以下は、あらかじめ用意した行列を使って相関行列を求めた例。参考に手計算でも求めている。組み込み関数であるvar関数は、標本分散ではなく不偏分散（平均値からの偏差の平方和を「標本数-1」で割った値）が求まることに注意。

> # 「情報量規準による統計解析入門」表12.1(p.154)
> # 1行1標本
> print(xx)
      [,1] [,2] [,3] [,4] [,5]
 [1,]  165   53   86   56   92
 [2,]  160   47   84   52   92
 [3,]  166   55   86   64   89
 [4,]  164   56   90   60   95
 [5,]  168   55   87   56   87
 [6,]  164   54   87   57   92
 [7,]  168   54   94   58   97
 [8,]  169   55   88   57   92
 [9,]  169   53   86   58   93
[10,]  166   56   84   57   90
> # 組み込み関数corにより求めた相関行列
> print(cor(xx))
            [,1]       [,2]      [,3]         [,4]         [,5]
[1,]  1.00000000  0.6173717 0.3620673  0.368208142 -0.043119609
[2,]  0.61737166  1.0000000 0.3341712  0.665307665 -0.107443062
[3,]  0.36206726  0.3341712 1.0000000  0.280441930  0.685251832
[4,]  0.36820814  0.6653077 0.2804419  1.000000000 -0.006370564
[5,] -0.04311961 -0.1074431 0.6852518 -0.006370564  1.000000000
> 
> # 以下、手計算により相関行列を求める、nに標本数(行列の行数)を代入
> n <- nrow(xx)
> # 各列の平均
> dmean <- apply(xx, 2, mean)
> # 各列の標本分散(平均からの偏差の平方和をnで割る、不偏分散ではない)
> dvar <- apply(xx, 2, function(x) (length(x) - 1) / length(x) * var(x))
> # 各列の標本分散の平方根
> dstde <- sqrt(dvar)
> # 列ごとに平均を引いて標本分散で割って標準化する
> xx0 <- sweep(xx, 2, dmean, FUN = "-")
> xx0 <- sweep(xx0, 2, dstde, FUN = "/")
> # 相関行列を求める
> rr <- t(xx0) %*% xx0 / n
> print(rr)
            [,1]       [,2]      [,3]         [,4]         [,5]
[1,]  1.00000000  0.6173717 0.3620673  0.368208142 -0.043119609
[2,]  0.61737166  1.0000000 0.3341712  0.665307665 -0.107443062
[3,]  0.36206726  0.3341712 1.0000000  0.280441930  0.685251832
[4,]  0.36820814  0.6653077 0.2804419  1.000000000 -0.006370564
[5,] -0.04311961 -0.1074431 0.6852518 -0.006370564  1.000000000

［R］後退代入法で連立方程式を解く

backsolve関数を使う。以下は連立方程式を行列で表記し、上三角行列の係数行列をA、定数項行列をy、求める解をxとおき、x = Ayと考えてxを求めた例。backsolve関数だけではなく、solve関数を使った例、逆行列を求めて正規方程式を解いた結果も示したが、得られた解は同じであることがわかる。

> set.seed(6)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxaa[lower.tri(mxaa)] <- 0
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    6    8    5    1
[2,]    0    9    1    3
[3,]    0    0    6    7
[4,]    0    0    0    3
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxy)
     [,1]
[1,]    5
[2,]    7
[3,]    2
[4,]    7
> backsolve(mxaa, mxy)
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333
> solve(mxaa, mxy)
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333
> print(solve(t(mxaa) %*% mxaa) %*% (t(mxaa) %*% mxy))
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333

backsolve関数は、与えられた行列の上三角の成分だけを使って解く。以下のとおり、すべての成分を使って解くsolve関数とは得られた解が違うことがわかる。

> set.seed(6)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    6    8    5    1
[2,]    9    9    1    3
[3,]    3    9    6    7
[4,]    4    7    9    3
> print(mxy)
     [,1]
[1,]    5
[2,]    7
[3,]    2
[4,]    7
> backsolve(mxaa, mxy)
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333
> solve(mxaa, mxy)
           [,1]
[1,]  1.8801598
[2,] -1.3715047
[3,]  0.8322237
[4,]  0.5299601

［R］重み付き非線形最小二乗法を行う

nls関数を使う。weightsオプションに重みを指定すると、その重みを考慮した計算を行う。

以下の例では、与えられた8組の測定値を使用して、非線形最小二乗法により回帰モデルy = b1 + b2 * x * exp(b3 * x)（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。

> # 計算の準備
> x <- c(1, 2, 3, 5, 6, 7, 8, 9)
> y <- c(1, 4, 10, 25, 32, 49, 64, 75)  # 測定値（図の黒丸）
> w <- rep(1, length(x))
> w[c(3, 5)] <- 100
> dtf <- data.frame(x, y, w)
> print(dtf)
  x  y   w
1 1  1   1
2 2  4   1
3 3 10 100
4 5 25   1
5 6 32 100
6 7 49   1
7 8 64   1
8 9 75   1
> f <- function(x, b) b[1] + b[2] * x * exp(b[3] * x)
> b0 <- c(1, 1, 1)
> 
> # 重みを考慮しない計算（図の青実線）
> r1 <- nls(y ~ f(x, b), start = list(b = b0))
> sm <- summary(r1)
> print(sm)       # nls関数による計算結果の概要
Formula: y ~ f(x, b)
Parameters:
   Estimate Std. Error t value Pr(>|t|)   
b1 -3.89078    2.46437  -1.579  0.17521   
b2  3.39494    0.74297   4.569  0.00601 **
b3  0.10847    0.02318   4.679  0.00544 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.502 on 5 degrees of freedom
Number of iterations to convergence: 21 
Achieved convergence tolerance: 3.21e-06
> print(sm$coef)  # 推定した回帰係数,その標準誤差,t値,p値
     Estimate Std. Error   t value    Pr(>|t|)
b1 -3.8907795 2.46436971 -1.578813 0.175211732
b2  3.3949412 0.74297094  4.569413 0.006005113
b3  0.1084657 0.02318143  4.678990 0.005438508
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 31.30051
> print(sum((y - fitted(r1)) ^ 2))
[1] 31.30051
> print(sm$df)                           # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 6.260101
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 2.502019
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            b1          b2            b3
b1  6.07311807 -1.56568510  0.0449342593
b2 -1.56568510  0.55200581 -0.0170136719
b3  0.04493426 -0.01701367  0.0005373789
> 
> # 重みを考慮した計算（図の赤実線）
> r2 <- nls(y ~ f(x, b), weights = w, start = list(b = b0))
> sm <- summary(r2)
> print(sm)       # nls関数による計算結果の概要
Formula: y ~ f(x, b)
Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b1 -0.03045    1.18131  -0.026 0.980435    
b2  2.06079    0.33571   6.139 0.001666 ** 
b3  0.15929    0.02080   7.659 0.000604 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.741 on 5 degrees of freedom
Number of iterations to convergence: 14 
Achieved convergence tolerance: 3.676e-06
> print(sm$coef)  # 推定した回帰係数,その標準誤差,t値,p値
      Estimate Std. Error     t value     Pr(>|t|)
b1 -0.03044686 1.18131366 -0.02577373 0.9804348365
b2  2.06079018 0.33571296  6.13854825 0.0016664488
b3  0.15929178 0.02079896  7.65864353 0.0006043864
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 69.95679
> print(sum(w * (y - fitted(r2)) ^ 2))
[1] 69.95679
> print(sm$df)                           # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 13.99136
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 3.740502
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            b1           b2            b3
b1  1.39550197 -0.371789063  0.0219612500
b2 -0.37178906  0.112703189 -0.0069215507
b3  0.02196125 -0.006921551  0.0004325965
> 
> # 図の描画
> plot(x, fitted(r1), type = "l", col = "blue", xlab = "x", ylab = "y")
> lines(x, fitted(r2), col = "red")
> points(x, y, pch = 20)

［R］パイプラインを使ってリストの要素を取り出す

magrittrパッケージのextract関数かextract2関数を使う。それぞれ|>演算子（%>%演算子）を使用したパイプラインにおいて、extract関数は[ ]の、extract2関数は[[ ]]のエイリアスとして機能する。

> chara <- c("アストラ", "ヴィタ", "リフ")
> seiyu <- c("遠藤綾", "日笠陽子", "瀬戸麻沙美")
> sakuhin <- c("ゼンレスゾーンゼロ", "崩壊3rd", "スノウブレイク")
> lis <- list(sakuhin = sakuhin, chara = chara, seiyu = seiyu)
> print(lis)
$sakuhin
[1] "ゼンレスゾーンゼロ" "崩壊3rd"            "スノウブレイク"    
$chara
[1] "アストラ" "ヴィタ"   "リフ"    
$seiyu
[1] "遠藤綾"     "日笠陽子"   "瀬戸麻沙美"
> lis[3]
$seiyu
[1] "遠藤綾"     "日笠陽子"   "瀬戸麻沙美"
> lis |> extract(3)
$seiyu
[1] "遠藤綾"     "日笠陽子"   "瀬戸麻沙美"
> lis |> extract(3) |> typeof()
[1] "list"
> lis[[3]]
[1] "遠藤綾"     "日笠陽子"   "瀬戸麻沙美"
> lis |> extract2(3)
[1] "遠藤綾"     "日笠陽子"   "瀬戸麻沙美"
> lis |> extract2(3) |> typeof()
[1] "character"

リストに[ ]を使うと、指定した要素だけからなるリストを返す。[[ ]]は指定した要素を取り出して、その要素の自体を返す。

［R］行列の行（列）ごとに演算を行う

apply関数を使う。第二引数に1を指定すると行ごと、2を指定すると列ごとに演算を行い、結果をベクトルで返す。第三引数に演算を行う関数を指定する。指定する関数は、組み込み関数でも自作関数でも、直接書いてもかまわない。

> mx <- matrix(1:8, 2, 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    1    3    5    7
[2,]    2    4    6    8
> apply(mx, 1, sum)
[1] 16 20
> apply(mx, 2, sum)
[1]  3  7 11 15
> f <- function(d) {sum(d ^ 2)}
> apply(mx, 2, f)
[1]   5  25  61 113
> apply(mx, 2, function(d) {sum(d ^ 3)})
[1]   9  91 341 855