Видео

今回の動画はちょっと難しいかな？　でも、機械学習を学ぶ基本！　まったりきいてね！

◎ GMMクラスタリングで異常検知 library(mclust) library(car) # データの準備 X <- Davis[-12, c("weight" , "height")] # GMMによるクラスタリングの実行 result <- Mclust(X) # 混合比を取り出す pi <- result$parameters$pro # 元のデータセットの再取得 X <- Davis[, c("weight" , "height")] # クラスタごとの正規分布の確率密度を計算 XX <- cdens(modelName = result$modelName, data = X, parameters = result$parameters) # 異常度を計算 a <- -log(as.matrix(XX) %*% as.matrix(pi)) # 異常度のプロット plot(a, type = "b", col = "blue", pch = 19, xlab = "Sample number", ylab = "Anomaly score", main = "Anomaly Score")

◎GMMでクラスタリング # 必要なパッケージの読み込み library(car) library(mclust) # 観測データを取得 X <- Davis[-12, c("weight", "height")] # Mclust関数を使って最適なモデルを選択 model <- Mclust(X) # BICで最適なコンポーネント数を取得 best_model <- model$G cat("BICが最小となるコンポーネント数:", best_model, " ") # 選択されたモデルの概要を出力 print(summary(model, parameters=TRUE)) # 結果のプロット plot(model)

# Rのプログラム # 正規性の検定(QQ プロット) # グループ1とグループ2のデータを抽出 group1 <- sleep$extra[sleep$group == 1]　 group2 <- sleep$extra[sleep$group == 2] sleep # 対になるデータの差を計算 diff <- group1 - group2 # グラフを1行2列に設定 par(mfrow = c(1, 2)) # 差のヒストグラムをプロット hist(diff, main = "Difference between Group1 and Group2", xlab = "Difference (Group1 - Group2)", col = "lightblue", border = "black") # 差のQQプロットをプロット qqnorm(diff, main = "QQ Plot of Differences", col = "blue") qqline(diff, col = "red", lwd = 2) # 理論的な正規分布線を追加   # データセット sleep の group 列に基づいてグループ1とグループ2の extra データを抽出 group1 <- sleep$extra[sleep$group == 1] group2 <- sleep$extra[sleep$group == 2] # ウィルコクソンの符号付き順位検定を実行（近似的なp値を使用） test_result <- wilcox.test(group1, group2, paired = FALSE, exact = FALSE) # 検定結果を表示 print(test_result)

わかりみサイエンスの原点にもどって、わかりやすくまとめてみました！

# ガンマ関数 n=0:5 factorial = gamma(n+1) data.frame(n, factorial) # 補正係数 library(gt) n <- 2:10 c=function(n){ return(sqrt((n-1)/2)*gamma((n-1)/2)/(gamma(n/2))) } cn=c(n) options(digits = 3) data.frame(n, cn) %>% gt()

08:48 の式の最後　(1-9)log2p ⇒ (1-9)log2(1-p) です。Rのプログラムはあっています。

# Rのプログラム # 必要なライブラリを読み込む library(ggplot2) # pの値を0から1まで0.01刻みで生成 p_values <- seq(0.01, 0.99, by = 0.01) # 0と1を除外 # エントロピーを計算 entropy_values <- -p_values * log2(p_values) - (1 - p_values) * log2(1 - p_values) # データフレームを作成 data <- data.frame(p = p_values, entropy = entropy_values) # グラフを描画 ggplot(data, aes(x = p, y = entropy)) + geom_line(linewidth = 2.0) + labs(x = "確率 P", y = "情報エントロピー H", size = 20) + theme_minimal() + theme(text = element_text(size = 20)) + theme(axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20))

# Rのプログラム ###### (1) 単回帰の概要 ＃ 00:01:02 # 表のRのプログラム # install.packages("Ecdat") # install.packages("gt") # パッケージの読み込み library(Ecdat) library(gt) # データの準備 N <- 1:30 Ice <- cbind(N, Icecream) head(Ice) # 必要な列の選択 Ice0 <- Ice[, c("N", "cons", "temp")] # データの分割 Ice1 <- Ice0[1:10, ] Ice2 <- Ice0[11:20, ] Ice3 <- Ice0[21:30, ] # データフレームの結合 tb <- data.frame(Ice1, Ice2, Ice3) # gt() を使ってきれいな表を作成 tb %>% gt() %>% tab_header( title = "Package of Icecreame", subtitle = "cons vs temp" ) %>% tab_source_note( source_note = md("1951年から1953年にアメリカで30回行われたアイスクリームの消費量 (消費量 vs 気温)") ) library(Ecdat) # 00:01:26 # Rで回帰分析 # install.packages("Ecdat") library(Ecdat) cons<-Icecream$cons ; cons temp<-Icecream$temp ; temp # 変数 x, yにデータを格納（コピペして汎用化するため） x <- temp 　　　　 y <- cons # きれいな図を描くための設定 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) 　 plot(x, y, xlab = "気温(temp)", ylab = "消費量(cons)", cex.lab = 2, cex.axis = 1.5, lwd = 4) # 散布図の描画 # 回帰分析表の出力 S_reg <- lm(y~x) summary(S_reg ) # 回帰直線の描画 abline(S_reg, lwd=4, lty=2, col="red") # 平均値の描画 mx = mean(x) my =mean(y) abline(v=mx, lwd=2) abline(h=my, lwd=2) ##### (2) 決定係数 # 00:12:00 # install.packages("Ecdat") library(Ecdat) x<-Icecream$temp y<-Icecream$cons # 回帰分析表の出力 S_reg <- lm(y~x) summary(S_reg)   ##### (3) 分散分析 # 00:17:40 library(Ecdat) y<-Icecream$cons x<-Icecream$temp S_reg <- lm(y~x) summary(S_reg ) R2 <- cor(x,y)^2 n<- length(x) (F0 <- R2/(1-R2)*(n-2)) (Fq<- qf(1-0.001, 1, n-2)) (p値<-pf(F0, 1, (n-2), lower.tail=FALSE)) # 00:18:04 # 例題で示した図 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) # packages <-c("Ecdat", "tidyverse") # install.packages(packages) library(Ecdat) y<-Icecream$cons x<-Icecream$temp R2 <- cor(x, y)^2 # 決定係数 n<- length(x)　# サンプルサイズ # F分布 f1=1 ; f2=28　　# 自由度 a=0.001 xmax=50 ymax=0.002 curve(df(x,1,28),0,xmax, xlab="F", ylab="df(F,1,28)", ylim=c(0, ymax),cex.lab = 2, cex.axis = 1.5, lwd = 4) # F値０～50 # F値 (F <- R2/(1-R2)*(n-2)) # 検定統計量 (up= qf(1-a, f1, f2)) # 棄却域境界の分位点 F(1,n-2)(0.1%有意水準) abline(v=up, lty=2, col="red", cex.lab = 2, cex.axis = 1.5, lwd = 4) abline(v=F, lty=2, col="blue",cex.lab = 2, cex.axis = 1.5, lwd = 4) # 色を塗る　 xvals <- seq(up,xmax, length=2000) dvals <- df(xvals,f1,f2) polygon(c(xvals,rev(xvals)), c(rep(0,2000),rev(dvals)), col="skyblue") ##### (4) 検出力とサンプルサイズ # 00:21:32 # Retrieve Cohen's suggested effect sizes library("pwr") cohen.ES(test = "f2", size = "small") cohen.ES(test = "f2", size = "medium") cohen.ES(test = "f2", size = "large") # 00:21:59 f <- function(x){ a <- sqrt(x/(1+x)) return(a) } a1=f(0.02) ; a2=f(0.15) ; a3=f(0.35) a1<-paste0("小=", round(a1, 2)) a2<-paste0("中=", round(a2, 2)) a3<-paste0("大=", round(a3, 2)) c(a1, a2, a3) # 00:23:56 # Power=0.8 サンプルサイズ # install.packages("pwr") library (pwr) v= round(pwr.f2.test(u=1, f2=0.02, sig.level=0.05, power=0.8)$v, 1) paste0 ("v=" ,v) # 00:24:37 # Power=0.9 # サンプルサイズ library(pwr) a <- 0.05 po <- 0.9 f <- function(x) { a <- ceiling(pwr.f2.test(u = 1, f2 = x, sig.level = a, power = po)$v) + 2 return(a) } a1 <- f(0.02) a2 <- f(0.15) a3 <- f(0.35) a1 <- paste0("小(0.02)=", a1) a2 <- paste0("中(0.15)=", a2) a3 <- paste0("大(0.35)=", a3) c(a1, a2, a3)

# プログラム更新しました # (1) 単回帰の概要 # 表のRのプログラム # install.packages("Ecdat") # install.packages("gt") # パッケージの読み込み library(Ecdat) library(gt) # データの準備 N <- 1:30 Ice <- cbind(N, Icecream) head(Ice) # 必要な列の選択 Ice0 <- Ice[, c("N", "cons", "temp")] # データの分割 Ice1 <- Ice0[1:10, ] Ice2 <- Ice0[11:20, ] Ice3 <- Ice0[21:30, ] # データフレームの結合 tb <- data.frame(Ice1, Ice2, Ice3) # gt() を使ってきれいな表を作成 tb %>% gt() %>% tab_header( title = "Package of Icecreame", subtitle = "cons vs temp" ) %>% tab_source_note( source_note = md("1951年から1953年にアメリカで30回行われたアイスクリームの消費量 (消費量 vs 気温)") ) library(Ecdat) # Rで回帰分析 # install.packages("Ecdat") library(Ecdat) cons<-Icecream$cons ; cons temp<-Icecream$temp ; temp # 変数 x, yにデータを格納（コピペして汎用化するため） x <- temp 　　　　 y <- cons # きれいな図を描くための設定 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) 　 plot(x, y, xlab = "気温(temp)", ylab = "消費量(cons)", cex.lab = 2, cex.axis = 1.5, lwd = 4) # 散布図の描画 # 回帰分析表の出力 S_reg <- lm(y~x) summary(S_reg ) # 回帰直線の描画 abline(S_reg, lwd=4, lty=2, col="red") # 平均値の描画 mx = mean(x) my =mean(y) abline(v=mx, lwd=2) abline(h=my, lwd=2)

# プログラム更新しました # (1) 単回帰の概要 # 表のRのプログラム # install.packages("Ecdat") # install.packages("gt") # パッケージの読み込み library(Ecdat) library(gt) # データの準備 N <- 1:30 Ice <- cbind(N, Icecream) head(Ice) # 必要な列の選択 Ice0 <- Ice[, c("N", "cons", "temp")] # データの分割 Ice1 <- Ice0[1:10, ] Ice2 <- Ice0[11:20, ] Ice3 <- Ice0[21:30, ] # データフレームの結合 tb <- data.frame(Ice1, Ice2, Ice3) # gt() を使ってきれいな表を作成 tb %>% gt() %>% tab_header( title = "Package of Icecreame", subtitle = "cons vs temp" ) %>% tab_source_note( source_note = md("1951年から1953年にアメリカで30回行われたアイスクリームの消費量 (消費量 vs 気温)") ) library(Ecdat) # Rで回帰分析 # install.packages("Ecdat") library(Ecdat) cons<-Icecream$cons ; cons temp<-Icecream$temp ; temp # 変数 x, yにデータを格納（コピペして汎用化するため） x <- temp 　　　　 y <- cons # きれいな図を描くための設定 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) 　 plot(x, y, xlab = "気温(temp)", ylab = "消費量(cons)", cex.lab = 2, cex.axis = 1.5, lwd = 4) # 散布図の描画 # 回帰分析表の出力 S_reg <- lm(y~x) summary(S_reg ) # 回帰直線の描画 abline(S_reg, lwd=4, lty=2, col="red") # 平均値の描画 mx = mean(x) my =mean(y) abline(v=mx, lwd=2) abline(h=my, lwd=2)

# プログラム更新しました # (2) 決定係数 # install.packages("Ecdat") library(Ecdat) x<-Icecream$temp y<-Icecream$cons # 回帰分析表の出力 S_reg <- lm(y~x) summary(S_reg)

# プログラム更新しました # (3) 分散分析 library(Ecdat) y<-Icecream$cons x<-Icecream$temp S_reg <- lm(y~x) summary(S_reg ) R2 <- cor(x,y)^2 n<- length(x) (F0 <- R2/(1-R2)*(n-2)) (Fq<- qf(1-0.001, 1, n-2)) (p値<-pf(F0, 1, (n-2), lower.tail=FALSE)) # 例題で示した図 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) # packages <-c("Ecdat", "tidyverse") # install.packages(packages) library(Ecdat) y<-Icecream$cons x<-Icecream$temp R2 <- cor(x, y)^2 # 決定係数 n<- length(x)　# サンプルサイズ # F分布 f1=1 ; f2=28　　# 自由度 a=0.001 xmax=50 ymax=0.002 curve(df(x,1,28),0,xmax, xlab="F", ylab="df(F,1,28)", ylim=c(0, ymax),cex.lab = 2, cex.axis = 1.5, lwd = 4) # F値０～50 # F値 (F <- R2/(1-R2)*(n-2)) # 検定統計量 (up= qf(1-a, f1, f2)) # 棄却域境界の分位点 F(1,n-2)(0.1%有意水準) abline(v=up, lty=2, col="red", cex.lab = 2, cex.axis = 1.5, lwd = 4) abline(v=F, lty=2, col="blue",cex.lab = 2, cex.axis = 1.5, lwd = 4) # 色を塗る　 xvals <- seq(up,xmax, length=2000) dvals <- df(xvals,f1,f2) polygon(c(xvals,rev(xvals)), c(rep(0,2000),rev(dvals)), col="skyblue")

# プログラム更新しました # Retrieve Cohen's suggested effect sizes library("pwr") cohen.ES(test = "f2", size = "small") cohen.ES(test = "f2", size = "medium") cohen.ES(test = "f2", size = "large") f <- function(x){ a <- sqrt(x/(1+x)) return(a) } a1=f(0.02) ; a2=f(0.15) ; a3=f(0.35) a1<-paste0("小=", round(a1, 2)) a2<-paste0("中=", round(a2, 2)) a3<-paste0("大=", round(a3, 2)) c(a1, a2, a3) # Power=0.8 サンプルサイズ # install.packages("pwr") library (pwr) v= round(pwr.f2.test(u=1, f2=0.02, sig.level=0.05, power=0.8)$v, 1) paste0 ("v=" ,v) # Power=0.9 # サンプルサイズ library(pwr) a <- 0.05 po <- 0.9 f <- function(x) { a <- ceiling(pwr.f2.test(u = 1, f2 = x, sig.level = a, power = po)$v) + 2 return(a) } a1 <- f(0.02) a2 <- f(0.15) a3 <- f(0.35) a1 <- paste0("小(0.02)=", a1) a2 <- paste0("中(0.15)=", a2) a3 <- paste0("大(0.35)=", a3) c(a1, a2, a3)

◎Rのプログラム # 1) 最尤推定法 data<-c(0,1,1) logLH<-function(x,q) { sum(x)*log(q)+(length(x)-sum(x))*log(1-q) } par(mar = c(6, 7, 5, 2)) plot(seq(0,1,0.01), logLH(data,seq(0,1,0.01)) , type="l", col="blue", main = "対数尤度関数の曲線" , xlab ="θ" , ylab = "対数尤度", lwd = 3, cex = 2, cex.main = 2, cex.lab = 2, cex.axis = 2 ) opt <- optimize(function(q) logLH(data,q),c(0,1),maximum = TRUE) print(str(opt)) abline(h=opt$objective, lty=2, col="blue") abline(v=opt$maximum,lty=2, col="blue") # 2) # Rによる最尤推定値の計算 data<-c(0,1,1) L_Be<-function(x,p) { p^sum(x)*(1-p)^(length(x)-sum(x)) } plot(seq(0,1,0.01),L_Be(data,seq(0,1,0.01)),type="l",col="blue") opt <- optimize(function(p) L_Be(data,p),c(0,1),maximum = TRUE) abline(h =opt$objective, lty=2, col="blue") abline(v =opt$maximum, lty=2, col="blue") print(opt) # ベータ分布＜ par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) a <- 50 b <- 50 x <- seq(0.01, 1.0, len = 500) 　　　　 y <- dbeta(x,a,b) 　　　　 plot(x, y, type = "l",col=a+b,ylim=c(0,10), ylab ="PDF", cex.lab = 1.5, cex.axis = 1.5 , lwd = 4)　　 m = a/(a+b) 　　　　　　　　　　　 var = a*b/((a+b)^2*(a+b+1))　　 sdd = sqrt(var)*100　　　　　　　 abline(v=m) legend("topright", legend = c(paste("Mean = ", round(mean(x), digits=2)) , paste("SD = ", round(sdd, digits=2) , "%")), cex = 2, bty = "n") # ベイズ更新の可視化のプログラム< data<-c(0,1,1) L_Be<-function(x,p) { p^sum(x)*(1-p)^(length(x)-sum(x)) } plot(seq(0,1,0.01),L_Be(data,seq(0,1,0.01)),type="l", xlab="", ylab="", ylim=c(0,8), xaxt="n", yaxt="n", col="blue") prior_beta<-function(p,a,b) dbeta(p,a,b)　　 a = 50 ; b = 50 beta<-function(x,p) prior_beta(p,a,b) 　 par(new=T) plot(seq(0,1,0.01),beta(data,seq(0,1,0.01)), type="l", xlab="", ylab="",ylim=c(0,8), xaxt="n", yaxt="n", col="red")　　 joint<-function(x,p) L_Be(x,p)*prior_beta(p,a,b)　　 par(new=T) plot(seq(0,1,0.01),joint(data,seq(0,1,0.01)),　　　　 type="l", ylim=c(0,8),xlab="", ylab="") optimize(function(p) L_Be(data,p),c(0,1),maximum = TRUE) optimize(function(p) L_Be(data,p)*prior_beta(p,a,b),c(0,1),maximum = TRUE) # 4 なぜα＝β＝ 50 ? a <- 50 　　　　　　　　　　　　　　 b <- 50 x <- seq(0.01, 1, len = 500) y <- dbeta(x,a,b) 　　　　 # グラフの表示 plot(x, y, type = "l",col="red",ylim=c(0,10))　　 m = a/(a+b)*100　　　　　　 var = a*b/((a+b)^2*(a+b+1))　 sdpt = sqrt(var)*100　　　　 paste("mean =", round(m, digits=2), "%") paste("sdv =", round(sdpt, digits=2) , "%") abline(v=m) install.packages("nleqslv") library(nleqslv) E=0.5　 V=0.05^2　 th<-function(x, E_th, V_th){ a=x[1] b=x[2] c(E_th-a/(a+b), V_th-(a*b)/(((a+b)^2)*(a+b+1))) } f<-function(x){ th(x, E, V) } an<-nleqslv(c(1,1), f) paste("α =" , round(an$x[1],digits=2), "β =", round(an$x[2], digits=2))

ハズレの扉を開けた後、選んだ扉を変えたら当たるかハズレるかの隔離では無いのですか•́ω•̀)? そこの所がイマイチ分かりません( ˘•ω•˘ ).｡oஇ

◎Rのプログラム # 1. ライブラリの読み込み library(kernlab) # 2. 乱数のシードを設定 set.seed(123) # 3. ランダムデータを生成 x = rnorm(1000) y = rnorm(1000) data = data.frame(type = 1, x, y) # 4. One-Class SVMモデルを作成 one_class_SVM = ksvm(type ~ ., data = data, type = "one-svc", kernel = "rbfdot", kpar = list(sigma = 0.1), cross = 10, nu = 0.1) # 5. 正常値と外れ値を予測 pre = predict(one_class_SVM) # 6. 予測結果の変換 c.pre = ifelse(pre == TRUE, 1, 2) # 7. 予測結果の結合 data.result = cbind(data, c.pre) # 8. 結果のプロット：正常値は青、外れ値は赤、外れ値は×で表示 plot(data.result[, 2:3], pch = ifelse(data.result$c.pre == 1, 21, 4), bg = ifelse(data.result$c.pre == 1, "blue", "red"), col = ifelse(data.result$c.pre == 1, "blue", "red"), cex = ifelse(data.result$c.pre == 1, 1, 1.5), lwd = ifelse(data.result$c.pre == 1, 1, 2), xlab = "X 値", ylab = "Y 値", main = "One-Class SVMによる異常検知", col.main = "blue", col.lab = "blue", col.axis = "blue") legend("topright", legend = c("正常値", "外れ値"), col = c("blue", "red"), pch = c(21, 4), pt.bg = c("blue", "red"), pt.cex = c(1, 1.5), lwd = c(1, 2))

◎Rのプログラム # その２ x<- c(30,36,47,48,50,52,55,55,55,60) y<- c(16,29,54,55,33,56,48,57,62,72) z<- c(48,69,108,114,67,124,91,106,116,132) d<-data.frame(x,y,z) View(d) d.lm<-lm(z~.,data=d) summary(d.lm) library(psych) pairs.panels(d) # その３ # AIC前の分析 y <- c(172, 156, 158, 168, 180, 170, 165, 175, 169, 155) # 身長 cm x1 <- c(75, 55, 60, 65, 80, 70, 60, 75, 70, 50) # 体重 kg x2 <- c(111, 83, 89, 98, 119, 104, 89, 113, 104, 74) # ウエスト cm # データフレームの作成 data <- data.frame("身長 cm" = y, "体重 kg" = x1, "ウエスト cm" = x2) # 1) 相関係数を求める coef <- round(cor(data), digits = 2) round(coef,digits = 2) # 変数間の相関をPlot pairs(data, pch=21, bg="red", cex=2) # 2) 回帰分析を実行 lm1 <- lm(y ~ x1 + x2, data=data) # 回帰分析結果を表示 round(coefficients(lm1), 2) # 回帰係数と切片の算出 summary(lm1) # 回帰分析の実行 # (参考)分析のもとになったデータと予測値、残差を一覧にする exp <- predict(lm1) # 元データに対する予測値 res <- residuals(lm1) # データと予測値の残差 view_lm1 <- data.frame(data[1], exp, res) # データフレームにまとめる print(round(view_lm1, digits =2 )) # 3)ステップワイズ(説明変数を減らしてAICを求める） step<-step(lm1) 　 # AIC後の分析 y <- c(172, 156, 158, 168, 180, 170, 165, 175, 169, 155) # 身長 cm x1 <- c(75, 55, 60, 65, 80, 70, 60, 75, 70, 50) # 体重 kg x2 <- c(111, 83, 89, 98, 119, 104, 89, 113, 104, 74) # ウエスト cm # AIC後の新たなモデルで回帰分析を実行 lm2 <- lm(y ~ x2, data=data) # 回帰分析結果を表示 round(coefficients(lm2), 2) # 回帰係数と切片の算出 summary(lm2) # 回帰分析の実行

◎Rのプログラム # その１# 自己共分散を計算する (acf(Nile, type = "covariance", lag.max = 5, plot = FALSE)) # 自己相関係数を計算する (acf(Nile, lag.max = 5, plot = FALSE)) plot(Nile) acf(Nile) acf(Nile, lag.max = 10) # その２ # プロットの設定 par(cex.lab = 1.5, cex.axis = 2.0, lwd = 3, cex.main = 4, oma = c(0, 1, 0, 0)) # Nileデータのプロット plot(Nile, main = "Nile Data", xlab = "Year", ylab = "Flow", type = "l") # 新しい描画領域を作成 dev.new() # グラフ描画領域を1行に2つ並べる par(mfrow = c(1, 2)) # プロットの設定 par(cex.lab = 1.5, cex.axis = 2.0, lwd = 4, cex.main = 3) # Nileデータの自己相関係数のコレログラムを描く acf(Nile, main = "ACF of Nile Data", xlab = "Lag", ylab = "ACF") # Nileデータの偏自己相関係数のコレログラムを描く acf(Nile,type="p", main = "PACF of Nile Data", xlab = "Lag", ylab = "PACF")

◎Rのプログラム # 5 サンプルサイズの設計 # install.packages("pwr") # はじめて使用する場合、先頭の # を消してください。 weight <-c(41, 36, 52, 38, 45, 42) score <-c(72, 81, 90, 72, 85, 77) r0=cor(weight, score) ; r0 # 相関係数 # 散布図 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) plot(weight, score, xlab = "体重(kg)", ylab = "点数(点)", cex.lab = 2, cex.axis = 1.5, lwd = 4) library("pwr") a=0.05 # 有意水準 # 検出力 pwr.r.test(n= 6, r=r0, sig.level=a, alternative=c("two.sided")) # サンプルサイズの設計 pwr.r.test(power =0.8, r=r0, sig.level=a, alternative=c("two.sided"))

◎Rのプログラム # 4 相関係数の検定 # -----簡単な方法 weight <-c(41, 36, 52, 38, 45, 42) score <-c(72, 81, 90, 72, 85, 77) # 散布図 　par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) 　plot(weight, score, xlab = "体重(kg)", ylab = "点数(点)", cex.lab = 2, cex.axis = 1.5, lwd = 4) cor(weight, score) # 相関係数 r cor.test(weight, score)　　　# 相関係数の検定 # -----詳細な方法 weight <-c(41, 36, 52, 38, 45, 42) score <-c(72, 81, 90, 72, 85, 77) # 検定統計量をもとめる n=length(weight) r=cor(weight,score,method="pearson") r # 相関係数 t=r*sqrt(n-2)/(sqrt(1-r^2)) # 相関係数 t 検定統計量 t # 相関係数 t 検定統計量 # 散布図とt分布を描画 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) # 自由度をもとめる(n-2) df=n-2 curve(dt(x,df),-5,5, cex.lab = 2, cex.axis = 1.5, lwd = 4) # t 分布曲線を書く # 図上に表記 qt(0.025,df) # 下方2.5%点の導出 qt(0.975,df) # 上方2.5%点の導出 abline(v=qt(0.025,df))　# 下方2.5%点の線 abline(v=qt(0.975,df)) # 上方2.5%点の線 abline(v=t,col=2, lwd = 4) # t 検定統計量の線 2*(1-pt(t, df)) # p値をもとめる　2*pt(t,df,lower.tail =FALSE)　でもよい

◎Rのプログラム #2 相関係数 library("corrr") # ------変数x (身長) cm xvar <- c(161.6, 153.9, 161.2, 172.0, 158.7, 163.3, 155.1, 159.8, 163.5, 147.9) # ------変数y (体重) kg yvar <- c(52.5, 46.7, 49.0, 58.5, 55.2, 48.3, 53.2, 47.0, 60.2, 41.7) # 散布図の描画 par(mar = c(5.5, 6.0, 4.1, 2)) par(mgp = c(4, 1.2, 0)) 　plot(xvar, yvar, xlab = "身長 cm", ylab = "体重 kg", cex.lab = 2, cex.axis = 1.5, lwd = 4) # ------関数を定義(2変数の標本共分散) # 2変数の標本共分散 (分母がn） sample_covariance <- function(x, y) { var(x, y) * (length(x)-1)/length(x) } sample_covariance(xvar, yvar) # 2変数の不偏共分散(分母がn-1) cov(xvar, yvar) # 2変数の相関係数 cor (xvar, yvar) # 参考 　　# (不偏）分散共分散行列で求める 　　sample_data <-data.frame(xvar, yvar) cov(sample_data) var(sample_data) # 相関行列を求める 　 cor(sample_data)

◎Rのプログラム #1 共分散と相関係数 # gtパッケージがインストールされていない場合、次の行のコメントを外してください # install.packages("gt") # 初回だけコメントを外す # 必要なライブラリを読み込む library(gt) # airqualityデータセットの最初の数行を表示 head(airquality) # airqualityデータセットの最初の数行を取得 airquality_head <- head(airquality) # gtパッケージで表を作成し、カスタマイズ airquality_table <- gt(data = airquality_head) %>% tab_header( title = "大気質データセット", subtitle = "最初の6行" ) %>% tab_style( style = list( cell_fill(color = "lightblue"), cell_text(weight = "bold") ), locations = cells_column_labels( columns = everything() ) ) %>% tab_options( table.border.top.color = "black", table.border.bottom.color = "black", table.border.top.width = px(3), table.border.bottom.width = px(3), row.striping.include_table_body = TRUE ) # 表を表示 print(airquality_table) # グラフを横に並べて表示するためにプロット領域を設定 par(mfrow = c(1, 3)) # プロットのための変数を割り当てる x <- airquality$Ozone # オゾン y <- airquality$Wind # 風速 # プロットの余白とラベルの位置を設定 par(mar = c(5.5, 6.0, 4.1, 2)) # 余白の広さを行数で指定.下，左，上，右の順. par(mgp = c(4, 1.2, 0)) # 余白の使い方.説明，ラベル，軸の位置を行で指定. # プロットを作成 plot(x, y, xlab = "Ozone", ylab = "Wind", cex.lab = 2, cex.axis = 1.5, lwd = 4) ######################## # 身長の単位が cmの場合のプロット xvar_cm <- c(161.6, 153.9, 161.2, 172.0, 158.7, 163.3, 155.1, 159.8, 163.5, 147.9) yvar <- c(52.5, 46.7, 49.0, 58.5, 55.2, 48.3, 53.2, 47.0, 60.2, 41.7) par(mar = c(5.5, 6.0, 4.1, 2)) # 余白の広さを行数で指定.下，左，上，右の順. par(mgp = c(4, 1.2, 0)) # 余白の使い方.説明，ラベル，軸の位置を行で指定. plot(xvar_cm, yvar, xlab = "身長 cm", ylab = "体重 kg", cex.lab = 2, cex.axis = 1.5, lwd = 4) # 身長の単位が mの場合のプロット xvar_m <- 0.01 * c(161.6, 153.9, 161.2, 172.0, 158.7, 163.3, 155.1, 159.8, 163.5, 147.9) yvar <- c(52.5, 46.7, 49.0, 58.5, 55.2, 48.3, 53.2, 47.0, 60.2, 41.7) par(mar = c(5.5, 6.0, 4.1, 2)) # 余白の広さを行数で指定.下，左，上，右の順. par(mgp = c(4, 1.2, 0)) # 余白の使い方.説明，ラベル，軸の位置を行で指定. plot(xvar_m, yvar, xlab = "身長 m", ylab = "体重 kg", cex.lab = 2, cex.axis = 1.5, lwd = 4) # 標本共分散（分母がn） sample_covariance <- function(x, y) { var(x, y) * (length(x) - 1) / length(x) } sample_covariance(xvar_cm, yvar) # 2変数の不偏共分散（分母がn-1） cov(xvar_cm, yvar) # 分散共分散行列（不偏） mx_cm <- matrix(c(xvar_cm, yvar), length(xvar_cm), 2) mx_cm cov(mx_cm) # 標本共分散（分母がn） sample_covariance(xvar_m, yvar) # 2変数の不偏共分散（分母がn-1） cov(xvar_m, yvar) # 分散共分散行列（不偏） mx_m <- matrix(c(xvar_m, yvar), length(xvar_m), 2) mx_m cov(mx_m)

00'57"" 総和記号のｋ⇒i です。

違った答えって、ルール変えたら確率変わるの当たり前やん・・・

後半、04:38以降、音声セトラブルのため、一部動画をCutしました。もうしわけありません。

◎Rのプログラム（標本分散と不偏分散） install.packages("gt") #パッケージをよみこむ(一度行えばよい） X <- c(184.2, 177.7, 168.0, 165.3, 159.1, 176.4, 176.0, 170.0, 177.3, 174.5, 164.6, 174.4, 174.8, 160.8, 162.1, 167.0, 167.3, 172.8, 168.1, 173.5) #----------きれいな表を書くgt()パッケージの読み込み library(gt) #----------hist()関数の戻り値を利用して度数分布を書く準備を行う h <- hist(X) h # hist関数の戻り値の表示 n = length(X) #サンプルサイズ rela_freq <- round(h$counts/n, digits = 2) #相対度数の算出 rank_n <- length(h$counts) # 階級の数 #----------階級の設定 class_names <- NULL # 階級の名前格納用 for(i in 1:rank_n) { class_names[i] <- paste(h$breaks[i], "～", h$breaks[i+1]) } #データフレームにまとめる frequency_table <- data.frame( 階級=class_names, 階級値=h$mids, 度数=h$counts, 相対度数= rela_freq, 累積度数=cumsum(h$counts), 累積相対度数= cumsum(rela_freq), 密度=h$density ) #データフレームを表示する frequency_table %>% gt() %>% tab_header(title = "度数分布表")

◎Rのプログラム # データ点の生成と散布図の描画 set.seed(123) x <- cbind(x = rnorm(10, 0, 1), y = rnorm(10, 0, 1)) x <- rbind(x, cbind(x = rnorm(10, 5, 1), y = rnorm(10, 5, 1))) x <- rbind(x, cbind(x = rnorm(10, 5, 1), y = rnorm(10, -5, 1))) plot(x, pch = 19, col = "blue", xlab = "X", ylab = "Y") # クラスタ数(k)の選択 k <- 3 # クラスタリングを実行 kmeans_result <- kmeans(x, centers = k) # k-平均法によるクラスタリング結果をプロット plot(x, col = kmeans_result$cluster, pch = 19) points(kmeans_result$centers, col = 1:k, pch = 8, cex = 2)

◎Rのプログラム(標準偏差UとUs) library(MASS) set.seed(111) # 母集団の生成 xx <- rnorm(250000) # サンプルサイズの設定 s_size <- 4 # 結果を格納するための空のベクトルを作成 U <- c() Us <- c() # サンプリングと統計量の計算 for (i in 1:2500) { s_xx <- sample(xx, size = s_size) # 標準偏差の計算 U <- append(U, sqrt(sum((s_xx - mean(s_xx))^2) / (s_size - 1))) Us <- append(Us, sqrt((s_size - 1) / 2) * gamma((s_size - 1)/2) / gamma(s_size/2) * sqrt(sum((s_xx - mean(s_xx))^2) / (s_size - 1))) } # 標準偏差Uのヒストグラム truehist(U, main=" U Square Root of Unbiased Sample Variance ", xlab="") abline(v=mean(U), lwd=5, col="blue") abline(v=1, lty=2, col="magenta", lwd=2) legend("topright", col=c("blue", "magenta"), lwd=c(5,2), lty=c(1,2), legend=c("Expected value of samples", "Population S.D.")) dev.new() # 標準偏差Usのヒストグラム truehist(Us, main="　Us Unbiased Standard Deviation ", xlab="") abline(v=mean(Us), lwd=5, col="blue") abline(v=1, lty=2, col="magenta", lwd=2) legend("topright", col=c("blue", "magenta"), lwd=c(5,2), lty=c(1,2), legend=c("Expected value of samples", "Population S.D."))

◎Rのプログラム（標本分散と不偏分散） library(MASS) set.seed(111) # 母集団の生成 xx <- rnorm(250000) # サンプルサイズの設定 s_size <- 4 # 結果を格納するための空のベクトルを作成 S2 <- c() U2 <- c() # サンプリングと統計量の計算 for (i in 1:2500) { s_xx <- sample(xx, size = s_size) #標本分散と不偏分散の計算 S2 <- append(S2, sum((s_xx - mean(s_xx))^2) / s_size) U2 <- append(U2, var(s_xx)) } #ヒストグラムの描画 # Sample Variance のヒストグラム truehist(S2, main="S^2 Sample Variance", xlab="") abline(v=mean(S2), lwd=5, col="blue") abline(v=1, lty=2, col="magenta", lwd=2) legend("topright", col=c("blue", "magenta"), lwd=c(5,2), lty=c(1,2), legend=c("Expected value of samples", "Population S.D.")) dev.new() # Unbiased Sample Variance のヒストグラム truehist(U2, main="U^2 Unbiased Sample Variance", xlab="") abline(v=mean(U2), lwd=5, col="blue") abline(v=1, lty=2, col="magenta", lwd=2) legend("topright", col=c("blue", "magenta"), lwd=c(5,2), lty=c(1,2), legend=c("Expected value of samples", "Population S.D."))

◎Rプログラム # install.packages("multcomp") # 最初に1度だけ # multcompパッケージを読み込む library(multcomp) # データセット data <- data.frame(category = factor(c(rep(1, 3), rep(2, 3), rep(3, 3)), labels = c("control", "B", "C")), response = c(10.2, 9.8, 10.6, 8.1, 9.5, 10.5, 12.7, 12.5, 11.9)) # Boxplotを描画 boxplot(response ~ category, data = data, main = "Production", xlab = "Group", ylab = "Response") dev.new() # ANOVAを実行 model <- lm(response ~ category, data = data) (anova_result <- summary(aov(model))) # Tukey's Multiple Comparisons (tukey_result <- glht(model, linfct = mcp(category = "Tukey"))) # Tukey-Kramer法 (tukey_kramer_result <- TukeyHSD(aov(model))) # Dunnett's Multiple Comparisons Dunnett_result <- glht(model, linfct = mcp(category = "Dunnett")) print(summary(Dunnett_result)) # グラフ描画 plot(tukey_result, cex.axis = 0.8, las = 1) dev.new() plot(tukey_kramer_result, cex.axis = 0.8, las = 1) dev.new() plot(Dunnett_result, cex.axis = 0.8)

お詫び　ダネットのtが一つ足りませんでした！サムネイルは修正しました。

◎Rのプログラム # install.packages("multcomp") # 最初に1度だけ # multcompパッケージを読み込む library(multcomp) # データセット data <- data.frame(category = factor(c(rep(1, 3), rep(2, 3), rep(3, 3)), labels = c("control", "B", "C")), response = c(10.2, 9.8, 10.6, 8.1, 9.5, 10.5, 12.7, 12.5, 11.9)) # Boxplotを描画 boxplot(response ~ category, data = data, main = "Production", xlab = "Group", ylab = "Response") dev.new() # ANOVAを実行 model <- lm(response ~ category, data = data) (anova_result <- summary(aov(model))) # Tukey-Kramer法 (tukey_kramer_result <- TukeyHSD(aov(model))) # Dunnett's Multiple Comparisons dunnett_result <- glht(model, linfct = mcp(category = "Dunnett")) print(summary(dunnett_result)) # グラフ描画 plot(tukey_kramer_result, cex.axis = 0.8, las = 1) dev.new() plot(dunnett_result, cex.axis = 0.8)

訂正　相関の連続、対応ありは　ノンパラメトリック⇒パラメトリック

#Pythonでグラムシュミットの正規直交化法 iimport numpy as np # 例題のベクトル a1 = np.array([0, 1], dtype=np.float64) a2 = np.array([1, 1], dtype=np.float64) # Gram-Schmidtの正規直交化 def gram_schmidt(vectors): num_vectors = len(vectors) ortho_basis = np.zeros_like(vectors, dtype=np.float64) # Gram-Schmidtのアルゴリズムに基づいて直交化 for i in range(num_vectors): ortho_basis[i] = vectors[i] for j in range(i): # 直交成分を計算して減算 ortho_basis[i] -= np.dot(vectors[i], ortho_basis[j]) / np.dot(ortho_basis[j], ortho_basis[j]) * ortho_basis[j] # 正規化 for i in range(num_vectors): ortho_basis[i] /= np.linalg.norm(ortho_basis[i]) return ortho_basis vectors = np.array([a1, a2], dtype=np.float64) ortho_basis = gram_schmidt(vectors) u1 = ortho_basis[0] u2 = ortho_basis[1] print("u1:", u1) print("u2:", u2)

分かりやすかったよ

すみません。再生リストを1から順番に並べていただくことはできないでしょうか？

JASPすごいですよ！

◎Rのプログラム 1 # 1 Duncanデータセット # 必要なパッケージをインストール(未インストールならば） required_packages <- c("car", "Hmisc") install_packages <- installed.packages() for (pkg in required_packages) { if (!(pkg %in% install_packages)) { print(paste(pkg, "is not installed. Installing...")) install.packages(pkg) } }

後半の標準正規分布はは、正確にはσ２乗ではなく、σイコール１の２乗ですね。

サムネの漢字は「誤る」じゃなくて「謝る」？？😮

今日もわかりみMAXでした！

◎Rのプログラム 2 # 2 #Boxplotとヒストグラムの表示 # 必要なパッケージリスト required_packages <- c("car", "Hmisc", "ggplot2", "cowplot") # インストールされていないパッケージだけをインストール installed_packages <- installed.packages() to_install <- required_packages[!(required_packages %in% installed_packages$Package)] if (length(to_install) > 0) { install.packages(to_install) } # パッケージをロード library(car) library(Hmisc) library(ggplot2) library(cowplot) # Duncan データセットを読み込む data(Duncan) # ボックスプロットを描画 boxplot_plot <- ggplot(Duncan, aes(x = type, y = income)) + geom_boxplot() + labs(title = "Income by Type", x = "Type", y = "Income") + theme( plot.title = element_text(size = 25), # タイトルの文字サイズを調整 axis.text = element_text(size = 20), # 軸のテキストの文字サイズを調整 axis.title = element_text(size = 22) # 軸のタイトルの文字サイズを調整 ) # ヒストグラムを描画 hist_plots <- lapply(c("bc", "prof", "wc"), function(t) { ggplot(Duncan[Duncan$type == t, ], aes(x = income)) + geom_histogram(fill = "lightblue", color = "black", bins = 20) + labs(title = paste("Income Histogram for", t), x = "Income", y = "Frequency") + theme( plot.title = element_text(size = 25), # タイトルの文字サイズを調整 axis.text = element_text(size = 20), # 軸のテキストの文字サイズを調整 axis.title = element_text(size = 22) # 軸のタイトルの文字サイズを調整 ) }) # グラフを組み合わせて表示（1行3列） final_plot <- plot_grid(boxplot_plot, plotlist = hist_plots, ncol = 4) # グラフを表示 print(final_plot)

◎Rのプログラム 1 # 1 Duncanデータセット # 必要なパッケージをインストール(未インストールならば） required_packages <- c("car", "Hmisc") install_packages <- installed.packages() for (pkg in required_packages) { if (!(pkg %in% install_packages)) { print(paste(pkg, "is not installed. Installing...")) install.packages(pkg) } } # インストールされたパッケージをロード library(car) library(Hmisc) # Duncan データセットを読み込み、最初の10行を表示 data(Duncan) head(Duncan, 10) # type別に所得の統計情報を計算し、統計情報を表示 summary_stats <- aggregate(income ~ type, data = Duncan, FUN = function(x) { c( mean = mean(x), sd = sd(x), IQR = IQR(x), quantiles = quantile(x, c(0, 0.35, 0.5, 0.75, 1)), n = length(x) ) }) print(summary_stats) # ボックスプロットを描画 boxplot(income ~ type, data = Duncan, main = "Income by Type", xlab = "Type", ylab = "Income") # Kruskal-Wallis検定の結果からη²(イータ・スクエア）を手動で計算 group_means <- tapply(Duncan$income, Duncan$type, mean) grand_mean <- mean(Duncan$income) SSB <- sum((group_means - grand_mean)^2 * table(Duncan$type)) SST <- sum((Duncan$income - grand_mean)^2) eta_squared <- SSB / SST # 結果を表示 print(paste("Eta-squared (Effect Size):", format(round(eta_squared, 3), nsmall = 3))) # インストールされたパッケージをロード library(car) library(Hmisc) # Kruskal-Wallis検定を実行 kruskal_result <- kruskal.test(income ~ type, data = Duncan) # 結果の表示 print(kruskal_result)

◎Rプログラム # car パッケージがインストールされているか確認し、なければインストールする if (!requireNamespace("car", quietly = TRUE)) { install.packages("car") } # car パッケージをロード library(car) # Duncanデータセットの読み込み data("Duncan") # 1列目をtype、2列目をincomeにしたデータフレームを作成 df <- data.frame(type = Duncan$type, income = Duncan$income) # データフレームをCSVファイルとして保存 write.csv(df, "output_file.csv", row.names = FALSE)

◎Rプログラム　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　# ベクトル CD を作成 CD <- c(12, 24, 36, 48, 60, 72) # ベクトル CD を行列に変換 matrix_CD <- matrix(CD, nrow = 2, ncol = 3, byrow = TRUE) # 全体の平均を計算 total_mean <- mean(CD) # 行ごとの平均を計算 row_means <- apply(matrix_CD, 1, mean) # 列ごとの平均を計算 col_means <- apply(matrix_CD, 2, mean) # 結果を表示 print(matrix_CD) cat("全体の平均: ", total_mean, " ") cat("1行目の平均: ", row_means[1], " ") cat("2行目の平均: ", row_means[2], " ") cat("1列目の平均: ", col_means[1], " ") cat("2列目の平均: ", col_means[2], " ") cat("3列目の平均: ", col_means[3], " ")

◎Rプログラム library(car) library(KernSmooth) x <- Davis[, c("weight", "height")] h <- c(dpik(x$weight), dpik(x$height)) est <- bkde2D(x, bandwidth = h, gridsize = c(10^3, 10^3)) d <- list(x = est$x1, y = est$x2, z = est$fhat) image(d, col = terrain.colors(7), xlim = c(35, 110), ylim = c(145, 200)) contour(d, add = TRUE) dev.new() n <- nrow(x) K <- matrix(-1, n, n) prefac <- (2 * pi * h)^(-0.5) for(nn in 1:n){ xnn <- x[nn,] x_x1 <- x[,1] - as.numeric(xnn[1]) K1 <- prefac[1] * exp(-0.5 * (1 / h[1])^2 * (x_x1 * x_x1)) x_x2 <- x[,2] - as.numeric(xnn[2]) K2 <- prefac[2] * exp(-0.5 * (1 / h[2])^2 * (x_x2 * x_x2)) K[, nn] <- K1 * K2 } aa <- colSums(K) - diag(K) lowerLim <- 10^(-20) aa[aa < lowerLim] <- lowerLim a <- (-1) * log(aa / (n - 1)) plot(a, xlab = "sample ID", ylab = "anomaly score")

最高に可愛くてわかりやすい！

◎Rプログラム # サンプルサイズを指定 n_values <- c(3, 5, 10, 15) # グラフの描画 par(mfrow = c(1, 1), cex.axis = 1.5, cex.lab = 1.5, cex.main = 1.8) # 1x1のグリッドにプロット # 軸の範囲を設定 xlim <- c(-3, 3) ylim <- c(0, 0.4) # 新しいプロット領域を設定 plot(1, type = "n", xlab = "t", ylab = " f(t) ", xlim = xlim, ylim = ylim) # 各サンプルサイズに対してt分布をプロット for (i in 1:length(n_values)) { n <- n_values[i] # 自由度を計算 df <- n - 1 # t分布の確率密度関数の計算 x <- seq(-3, 3, length.out = 100) y <- dt(x, df) # グラフの描画 lines(x, y, col = rainbow(length(n_values))[i], lwd = 3) # Increase line width to 3 } # グラフにタイトルを追加 title(main = "t分布") # サンプルサイズごとに凡例を追加 legend("topright", legend = paste("φ ＝ n－1 ＝ ", n_values - 1), col = rainbow(length(n_values)), lwd = 3, cex = 1, bty = "n") # Reduce legend text size to 1

◎Rプログラム # install.packages("class") # 最初に１度だけインストール # 1)ランダムな正常データの生成 set.seed(123) normal_data <- matrix(rnorm(200), ncol = 2) # 2)異常データの生成 anomaly_data <- matrix(c(5, 5, -5, -5), ncol = 2) # 3) k近傍法による異常検知 library(class) k <- 3 # 近傍点の数 all_data <- rbind(normal_data, anomaly_data) labels <- c(rep("normal", nrow(normal_data)), rep("anomaly", nrow(anomaly_data))) knn_result <- knn(all_data, normal_data, labels, k = k) # 4)結果のプロット plot(all_data, col = ifelse(knn_result == "normal", "blue", "red"), pch = 19, main = "k-NN Anomaly Detection") points(anomaly_data, col = "red", pch = 19) legend("topright", legend = c("Normal", "Anomaly"), col = c("blue", "red"), pch = 19)