« 2023年9月 | トップページ | 2023年11月 »
数値型ベクトルで連番を作成し、それをsprintf関数で加工して文字列に変換すればよい。連番は、1,2,3,…と001,002,003,…のような2種類があるが、sprintf関数で指定する書式文字列をうまく利用すれば、どちらも簡単に作ることができる。cat関数を使えば、ファイルに簡単に出力できる。
> n <- 1:3
> sprintf("%d", n)
[1] "1" "2" "3"
> sprintf("%04d", n)
[1] "0001" "0002" "0003"
> s <- sprintf("file%04d.jpg", n)
> print(s)
[1] "file0001.jpg" "file0002.jpg" "file0003.jpg"
> cat(s, file = "temp.csv", sep = "\n")
temp.csvの中身
file0001.jpg
file0002.jpg
file0003.jpg
seq関数を使う。第1引数には初項を指定する。第2引数は第3引数の指定により異なる。第3引数はデフォルトではbyオプションであり、byオプションは公差の指定で、初項から加算して第2引数を超えない値を末項にした数列を作成する。第3引数にlengthオプションを指定した場合は、第1引数と第2引数をそれぞれ初項、末項にして、指定した値が項数となる数列を作成する。
第1引数と第2引数は自動で比較され、byオプションで公差を明示しないと自動的に公差の正負が決められてしまうことに注意。第3引数を指定しないと自動的に公差は1,-1,0のいずれかになる。当然、実数も指定することができる。
> seq(1, 5)
[1] 1 2 3 4 5
> seq(1, 10, 3)
[1] 1 4 7 10
> seq(1, 12, 3)
[1] 1 4 7 10
> seq(1, 12, by = 3)
[1] 1 4 7 10
> seq(1, 12, length = 3)
[1] 1.0 6.5 12.0
> seq(1, -2)
[1] 1 0 -1 -2
> seq(-5, -5)
[1] -5
> seq(1, 1, length = 10)
[1] 1 1 1 1 1 1 1 1 1 1
> seq(1.30, 1.22, by = -0.02)
[1] 1.30 1.28 1.26 1.24 1.22
> seq(1.30, 1.22, length = 4)
[1] 1.300000 1.273333 1.246667 1.220000
> dtf <- read.csv("table5_1.csv", header = TRUE)
> n <- nrow(dtf)
> xxi <- dtf$xxi
> ddi <- dtf$ddi
> yyi <- dtf$yyi
> k <- 3
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi, ddi), ncol = 3)
> mxbh <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> mxyh <- mxxx %*% mxbh
> mxu <- mxy - mxyh
> ssy2 <- sum(t(mxy) %*% mxy) - sum(yyi) ^ 2 / n
> ssyh2 <- sum(t(mxyh) %*% mxyh) - sum(yyi) ^ 2 / n
> sse2 <- as.vector(t(mxu) %*% mxu)
> s2 <- sse2 / (n - k)
> s <- sqrt(s2)
> mxvv <- s2 * solve(t(mxxx) %*% mxxx)
> mxs <- sqrt(diag(mxvv))
> rr2 <- ssyh2 / ssy2
> tval <- as.vector(mxbh) / as.vector(mxs)
> # 推定された回帰モデルのパラメーター α^,β^,γ^
> print(as.vector(mxbh))
[1] 0.5 1.1 15.8
> # 決定係数 R^2
> print(rr2)
[1] 0.9936886
> # 各パラメーターの t値
> print(tval)
[1] 0.4303315 5.1854497 13.8309443
> # 図の出力
> png("fig5_1.png", width = 512, height = 512)
> plot(xxi, yyi, asp = 1., type = "n", xlim = c(0, 12), ylim = c(0, 30), xlab = "X", ylab = "Y")
> x <- seq(0, 12, length = 100)
> lines(x, mxbh[1] + mxbh[2] * x, lty = "solid")
> lines(x, mxbh[1] + mxbh[2] * x + mxbh[3], lty = "dashed")
> points(xxi, yyi, pch = 20)
> dev.off()
null device
1
System.Text.RegularExpressions名前空間のRegex.IsMatchメソッドを使う。第1引数に検索対象の文字列を、第2引数に正規表現を指定する。マッチすれば（指定した正引き表現に一致する箇所が見つかれば）trueを返し、そうではない場合はfalseが返される。
> using System.Text.RegularExpressions;
> Regex.IsMatch("ABC", "A")
true
> Regex.IsMatch("ABC", "ABCD")
false
> Regex.IsMatch("ABC", "^A")
true
> Regex.IsMatch("ABC", "A$")
false
> Regex.IsMatch("\t", "^\s+$")
(1,23): error CS1009: 認識できないエスケープ シーケンスです
> Regex.IsMatch("\t", "^\\s+$")
true
> Regex.IsMatch("\t", @"^\s+$")
true
> Regex.IsMatch("\t\t", @"^\s+$")
true
> Regex.IsMatch(" ", @"^\s+$")
true
> Regex.IsMatch(" AB ", @"^\s+$")
false
PNG形式であればpng関数、JPEG形式であればjpeg関数、TIFF形式であればtiff関数を使う。各関数は、いずれも図を描画する際の出力先をファイルに変更する関数であり、図の描画が終わったら、dev.off関数を必ず実行しなければならない。dev.off関数を実行することで、ファイルへの出力が完了する。
出力時の画像の大きさを指定することができる。widthオプションには横方向の大きさ（単位はピクセル）、heightオプションには縦方向の大きさ（同）を指定する。
> x <- seq(0, 6 * pi, length = 100)
> y1 <- sin(x)
> y2 <- sin(8 * x)
> png("pngfunc.png", width = 512, height = 256)
> plot(x, y1 + y2, asp = 1., type = "l")
> dev.off()
null device
1
> jpeg("jpegfunc.jpg", width = 512, height = 256)
> plot(x, y1 + y2, asp = 1., type = "l")
> dev.off()
null device
1
> tiff("tifffunc.tif", width = 512, height = 256)
> plot(x, y1 + y2, asp = 1., type = "l")
> dev.off()
null device
1
公式ページからダウンロードした入力ファイル「Weight.dat」をカレントディレクトリに置いておくこと。
> library(KFAS)
> dtf <- read.table("Weight.dat", header = FALSE)
> yt <- dtf$V1
> # ローカルレベルモデル 灰実線
> mod1tr <- SSModel(yt ~ SSMtrend(1, Q = NA), H = NA)
> fit1tr <- fitSSM(mod1tr, numeric(2), method = "BFGS")
> kfs1tr <- KFS(fit1tr$model)
> aht1tr <- kfs1tr$alphahat[, "level"]
> # 2次のトレンドモデル 黒実線
> mod2tr <- SSModel(yt ~ SSMtrend(2, Q = c(list(0), list(NA))), H = NA)
> fit2tr <- fitSSM(mod2tr, numeric(2), method = "BFGS")
> kfs2tr <- KFS(fit2tr$model)
> aht2tr <- kfs2tr$alphahat[, "level"]
> ahts2tr <- kfs2tr$alphahat[, "slope"]
> # ローカル線形トレンドモデル 黒破線
> modllt <- SSModel(yt ~ SSMtrend(2, Q = c(list(NA), list(NA))), H = NA)
> fitllt <- fitSSM(modllt, numeric(3), method = "BFGS")
> kfsllt <- KFS(fitllt$model)
> ahtllt <- kfsllt$alphahat[, "level"]
> ahtsllt <- kfsllt$alphahat[, "slope"]
> #
> png("kalmanfil_p078_1.png", width = 768, height = 512)
> plot(yt, lty = 3, type = "o", xlab = "経過日数", ylab = "水準成分")
> lines(aht1tr, lwd = 2, col = 8) # ローカルレベルモデル(灰実線)
> lines(aht2tr, lwd = 2) # 2次のトレンドモデル(黒実線)
> lines(ahtllt, lwd = 2, lty = 2) # ローカル線形トレンドモデル(黒破線)
> dev.off()
null device
1
> #
> png("kalmanfil_p078_2.png", width = 768, height = 512)
> plot(ahts2tr, lwd = 2, xlab = "経過日数", ylab = "傾き成分") # 2次のトレンドモデル(黒実線)
> lines(ahtsllt, lwd = 2, lty = 2) # ローカル線形トレンドモデル(黒破線)
> dev.off()
null device
1
> #
> r <- qv <- seta12 <- seta22 <- seps2 <- mll <- aic <- mse <- double(3)
> r <- c(2, 2, 3)
> qv <- c(1, 2, 2)
> seta12[1] <- as.vector(fit1tr$model$Q[, , 1])
> seta12[2] <- as.vector(diag(fit2tr$model$Q[, , 1])[1])
> seta22[2] <- as.vector(diag(fit2tr$model$Q[, , 1])[2])
> seta12[3] <- as.vector(diag(fitllt$model$Q[, , 1])[1])
> seta22[3] <- as.vector(diag(fitllt$model$Q[, , 1])[2])
> seps2[1] <- as.vector(fit1tr$model$H[, , 1])
> seps2[2] <- as.vector(fit2tr$model$H[, , 1])
> seps2[3] <- as.vector(fitllt$model$H[, , 1])
> mll[1] <- kfs1tr$logLik - sum(kfs1tr$Finf > 0) * log(2 * pi) / 2
> mll[2] <- kfs2tr$logLik - sum(kfs2tr$Finf > 0) * log(2 * pi) / 2
> mll[3] <- kfsllt$logLik - sum(kfsllt$Finf > 0) * log(2 * pi) / 2
> aic <- -2 * mll + 2 * (r + qv)
> idx <- 3:60
> n <- length(idx)
> mse[1] <- sum(kfs1tr$v[idx] ^ 2) / n
> mse[2] <- sum(kfs2tr$v[idx] ^ 2) / n
> mse[3] <- sum(kfsllt$v[idx] ^ 2) / n
> dtf <- data.frame(r, q = qv, seta12, seta22, seps2, mll, aic, mse)
> # 表3.2
> print(dtf)
r q seta12 seta22 seps2 mll aic mse
1 2 1 0.07078040 0.000000e+00 0.1491052 -48.58750 103.1750 0.2984077
2 2 2 0.00000000 7.097290e-04 0.2432792 -54.89776 117.7955 0.3581985
3 3 2 0.08655912 5.260979e-07 0.1371066 -51.74947 113.4989 0.3234503
公式ページからダウンロードした入力ファイル「Weight.dat」をカレントディレクトリに置いておくこと。太実線が予測値、細実線が95%信頼区間。破線が95%予測区間。
> dtf <- read.table("Weight.dat", header = FALSE)
> yt <- yt0 <- dtf$V1
> yt <- yt0[1:50]
> n0 <- length(yt0)
> n <- length(yt)
> np <- n0 - n
> t0 <- 1:n0
> t <- 1:n
> tp <- setdiff(t0, t)
> mod <- SSModel(yt ~ SSMtrend(1, Q = NA), H = NA)
> fit <- fitSSM(mod, numeric(2), method = "BFGS")
> kfs <- KFS(fit$model)
> at <- kfs$a[c(-1)]
> ppt <- kfs$P[, , -1] - as.vector(fit$model$Q)
> seta2 <- as.vector(fit$model$Q)
> seps2 <- as.vector(fit$model$H)
> abt <- rep(at[n], np)
> j <- 1:(np)
> ppbt <- rep(ppt[n], np) + j * seta2
> ffbt <- ppbt + seps2
> alp <- 0.05
> abtcon1 <- abt + qnorm(alp / 2) * sqrt(ppbt)
> abtcon2 <- abt + qnorm(1 - alp / 2) * sqrt(ppbt)
> abtpre1 <- abt + qnorm(alp / 2) * sqrt(ffbt)
> abtpre2 <- abt + qnorm(1 - alp / 2) * sqrt(ffbt)
> plot(t0, yt0, type = "n", xlab = "経過日数", ylab = "体重(kg)", ylim = c(83, 87.5))
> lines(t0, yt0, lty = "dotted", col = "gray")
> points(t0, yt0, pch = 21, bg = "white", col = "gray", cex = 1.3)
> lines(tp, abtpre1, lty = "dashed", lwd = 3)
> lines(tp, abtpre2, lty = "dashed", lwd = 3)
> lines(tp, abtcon1, lty = "solid", lwd = 1)
> lines(tp, abtcon2, lty = "solid", lwd = 1)
> lines(tp, abt, lty = "solid", lwd = 2)
> points(t, yt, pch = 21, bg = "white", cex = 1.3)
-encodingオプションを使う。以下は、それぞれUTF-8、シフトJISを指定してコンパイルした例。
>javac -encoding utf-8 PrintText.java
>javac -encoding sjis PrintText.java
コンパイラでサポートされていない文字コードを指定すると、エラーが発生する。
>javac -encoding euc PrintEncode1.java
エラー: サポートされていないエンコーディングです: euc
エラー1個
公式ページからダウンロードした入力ファイル「Weight.dat」をカレントディレクトリに置いておくこと。細実線が95%信頼区間。破線が95%予測区間。
> library(KFAS)
> dtf <- read.table("Weight.dat", header = FALSE)
> yt <- yt0 <- dtf$V1
> n <- nrow(dtf)
> t <- 1:n
> mv <- 21:40
> yt[mv] <- NA
> mod <- SSModel(yt ~ SSMtrend(1, Q = NA), H = NA)
> fit <- fitSSM(mod, numeric(2), method = "BFGS")
> kfs <- KFS(fit$model)
> aht <- kfs$alphahat
> vvt <- as.vector(kfs$V)
> seps2 <- as.vector(fit$model$H)
> alp <- 0.05
> ahtcon1 <- aht + qnorm(alp / 2) * sqrt(vvt)
> ahtcon2 <- aht + qnorm(1 - alp / 2) * sqrt(vvt)
> ahtpre1 <- aht + qnorm(alp / 2) * sqrt(vvt + seps2)
> ahtpre2 <- aht + qnorm(1 - alp / 2) * sqrt(vvt + seps2)
> plot(t, yt, type = "n", xlab = "経過日数", ylab = "体重(kg)")
> lines(t, yt0, lty = "dotted", col = "gray")
> points(t, yt0, pch = 21, bg = "white", col = "gray", cex = 1.3)
> lines(t[mv], ahtpre1[mv], lty = "dashed", lwd = 3)
> lines(t[mv], ahtpre2[mv], lty = "dashed", lwd = 3)
> lines(t, ahtcon1, lty = "solid", lwd = 1)
> lines(t, ahtcon2, lty = "solid", lwd = 1)
> lines(t, aht, lty = "solid", lwd = 2)
> lines(t, yt, lty = "dotted")
> points(t, yt, pch = 21, bg = "white", cex = 1.3)
$文字による文字列補間を使うと、String.Formatメソッドによる変数の展開を、簡単に表すことができる。以下、例。
> string mei = "Kaname";
> string sei = "Buccaneer";
> Console.WriteLine(String.Format("{0} {1}", mei, sei));
Kaname Buccaneer
> Console.WriteLine($"{mei} {sei}");
Kaname Buccaneer
> Console.WriteLine($"{mei, 10} {sei}");
Kaname Buccaneer
> Console.WriteLine($"{mei, 2} {sei}");
Kaname Buccaneer
文字列以外にも例えば浮動小数点数も同じように扱うことができる。
> double d = 1.2345;
> Console.WriteLine(String.Format("{0:F5}", d));
1.23450
> Console.WriteLine($"{d:F5}");
1.23450
> Console.WriteLine($"{d, 9:F5}");
1.23450
> Console.WriteLine($"{d:F5} {d:F3}");
1.23450 1.235
以下はWindows版における例。以下をPrintEncode1.javaと保存して実行する。文字コードはシフトJISでもUTF-8でもどちらでもよい（1バイト文字しかないためどちらにしても同一のファイルになる）。
class PrintEncode1 {
public static void main(String[] args) {
System.out.println(System.getProperty("file.encoding"));
}
}
コンパイルおよび実行する。
>javac -version
javac 21
>java -version
openjdk version "21" 2023-09-19
OpenJDK Runtime Environment (build 21+35-2513)
OpenJDK 64-Bit Server VM (build 21+35-2513, mixed mode, sharing)
>javac PrintEncode1.java
>java PrintEncode1
UTF-8
上記のバージョンでは、デフォルトの文字コードはUTF-8であることがわかる。
このバージョンではデフォルトはUTF-8であるが、ソースコードにシフトJISを使うことは可。以下をシフトJISでPrintEncode2.javaと保存して実行してみる。コンパイルには-encodingオプションを使いsjisを指定する。
class PrintEncode2 {
public static void main(String[] args) {
System.out.println("ABC123あいう");
}
}
コンパイルおよび実行する。
>javac -encoding sjis PrintEncode2.java
>java PrintEncode2
ABC123あいう
ソースコードはソフトJISだが、正しく動作していることがわかる。なお、UTF-8で強制的にコンパイルしようとすると、エラーが発生する。
>javac -encoding utf-8 PrintEncode2.java
PrintEncode2.java:3: エラー: この文字(0x82)は、エンコーディングutf-8にマップできません
System.out.println("ABC123??????");
^
PrintEncode2.java:3: エラー: この文字(0xA0)は、エンコーディングutf-8にマップできません
System.out.println("ABC123??????");
^
（以下、表示省略）
公式ページからダウンロードした入力ファイル「Weight.dat」をカレントディレクトリに置いておくこと。
> dtf <- read.table("Weight.dat", header = FALSE)
> yt <- dtf$V1
> n <- nrow(dtf)
> tn <- 1:n
> #
> mod <- SSModel(yt ~ SSMtrend(1, Q = NA), H = NA)
> fit <- fitSSM(mod, numeric(2), method = "BFGS")
> kfs <- KFS(fit$model)
> aht <- kfs$alphahat
> vvt <- as.vector(kfs$V)
> alp <- 0.05
> ahtcon1 <- aht + qnorm(alp / 2) * sqrt(vvt)
> ahtcon2 <- aht + qnorm(1 - alp / 2) * sqrt(vvt)
> #
> png("kalmanfil_p033.png", width = 768, height = 512)
> plot(tn, yt, type = "n", ylim = c(83, 87), xlab = "経過日数", ylab = "体重(kg)")
> lines(tn, ahtcon1, lty = "solid", lwd = 1)
> lines(tn, ahtcon2, lty = "solid", lwd = 1)
> lines(tn, aht, lty = "solid", lwd = 2)
> lines(tn, yt, lty = "dotted")
> points(tn, yt, pch = 21, bg = "white", cex = 1.3)
> dev.off()
null device
1
公式ページからダウンロードした入力ファイル「Weight.dat」をカレントディレクトリに置いておくこと。
> dtf <- read.table("Weight.dat", header = FALSE)
> yt <- dtf$V1
> n <- nrow(dtf)
> tn <- 1:n
> #
> mod <- SSModel(yt ~ SSMtrend(1, Q = NA), H = NA)
> fit <- fitSSM(mod, numeric(2), method = "BFGS")
> kfs <- KFS(fit$model)
> at <- kfs$a[c(-1)]
> ppt <- kfs$P[, , -1] - as.vector(fit$model$Q)
> alp <- 0.05
> atcon1 <- at + qnorm(alp / 2) * sqrt(ppt)
> atcon2 <- at + qnorm(1 - alp / 2) * sqrt(ppt)
> #
> png("kalmanfil_p029.png", width = 768, height = 512)
> plot(tn, yt, type = "n", ylim = c(83, 87), xlab = "経過日数", ylab = "体重(kg)")
> lines(tn, atcon1, lty = "solid", lwd = 1)
> lines(tn, atcon2, lty = "solid", lwd = 1)
> lines(tn, at, lty = "solid", lwd = 2)
> lines(tn, yt, lty = "dotted")
> points(tn, yt, pch = 21, bg = "white", cex = 1.3)
> dev.off()
null device
1
範囲の計算には、max関数とmin関数を組み合わせて使っても結果は同じ。
> x <- c(51, 49, 50, 57, 43) # 測定値
> print(max(x)) # 最大値
[1] 57
> print(min(x)) # 最小値
[1] 43
> print(range(x)[2] - range(x)[1]) # 範囲
[1] 14
Rに標準で搭載のvar関数は不偏分散を求めるものであり、偏差平方和を標本数で割ったものではない。そこで、最初にvarpという自作関数を使用して、途中で使用している。
> varp <- function(val) {
+ return((length(val) - 1) / length(val) * var(val))
+ }
> dtf <- read.csv("table12_1.csv", header = TRUE)
> mxxx <- as.matrix(dtf[, 2:6])
> n <- nrow(mxxx)
> p <- ncol(mxxx)
> dmean <- apply(mxxx, 2, mean) # 各列の平均
> dvarp <- apply(mxxx, 2, varp) # 各列の分散(不偏分散ではない)
> dstde <- sqrt(dvarp) # 各列の標準偏差
> mxtemp <- sweep(mxxx, 2, dmean, FUN = "-")
> mxxx0 <- sweep(mxtemp, 2, dstde, FUN = "/")
> mxrr <- t(mxxx0) %*% mxxx0 / n
> eigval <- eigen(mxrr)$values
> d1 <- eigval / sum(eigval)
> d2 <- double(p)
> for (i in 1:p) d2[i] <- sum(d1[1:i])
> eigvec <- eigen(mxrr)$vectors
> rownames(eigvec) <- c("身長", "体重", "バスト", "ウェスト", "ヒップ")
> colnames(eigvec) <- c("x1", "x2", "x3", "x4", "x5")
> dtf <- t(data.frame(固有値 = eigval, 寄与率 = d1, 累積寄与率 = d2))
> colnames(dtf) <- c("x1", "x2", "x3", "x4", "x5")
> mxxxa <- mxxx0 %*% eigvec
> colnames(mxxxa) <- c("第1", "第2", "第3", "第4", "第5")
> # 相関行列
> print(mxrr)
height weight bust waist hip
height 1.00000000 0.6173717 0.3620673 0.368208142 -0.043119609
weight 0.61737166 1.0000000 0.3341712 0.665307665 -0.107443062
bust 0.36206726 0.3341712 1.0000000 0.280441930 0.685251832
waist 0.36820814 0.6653077 0.2804419 1.000000000 -0.006370564
hip -0.04311961 -0.1074431 0.6852518 -0.006370564 1.000000000
> # 固有ベクトル
> print(eigvec)
x1 x2 x3 x4 x5
身長 -0.4891523 -0.1816188 0.71013804 -0.4612335 0.1034999
体重 -0.5504555 -0.2906559 -0.06017657 0.6881357 0.3679210
バスト -0.4417133 0.5304590 0.08611513 0.2402352 -0.6770332
ウェスト -0.4869700 -0.2144097 -0.68751116 -0.4820622 -0.1087803
ヒップ -0.1596195 0.7451010 -0.10952321 -0.1537397 0.6194472
> # 固有値・他
> print(dtf)
x1 x2 x3 x4 x5
固有値 2.3741918 1.5421066 0.6417648 0.26079368 0.18114304
寄与率 0.4748384 0.3084213 0.1283530 0.05215874 0.03622861
累積寄与率 0.4748384 0.7832597 0.9116127 0.96377139 1.00000000
> # 第1～5主成分
> print(mxxxa)
第1 第2 第3 第4 第5
[1,] 0.77600803 0.06769004 0.09345140 0.07497037 0.21346159
[2,] 4.00845539 1.03038774 -0.20956461 -0.22943730 -0.24070481
[3,] -1.01409374 -1.65269885 -1.45873809 -0.69958763 -0.43883604
[4,] -1.17933206 1.06949076 -1.19162979 0.58640025 0.19802651
[5,] -0.07964604 -1.56295274 1.07839702 0.47967860 -0.76054108
[6,] 0.41384431 0.13304684 -0.40348070 0.44476269 0.04544139
[7,] -1.88001803 2.48238435 0.43763819 -0.10261088 -0.36915114
[8,] -0.88302667 -0.13708737 0.93628588 -0.05892987 0.14806400
[9,] -0.35297587 -0.07669839 0.64665696 -1.00651873 0.52352423
[10,] 0.19078467 -1.35356239 0.07098373 0.51127250 0.68071535
一部の固有ベクトルと主成分は、書籍p.164の表12.5の値と符号が異なるが、固有ベクトルは0以外の実数をかけたものもまた固有ベクトルであることから、問題はない。
統計学における「範囲」とは、データのとる値の最大値から最小値を引いた差のこと。Rにはこれを直接求める関数は標準で搭載されていないが、それに近い動作をする関数はある。
range関数は、引数に与えたベクトルの最小値と最大値をベクトルで返す。これを利用して範囲を求めることができる。ただしmax関数とmin関数を組み合わせた計算と手順に大差はない。
> d <- c(161, 165, 160, 164, 165, 163, 164, 162, 164, 170)
> range(d)
[1] 160 170
> print(range(d)[2] - range(d)[1]) # 範囲
[1] 10
> print(max(d) - min(d)) # 範囲
[1] 10
eigen関数を使う。
> # 例1
> d <- c(1, 3, -2, -4)
> mx <- matrix(d, 2, 2)
> print(mx)
[,1] [,2]
[1,] 1 -2
[2,] 3 -4
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1] -2 -1
$vectors
[,1] [,2]
[1,] 0.5547002 0.7071068
[2,] 0.8320503 0.7071068
> # 固有値（合計2つ）
> print(eig$value)
[1] -2 -1
> # 固有値 -2 のときの固有ベクトル
> print(eig$vector[, 1])
[1] 0.5547002 0.8320503
> # 固有値 -1 のときの固有ベクトル
> print(eig$vector[, 2])
[1] 0.7071068 0.7071068
>
> # 例2
> d <- c(1, 1, 1, 0, 1, 0, 2, 1, 0)
> mx <- matrix(d, 3, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 1 0 2
[2,] 1 1 1
[3,] 1 0 0
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1] 2 1 -1
$vectors
[,1] [,2] [,3]
[1,] 0.5345225 0 -0.7071068
[2,] 0.8017837 1 0.0000000
[3,] 0.2672612 0 0.7071068
> # 固有値（合計3つ）
> print(eig$value)
[1] 2 1 -1
> # 固有値 2 のときの固有ベクトル
> print(eig$vector[, 1])
[1] 0.5345225 0.8017837 0.2672612
> # 固有値 1 のときの固有ベクトル
> print(eig$vector[, 2])
[1] 0 1 0
> # 固有値 -1 のときの固有ベクトル
> print(eig$vector[, 3])
[1] -0.7071068 0.0000000 0.7071068
固有値は数値的に求めているため、基本的に戻り値は近似値である。そのため、解析的には整数で得られる場合でも戻り値が複素数になる場合がある。この場合、固有ベクトルも複素数になる。以下の例では、解析的には固有値は3と-1（二重根）と求められるが、戻り値は複素数になっている。虚部が無視できるくらい小さいため、虚部を切り捨てて実数ととして扱ってかまわないが、ケースバイケースであることに注意。
> d <- c(1, 2, -1, 2, 1, -1, -4, 4, -1)
> mx <- matrix(d, 3, 3)
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1] 3+0i -1+0i -1-0i
$vectors
[,1] [,2] [,3]
[1,] -0.9045340+0i -7.071068e-01+0.000000e+00i -7.071068e-01+0.000000e+00i
[2,] -0.3015113+0i 7.071068e-01-0.000000e+00i 7.071068e-01+0.000000e+00i
[3,] 0.3015113+0i 0.000000e+00+2.637113e-09i 0.000000e+00-2.637113e-09i
p.41の表に掲載の値をtable_p041.csvに保存しておく。
> library(tidyverse)
> tbl <- read_csv("table_p041.csv", progress = FALSE, show_col_types = FALSE)
> x <- tbl$x
> n <- length(tbl$x)
> xm <- sum(x) / n
> s2 <- sum((x - xm) ^ 2) / n
> s <- sqrt(s2)
> # 分散
> print(s2)
[1] 20
> # 標準偏差
> print(s)
[1] 4.472136
以下をファイルtable_p041.csvに保存する。
no, x
1, 51
2, 49
3, 50
4, 57
5, 43
計算する。
> library(tidyverse)
> tbl <- read_csv("table_p041.csv", progress = FALSE, show_col_types = FALSE)
> # 資料
> print(as.data.frame(tbl))
no x
1 1 51
2 2 49
3 3 50
4 4 57
5 5 43
> cat("偏差\n")
偏差
> no <- tbl$no
> x <- tbl$x
> n <- length(x)
> xm <- sum(x) / n
> print(data.frame(番号 = no, 偏差 = x - xm))
番号 偏差
1 1 1
2 2 -1
3 3 0
4 4 7
5 5 -7
NuGetを利用して.NET用のImageMagickのパッケージをインストールすればよい。とりあえず、以下の3つのパッケージをインストールすればよい。
Magick.NET.Core
Magick.NET.SystemDrawing
Magick.NET-Q16-AnyCPU
以下、パッケージのインストール手順。
メニュー「ツール」→「NuGetパッケージマネージャー」→「パッケージマネージャーコンソール」。中央下に「パッケージマネージャーコンソール」ウィンドウが表示される。次のコマンドを入力すると、Magick関係のパッケージが表示される。
PM> Find-Package Magick
Id Versions Description
-- -------- -----------
Magick.NET.Core {○○○} A .NET API to the ImageMagick image-processing library for Desktop and Web.
Magick.NET-Q16-AnyCPU {○○○} A .NET API to the ImageMagick image-processing library for Desktop and Web.
（以下、表示省略）
上記の3つのパッケージをインストールする。
PM> Install-Package Magick.NET.Core
'.NETFramework,Version=○○○' を対象とするプロジェクト '○○○' に関して、パッケージ 'Magick.NET.Core.○○○' の依存関係情報の収集を試行しています
依存関係情報の収集に ○ ms かかりました
（表示省略）
'Magick.NET.Core ○○○' が ○○○ に正常にインストールされました
NuGet の操作の実行に ○ sec かかりました
経過した時間: ○○○
PM> Install-Package Magick.NET.SystemDrawing
（表示省略）
PM> Install-Package Magick.NET-Q16-AnyCPU
（表示省略）
これでインストールされたはず。ソリューションエクスプローラーのプロジェクト名のところで右クリック。表示されたコンテキストメニューで「NuGetパッケージの管理」をクリック。中央上部に「NuGetパッケージマネージャ」が表示される。ここに、3つのパッケージが表示されていれば、パッケージのインストールは完了。NuGetパッケージマネージャを閉じる（タブのラベルの右端の×をクリック）。これで使えるようになったはず。
動作確認をしてみる。実際に使うときは、以下の宣言が必要。
Imports ImageMagick
例えば、Visual Basicで標準では扱えないwebpファイル（以下の例ではD:\image.webp）は、次のようにすればBitmapオブジェクトに変換することができる。
Dim bitmapobj As System.Drawing.Bitmap
Dim magimg As ImageMagick.MagickImage
magimg = New ImageMagick.MagickImage("D:\image.webp")
bitmapobj = magimg.ToBitmap()
ApplicationクラスのExecutablePathプロパティとPathクラスのGetFileNameメソッドを組み合わせて使う。以下のソースコードをコンパイルして実行してみる。
using System;
using System.IO;
using System.Windows.Forms;
class Exename {
static void Main() {
Console.WriteLine(Path.GetFileName(Application.ExecutablePath));
}
}
実行してみる。
>exename.exe
exename.exe
ApplicationクラスのStartupPathプロパティを使う。以下のソースコードをコンパイルして実行してみる。
using System;
using System.Windows.Forms;
class Supath {
static void Main() {
Console.WriteLine(Application.StartupPath);
}
}
実行してみる。
D:\work>supath.exe
D:\work
D:\work>cd ..
D:\>work\supath.exe
D:\work
source関数を実行したときに、ファイルの文字コードを正しく指定していない場合、このエラーメッセージが表示される。encodingオプションに文字コードを正しく指定すればよい。以下の例では、PowerShellで文字コードをシフトJISにしたスクリプトを作成し、それを実行してみる。
PS > @'
>> cat("ABC123あいう\n")
>> '@ | Out-File -Encoding default temp.R
PS > Format-Hex .\temp.R
パス: ○○○
00 01 02 03 04 05 06 07 08 09 0A 0B 0C 0D 0E 0F
00000000 63 61 74 28 22 41 42 43 31 32 33 82 A0 82 A2 82 cat("ABC123 ¢
00000010 A4 5C 6E 22 29 0D 0A ¤\n")..
> Get-Content .\temp.R
cat("ABC123あいう\n")
Rでスクリプトを動作させてみる。以下のRの環境下では、Sys.getlocale関数の戻り値のとおり、文字コードはUTF-8である。
> Sys.getlocale()
[1] "LC_COLLATE=Japanese_Japan.utf8;LC_CTYPE=Japanese_Japan.utf8;LC_MONETARY=Japanese_Japan.utf8;LC_NUMERIC=C;LC_TIME=Japanese_Japan.utf8"
> source("temp.R")
エラー: invalid multibyte character in parser (temp.R:1:12)
> source("temp.R", encoding = "SJIS")
ABC123あいう
Windowsの日本語版Rの起動時の文字コードは長らくシフトJISであったが、バージョン4.2以降はUTF-8に変わった。R起動時に使用できる文字コードを自動的にシフトJISに変更したい場合は、.RProfileファイルを使用する。
バージョン4.2以降は、デフォルトの文字コードはUTF-8。0x82a0（シフトJISで「あ」）を表示してみる。
'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。
> rawToChar(as.raw(c(0x82, 0xa0)))
[1] "\x82\xa0"
表示されない。文字コードはシフトJIS以外であることがわかる。
Rはスタートメニューのアイコンから起動する場合は、フォルダーC:\Users\○○○\Documentsの.RProfileを自動的に読み込む。テキストエディタで作成し、以下の1行を追加してRを再起動してみる。
Sys.setlocale(locale = "Japanese_Japan.932")'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。
[1] "LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932"
警告メッセージ:
Sys.setlocale(locale = "Japanese_Japan.932") で:
using locale code page other than 65001 ("UTF-8") may cause problems
>
コンソールの最初の表示が変わる。文字コードを確認してみる。
> rawToChar(as.raw(c(0x82, 0xa0)))
[1] "あ"
環境下の文字コードがシフトJISであることがわかる。
.RDataファイルをダブルクリックしてRを起動している場合は、その.RDataファイルと同じフォルダーに置かれた.RProfileを自動的に読み込むため、テキストエディタで作成し、上記と同様にその1行を追加してRを再起動すればよい。
par関数に"usr"パラメーターを指定して実行する。実行すると数値型ベクトルで（横軸方向左端、横軸方向右端、縦軸方向下端、縦軸方向上端）の値を返す。以下、実行例。
> x <- 1:10
> y <- log(x)
> plot(x, y)
> par("usr")
[1] 0.6400000 10.3600000 -0.0921034 2.3946885
> plot(x, y, asp = 1.0)
> par("usr")
[1] 0.640000 10.360000 -3.202457 5.505042
> plot(x, y, asp = 1.0, xaxs = "i", yaxs = "i")
> par("usr")
[1] 1.000000 10.000000 -2.879957 5.182542
plot関数はaspオプションを指定しないと、与えたプロット点に応じてアスペクト比が決まり描画範囲も決まる。例えばaspオプションに1を与えると、描画範囲が短いほうの軸はアスペクト比に応じて強制的に描画範囲が伸ばされる。上の2番目の例では縦軸方向の描画範囲が伸ばされている。
xaxsとyaxsオプションを指定しないと（デフォルトの"r"）、それぞれのプロット点による描画範囲の4%だけ、両方向に描画範囲が広がる（横軸であれば10-1=9より9×0.04=0.36なので、左端は1-0.36=0.64）。この描画範囲を与えたプロット点の範囲ぴったりにしたい場合は、xaxsとyaxsオプションにそれぞれ"i"を与えればよい。上の3番目の例では、図もpar関数の戻り値もそのようになっている。
簡単に求めるには、以下のようにseq関数とlength関数を組み合わせればよい。差が1年未満であれば0、1年～2年未満であれば1、…が得られる。
> length(seq(as.Date("2020-02-15"), as.Date("2020-02-15"), by = "year")) - 1
[1] 0
> length(seq(as.Date("2020-02-15"), as.Date("2020-02-16"), by = "year")) - 1
[1] 0
> length(seq(as.Date("2020-02-15"), as.Date("2021-02-14"), by = "year")) - 1
[1] 0
> length(seq(as.Date("2020-02-15"), as.Date("2021-02-15"), by = "year")) - 1
[1] 1
> length(seq(as.Date("2020-02-15"), as.Date("2021-02-16"), by = "year")) - 1
[1] 1
> length(seq(as.Date("2020-02-15"), as.Date("2031-02-14"), by = "year")) - 1
[1] 10
> length(seq(as.Date("2020-02-15"), as.Date("2031-02-15"), by = "year")) - 1
[1] 11
> length(seq(as.Date("2020-02-15"), as.Date("2031-02-16"), by = "year")) - 1
[1] 11
readrパッケージのwrite_delim関数は、デフォルトでは出力の経過を画面表示しようとするため、その影響で、使用後は画面表示が乱れてカーソルがおかしなところに飛んでしまう。この画面表示の乱れが起きないようにするには、progressオプションにFALSEを指定する。
> library(tidyverse)
> tbl <- tibble(no = 1:3, name = c("ABC", "abc", "123"))
> write_delim(tbl, "temp.csv", delim = ",")
←カーソルがコマンドラインのプロンプトから離れた箇所に表示され、表示が乱れる。
一度Rを終了し、再起動して次のコマンドを実行する。
> library(tidyverse)
> tbl <- tibble(no = 1:3, name = c("ABC", "abc", "123"))
> write_delim(tbl, "temp.csv", delim = ",", progress = FALSE)
←カーソルがコマンドラインのプロンプトのすぐ右側に表示され、表示が乱れない。
以下を.RProfileに記述しておけば、R起動時に自動的に読み込み、デフォルトでそのような設定になる。
options(
readr.show_progress = FALSE
)
R起動後でも、以下のコマンドを実行すれば同じようにデフォルトでそのような設定になる。
> options(readr.show_progress = FALSE)
なお、これらはwrite_csv関数などでも同様である。
sum関数とis.na関数を組み合わせると簡単に計算できる。
> d <- c(1, 2, NA, 4, NA)
> is.na(d)
[1] FALSE FALSE TRUE FALSE TRUE
> sum(is.na(d))
[1] 2
> sum(!is.na(d))
[1] 3
sum関数は論理型ベクトルを与えると、TRUEを1、FALSEを0と見なして計算するため（「Logical true values are regarded as one, false values as zero.」）、is.na関数でNAか否かの論理型ベクトルに変換してそれをsum関数に与えれば、ベクトルに含まれるNAの個数を簡単に求めることができる。
mutate関数を使う。なお、簡潔な表示にするため、tibbleの表示はデータフレームに変換して行っている。
> library(tidyverse)
> tbl <- read_csv("meibo.csv", show_col_types = FALSE, progress = FALSE)
> as.data.frame(tbl)
no sei mei age weight
1 1 あいう えお 11 101.1
2 2 かきく けこ 22 202.2
3 3 さしす せそ 22 303.3
4 4 あいう けこ 33 404.4
> as.data.frame(mutate(tbl, age2 = age * 2, seimei = paste0(sei, mei)))
no sei mei age weight age2 seimei
1 1 あいう えお 11 101.1 22 あいうえお
2 2 かきく けこ 22 202.2 44 かきくけこ
3 3 さしす せそ 22 303.3 44 さしすせそ
4 4 あいう けこ 33 404.4 66 あいうけこ
> as.data.frame(tbl %>% mutate(seibetsu = c("F", "M", "F", "F")))
no sei mei age weight seibetsu
1 1 あいう えお 11 101.1 F
2 2 かきく けこ 22 202.2 M
3 3 さしす せそ 22 303.3 F
4 4 あいう けこ 33 404.4 F
