［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 * 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] * exp(b[3] * x)
> b0 <- c(1, 1, 1)
> 
> # 重みを考慮しない計算（図の青実線）
> r1 <- nls(y ~ f(x, b), start = list(b = b0))
> sm <- summary(r1)
> print(r1)       # nls関数による計算結果の概要
Nonlinear regression model
  model: y ~ f(x, b)
   data: parent.frame()
      b1       b2       b3 
-25.4828  21.1533   0.1756 
 residual sum-of-squares: 32.98
Number of iterations to convergence: 21 
Achieved convergence tolerance: 3.185e-06
> print(sm$coef)  # 推定した回帰係数,その標準誤差,t値,p値
      Estimate  Std. Error   t value    Pr(>|t|)
b1 -25.4827638 10.22340189 -2.492591 0.054986347
b2  21.1532994  8.08590380  2.616071 0.047321602
b3   0.1755871  0.03299953  5.320898 0.003137863
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 32.98248
> print(sum((y - fitted(r1)) ^ 2))
[1] 32.98248
> print(sm$df)                           # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 6.596497
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 2.568365
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            b1          b2           b3
b1 104.5179462 -81.8524151  0.329318514
b2 -81.8524151  65.3818402 -0.265714306
b3   0.3293185  -0.2657143  0.001088969
> 
> # 重みを考慮した計算（図の赤実線）
> r2 <- nls(y ~ f(x, b), weights = w, start = list(b = b0))
> sm <- summary(r2)
> print(r2)       # nls関数による計算結果の概要
Nonlinear regression model
  model: y ~ f(x, b)
   data: parent.frame()
      b1       b2       b3 
-10.9272  10.1219   0.2413 
 weighted residual sum-of-squares: 70.58
Number of iterations to convergence: 16 
Achieved convergence tolerance: 5.171e-07
> print(sm$coef)  # 推定した回帰係数,その標準誤差,t値,p値
     Estimate Std. Error   t value     Pr(>|t|)
b1 -10.927232 3.53132490 -3.094372 0.0270250777
b2  10.121862 2.44947466  4.132258 0.0090650345
b3   0.241281 0.02632775  9.164514 0.0002593529
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 70.57934
> print(sum(w * (y - fitted(r2)) ^ 2))
[1] 70.57934
> print(sm$df)                           # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 14.11587
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 3.757109
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            b1          b2            b3
b1 12.47025554 -8.57361848  0.0908751718
b2 -8.57361848  5.99992612 -0.0642324913
b3  0.09087517 -0.06423249  0.0006931503
> 
> # 図の描画
> 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

［R］LU分解を行う

matrixcalcパッケージのlu.decomposition関数を使う。n次の正方行列Aが与えられたとき、下三角行列Lと上三角行列Uを用いてA = L Uと分解することを、行列AのLU分解という。戻り値はLが下三角行列、Uが上三角行列。lu.decomposition関数はピボット選択は行わないため、戻り値に置換行列はない。以下は4次の正方行列のパスカル行列をLU分解した例。最後に、元の行列に戻るかどうか計算している。

> n <- 4
> aa <- matrix(1, n, n)
> for (i in 2:n) for (j in 2:n) aa[i, j] <- aa[i, j - 1] + aa[i - 1, j]
> print(aa)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    2    3    4
[3,]    1    3    6   10
[4,]    1    4   10   20
> library(matrixcalc)
> r <- lu.decomposition(aa)
> print(r$L)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    1    0    0
[3,]    1    2    1    0
[4,]    1    3    3    1
> print(r$U)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    0    1    2    3
[3,]    0    0    1    3
[4,]    0    0    0    1
> print(r$L %*% r$U)
     [,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］正規表現でIPv4によるIPアドレス表記の文字列かどうか判定する

正規表現パターンに「^([0-9]+\\.){3}[0-9]+$」を指定すればよい。あらかじめ文字列の先頭と末尾の空白は取り除いておくこと。

> s <- c("12.34.56.78", "123.456.789.123", "123.456", "a.b.c.d", "ab.cd.ef.gh")
> grep("^([0-9]+\\.){3}[0-9]+$", s, value = TRUE)
[1] "12.34.56.78"     "123.456.789.123"

［R］下三角行列を作成する

upper.tri関数に行列を指定すると、対角成分より上の成分（各成分を(i,j)と表記すればi<jの成分）だけTRUEを返す。これを利用してその成分に0を代入すれば、既存の行列を簡単に下三角行列（対角成分より上の成分はすべて0の行列）に変換することができる。

> set.seed(5)
> mx <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    2    1    9    3
[2,]    7    7    1    6
[3,]    9    5    3    3
[4,]    3    8    5    2
> print(upper.tri(mx))
      [,1]  [,2]  [,3]  [,4]
[1,] FALSE  TRUE  TRUE  TRUE
[2,] FALSE FALSE  TRUE  TRUE
[3,] FALSE FALSE FALSE  TRUE
[4,] FALSE FALSE FALSE FALSE
> mx[upper.tri(mx)] <- 0
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    2    0    0    0
[2,]    7    7    0    0
[3,]    9    5    3    0
[4,]    3    8    5    2

［R］上三角行列を作成する

lower.tri関数に行列を指定すると、対角成分より下の成分（各成分を(i,j)と表記すればi>jの成分）だけTRUEを返す。これを利用してその成分に0を代入すれば、既存の行列を簡単に上三角行列（対角成分より下の成分はすべて0の行列）に変換することができる。

> set.seed(4)
> mx <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    6    8    9    1
[2,]    1    3    1    9
[3,]    3    7    7    4
[4,]    3    9    3    5
> print(lower.tri(mx))
      [,1]  [,2]  [,3]  [,4]
[1,] FALSE FALSE FALSE FALSE
[2,]  TRUE FALSE FALSE FALSE
[3,]  TRUE  TRUE FALSE FALSE
[4,]  TRUE  TRUE  TRUE FALSE
> mx[lower.tri(mx)] <- 0
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    6    8    9    1
[2,]    0    3    1    9
[3,]    0    0    7    4
[4,]    0    0    0    5

［R］ディレクトリを作成する

dir.create関数を使う。以下は、一時ディレクトリ（テンポラリーディレクトリ）に、実際にディレクトリを作成した例。

> print(tempdir())
[1] "C:\\Users\\○○○\\AppData\\Local\\Temp\\×××"
> list.files(tempdir(), include.dir = TRUE)
character(0)
> dir.create(file.path(tempdir(), "鈴木みのり"))
> list.files(tempdir(), include.dir = TRUE)
[1] "鈴木みのり"
> dir.create(file.path(tempdir(), "鈴木みのり", "セナディア"))
> list.files(tempdir(), include.dir = TRUE)
[1] "鈴木みのり"
> list.files(file.path(tempdir(), "鈴木みのり"), include.dir = TRUE)
[1] "セナディア"

［R］ウェブサイトの応答ヘッダーの取得に失敗する

curlGetHeaders関数は、ウェブサイトの応答ヘッダーを得ることができるが、外部の一般のウェブサイト（グーグルやヤフー）は読み取れるのにlocalhostやイントラのウェブサイトが読み取れないことがある。この場合、内部のウェブサイトはプロキシを経由しないで直接接続する設定になっていると考えられる。そのため、特定のURLではプロキシを使わないような設定にすればよい。

プロキシを使わずに接続するドメインを、Sys.setenv関数を使って環境変数no_proxyに設定すればよい。複数ある場合はコンマで区切って指定する。

> curlGetHeaders("https://www.google.co.jp") |> head(1)
[1] "HTTP/1.1 200 Connection established\r\n"
> curlGetHeaders("http://localhost") |> head(1)
[1] "HTTP/1.0 403 Forbidden\r\n"
> curlGetHeaders("http://intra.company.co.jp") |> head(1)
[1] "HTTP/1.0 502 Bad Gateway\r\n"
> Sys.setenv("no_proxy" = "localhost,intra.company.co.jp,123.456.78.9")
> curlGetHeaders("http://localhost") |> head(1)
[1] "HTTP/1.1 200 OK\r\n"
> curlGetHeaders("http://intra.company.co.jp") |> head(1)
[1] "HTTP/1.1 200 OK\r\n"
> curlGetHeaders("http://123.456.78.9") |> head(1)
[1] "HTTP/1.1 200 OK\r\n"

URLにIPアドレスを直接指定する場合もあるだろうが、そのような場合read_htmlはまったく応答せずしばらくR自体が停止したような状態になってしまうが、上記の例のとおりにそのまま書き込めば、応答するようになる。