[艦これで学ぶ統計学] 第7回 非線形回帰分析

第4回から続く回帰分析の解説の続きです。

前回重回帰分析を行いましたが、搭載スロット数が2, 3のときに推定式がずれているように見受けられました。そこでよりよい推定式を考えてみます。ここで搭載スロットごとに個別に夜偵の判定を行うという仮説を立ててみます。このとき、すべての夜偵が不発する確率は以下のようになります。

夜偵不発率 = スロット当たりの夜偵不発率 ^ 搭載スロット数

よって、夜偵発動率は以下のようになります。

1 - 夜偵発動率 = (1 - スロットあたりの夜偵発動率) ^ 搭載スロット数

夜偵発動率 = 1 - (1 - スロットあたりの夜偵発動率) ^ 搭載スロット数

この仮説をもとに推定式を立ててみます。

夜偵発動率 = 1 - (1 - (a + b × √レベル)) ^ 搭載スロット数

この推定式をもとに回帰分析を行ってみます。このような複雑な推定式の回帰分析を非線形回帰分析といいます。Rで非線形回帰分析を行うには、Rに以下のように入力します。

data <- read.csv(fileEncoding = "UTF-8", "C:\\Users\\Master\\Downloads\\艦これで学ぶ統計学 - 第6回 重回帰分析.csv")

result <- nls(発動可否 ~ 1 - (1 - (a + b * sqrt(レベル))) ^ 搭載スロット数, data, list(a = 0, b = 0))

summary(result)

result <- nls(発動可否 ~ 1 - (1 - (a + b * sqrt(レベル))) ^ 搭載スロット数, data, list(a = 0, b = 0))とすることで非線形回帰分析を行うことができます。今までの回帰分析と違い、推定式のa, bも式の中に書く必要があります。

この結果から推定式は以下のように求められます。

夜偵発動率 = 1 - (1 - (-0.0382966±0.0059711) + (0.0716064±0.0008197) × √レベル) ^ 搭載スロット数

この推定式のグラフを書いてみます。前回の重回帰分析と同様の方法で書くことができます。Rに以下のように入力します。

result.x <- 1:150

result.newdata.1 <- data.frame(レベル = result.x, 搭載スロット数 = 1)

result.y.1 <- predict(result, result.newdata.1)

plot(result.x, result.y.1, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab = "夜偵発動率", type = "l", lty = 1, col = 1)

par(new = TRUE)

result.newdata.2 <- data.frame(レベル = result.x, 搭載スロット数 = 2)

result.y.2 <- predict(result, result.newdata.2)

plot(result.x, result.y.2, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab = "夜偵発動率", type = "l", lty = 2, col = 2)

par(new = TRUE)

result.newdata.3 <- data.frame(レベル = result.x, 搭載スロット数 = 3)

result.y.3 <- predict(result, result.newdata.3)

plot(result.x, result.y.3, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab = "夜偵発動率", type = "l", lty = 3, col = 3)

par(new = TRUE)

legend("bottomright", legend = c("搭載スロット数1", "搭載スロット数2", "搭載スロット数3"), lty = c(1, 2, 3), col = c(1, 2, 3))

error.bar <- function (x, y, ...) {

    x.mean <- mean(x)

    y.mean <- mean(y)

    y.se <- sqrt(var(y) / length(y))

    points(x.mean, y.mean, ...);

    arrows(x.mean, y.mean - y.se, x.mean, y.mean + y.se, angle = 90, code = 3, length = 0.05, ...)

}

data.1.1 <- subset(data, 搭載スロット数 == 1 & レベル == 1)

error.bar(data.1.1$レベル, data.1.1$発動可否, col = 1)

data.1.50.55 <- subset(data, 搭載スロット数 == 1 & 50 <= レベル & レベル <= 55)

error.bar(data.1.50.55$レベル, data.1.50.55$発動可否, col = 1)

data.1.99 <- subset(data, 搭載スロット数 == 1 & レベル == 99)

error.bar(data.1.99$レベル, data.1.99$発動可否, col = 1)

data.1.116.125 <- subset(data, 搭載スロット数 == 1 & 116 <= レベル & レベル <= 125)

error.bar(data.1.116.125$レベル, data.1.116.125$発動可否, col = 1)

data.1.150 <- subset(data, 搭載スロット数 == 1 & レベル == 150)

error.bar(data.1.150$レベル, data.1.150$発動可否, col = 1)

data.2.1 <- subset(data, 搭載スロット数 == 2 & レベル == 1)

error.bar(data.2.1$レベル, data.2.1$発動可否, col = 2)

data.2.99 <- subset(data, 搭載スロット数 == 2 & レベル == 99)

error.bar(data.2.99$レベル, data.2.99$発動可否, col = 2)

data.2.150 <- subset(data, 搭載スロット数 == 2 & レベル == 150)

error.bar(data.2.150$レベル, data.2.150$発動可否, col = 2)

data.3.1 <- subset(data, 搭載スロット数 == 3 & レベル == 1)

error.bar(data.3.1$レベル, data.3.1$発動可否, col = 3)

搭載スロット数が2, 3のときもあてはまっているようにみえます。

AICも計算してみます。Rに以下のように入力します。

AIC(result)

前回計算したAICよりも小さくなっているため当てはまりがよくなっていることがわかります。

第4回から続いた回帰分析の解説もここまでになります。皆さんの艦これプレイの助けになれば幸いです。

最後に今回使用したRのコードを載せておきます。

data <- read.csv(fileEncoding = "UTF-8", "C:\\Users\\Master\\Downloads\\艦これで学ぶ統計学 - 第6回 重回帰分析.csv")

result <- nls(発動可否 ~ 1 - (1 - (a + b * sqrt(レベル))) ^ 搭載スロット数, data, list(a = 0, b = 0))

summary(result)

result.x <- 1:150

result.newdata.1 <- data.frame(レベル = result.x, 搭載スロット数 = 1)

result.y.1 <- predict(result, result.newdata.1)

plot(result.x, result.y.1, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab = "夜偵発動率", type = "l", lty = 1, col = 1)

par(new = TRUE)

result.newdata.2 <- data.frame(レベル = result.x, 搭載スロット数 = 2)

result.y.2 <- predict(result, result.newdata.2)

plot(result.x, result.y.2, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab = "夜偵発動率", type = "l", lty = 2, col = 2)

par(new = TRUE)

result.newdata.3 <- data.frame(レベル = result.x, 搭載スロット数 = 3)

result.y.3 <- predict(result, result.newdata.3)

plot(result.x, result.y.3, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab = "夜偵発動率", type = "l", lty = 3, col = 3)

par(new = TRUE)

legend("bottomright", legend = c("搭載スロット数1", "搭載スロット数2", "搭載スロット数3"), lty = c(1, 2, 3), col = c(1, 2, 3))

error.bar <- function (x, y, ...) {

    x.mean <- mean(x)

    y.mean <- mean(y)

    y.se <- sqrt(var(y) / length(y))

    points(x.mean, y.mean, ...);

    arrows(x.mean, y.mean - y.se, x.mean, y.mean + y.se, angle = 90, code = 3, length = 0.05, ...)

}

data.1.1 <- subset(data, 搭載スロット数 == 1 & レベル == 1)

error.bar(data.1.1$レベル, data.1.1$発動可否, col = 1)

data.1.50.55 <- subset(data, 搭載スロット数 == 1 & 50 <= レベル & レベル <= 55)

error.bar(data.1.50.55$レベル, data.1.50.55$発動可否, col = 1)

data.1.99 <- subset(data, 搭載スロット数 == 1 & レベル == 99)

error.bar(data.1.99$レベル, data.1.99$発動可否, col = 1)

data.1.116.125 <- subset(data, 搭載スロット数 == 1 & 116 <= レベル & レベル <= 125)

error.bar(data.1.116.125$レベル, data.1.116.125$発動可否, col = 1)

data.1.150 <- subset(data, 搭載スロット数 == 1 & レベル == 150)

error.bar(data.1.150$レベル, data.1.150$発動可否, col = 1)

data.2.1 <- subset(data, 搭載スロット数 == 2 & レベル == 1)

error.bar(data.2.1$レベル, data.2.1$発動可否, col = 2)

data.2.99 <- subset(data, 搭載スロット数 == 2 & レベル == 99)

error.bar(data.2.99$レベル, data.2.99$発動可否, col = 2)

data.2.150 <- subset(data, 搭載スロット数 == 2 & レベル == 150)

error.bar(data.2.150$レベル, data.2.150$発動可否, col = 2)

data.3.1 <- subset(data, 搭載スロット数 == 3 & レベル == 1)

error.bar(data.3.1$レベル, data.3.1$発動可否, col = 3)

AIC(result)