Tokyo.R #70 で発表しました
Tokyo.R #70 「ベイズ・統計モデリング回」(ちゃんと名前ついてた訳ではないけど)にて
初心者セッション2を担当しました。
うまく喋れた気が全くしないので、まだまだ要精進です。
資料はこちらになります。
https://www.slideshare.net/secret/bKnQdsOodauEUu
Rコードを資料に記載できなかったので、
こちらに記載しておきます。
※実データの部分は、諸般の事情で抜いてます。
あえてset.seed()は設定していないので、
何回も乱数生成→プロットを繰り返して、
「乱数生成って本当に乱数生成してんだな」ってところを感じていただければ。(何を言ってるんだ)
# パッケージのインストール
library(tidyverse)
library(extraDistr)
library(metRology)
# 標準正規分布 ----------
# 乱数生成
X <- rnorm(5000, mean = 0, sd = 1)
# データフレーム化
df <- data.frame(X = X)
# プロット
ggplot(df_unif, aes(x = X, y = ..density..)) +
geom_histogram(bins = 50, alpha = 0.8) +
stat_density(position='identity', geom='line',
colour = "#D02090") +
theme_bw() +
guides(fill = "none")
# 一様分布 ----------
# 乱数生成
X_unif <- runif(5000, min = 0, max = 1)
# データフレーム化
df_unif <- data.frame(X = X_unif)
# データの可視化
# プロット
ggplot(df_unif) +
geom_histogram(aes(x = X), bins = 50, alpha = 0.8) +
theme_bw() +
guides(fill = "none")
# ベルヌーイ分布 ----------
# 乱数生成
X_bern <- extraDistr::rbern(5000, prob = 0.5)
# データフレーム化
df_bern <- data.frame(X = X_bern)
# プロット
ggplot(df_bern) +
geom_histogram(aes(x = X), alpha = 0.8) +
theme_bw() +
guides(fill = "none")
# 二項分布 ----------
# 乱数生成
X_binom <- rbinom(5000, size = 100, prob = 0.5)
# データフレーム化
df_binom <- data.frame(X = X_binom)
# プロット
ggplot(df_binom) +
geom_histogram(aes(x = X), alpha = 0.8) +
theme_bw() +
guides(fill = "none")
# 正規分布 ----------
# 乱数生成
X_norm <- rnorm(5000,
mean = 8471, # 実データの平均
sd = 3765.519) # 実データの標準偏差
# データフレーム化
df_norm <- data.frame(X = X_norm)
# プロット
ggplot(df_norm, aes(x = X, y = ..density..)) +
geom_histogram(bins = 30, alpha = 0.8) +
stat_density(position='identity', geom='line',
colour = "#D02090") +
theme_bw() +
guides(fill = "none")
# Studentのt分布 ----------
# 乱数生成
X_t <- metRology::rt.scaled(5000, df = 5,
mean = 8471, # 実データの平均
sd = 3765.519) # 実データの標準偏差
# データフレーム化
df_t <- data.frame(X = X_t)
# プロット
ggplot(df_t, aes(x = X, y = ..density..)) +
geom_histogram(bins = 30, alpha = 0.8) +
stat_density(position='identity', geom='line',
colour = "#D02090") +
theme_bw() +
guides(fill = "none")
# コーシー分布 ----------
# 乱数生成
X_cauchy <- rcauchy(123,
location = 8471, # 位置パラメータ(実データの平均)
scale = 3) # 尺度パラメータ(適当)
# データフレーム化
df_cauchy <- data.frame(X = X_cauchy)
# プロット
dens <- density(df_cauchy$X)
ggplot(df_cauchy, aes(x = X)) +
geom_histogram(aes(y=..density..), bins = 30, alpha = 0.8) +
geom_density(colour = "#D02090") +
theme_bw() +
xlim(range(dens$x)) +
guides(fill = "none")
# 対数正規分布 ----------
# 実データの平均・標準偏差に合わせて乱数を生成する
sdlog <- sqrt(log((3765.519/8471)^2 + 1))
meanlog <- log(8471) - (sdlog^2) / 2
X_rlnorm <- rlnorm(5000, meanlog = meanlog, sdlog = sdlog)
# データフレーム化
df_rlnorm <- data.frame(X = X_rlnorm)
# プロット
ggplot(df_rlnorm, aes(x = X, y = ..density..)) +
geom_histogram(bins = 30, alpha = 0.8) +
stat_density(position='identity', geom='line',
colour = "#D02090") +
theme_bw() +
guides(fill = "none")
# ポアソン分布 ----------
# 乱数生成
X_pois <- rpois(5000,
lambda = 3.725) # 実データの平均値
# データフレーム化
df_pois <- data.frame(X = X_pois)
ggplot(df_pois, aes(x = X)) +
geom_histogram(bins = 30, alpha = 0.8) +
theme_bw() +
xlim(0, 30) +
guides(fill = "none")
# ゼロ過剰ポアソン分布 ----------
# 乱数生成
# ベルヌーイ分布で、0の多い0/1の乱数を作成
# rbinom()でsize=1を指定してもできます
# 170/1000は、実データの比率
X_bern <- rbinom(n = 5000, size = 1, prob = 170/1000)
# ポアソン分布に従う乱数生成
X_zero_pois <- c(rpois(2000, lambda = 3.725), X_bern[X_bern == 0])
# データフレーム化
df_zero_pois <- data.frame(X = X_zero_pois)
# プロット
ggplot(df_zero_pois, aes(x = X)) +
geom_histogram(bins = 10, alpha = 0.8) +
theme_bw() +
guides(fill = "none")
# 指数分布 ----------
# 乱数生成
X_exp <- rexp(5000, rate = 170/1000)
# データフレーム化
df_exp <- data.frame(X = X_exp)
# プロット
ggplot(df_exp, aes(x = X, y = ..density..)) +
geom_histogram(bins = 20, alpha = 0.8) +
stat_density(position='identity', geom='line',
colour = "#D02090") +
theme_bw() +
guides(fill = "none")
# ガンマ分布 -----------
# 乱数生成
shape = 3 # 3回実施する
rate = 30 / 7 # 1週間に1回 / 1ヶ月に3回
X_gamma <- rgamma(5000, shape = shape, rate = rate)
# データフレーム化
df_gamma <- data.frame(X = X_gamma)
# プロット
ggplot(df_gamma, aes(x = X, y = ..density..)) +
geom_histogram(bins = 20, alpha = 0.8) +
stat_density(position='identity', geom='line',
colour = "#D02090") +
theme_bw() +
guides(fill = "none")
# ある番組の視聴人数のシミュレーション ----------
N <- 7 # 1週間(7日)
X <- 1:N # 1周期のタイムポイントの数
Intercept <- 50 # 一番最初の値 (※適当)
gamma <- 0.7 # 指数関数の係数(0-1の値) (※適当)
r <- 10 # リフトする係数 (※適当)
# 1周期目
mu_1 = c(Intercept * gamma ^ X)
# 2周期目
mu_2 = c(mu_1[N] * r * gamma ^ X) # 1周期目の最後の値を2倍して、それが切片
# 3周期目
mu_3 = c(mu_2[N] * r * gamma ^ X)
# 4周期目
mu_4 = c(mu_3[N] * r * gamma ^ X)
# 4周期分をがっちゃんこ
mu <- c(mu_1, mu_2, mu_3, mu_4)
# 指数関数で"Xが上がるほど切片が減衰する"のデータをYに
Y <- rnorm(4 * N, mean = mu, sd = 0.5)
# プロット
X_label <- 1:28 # 1:7を4回くりかえしたので、X軸として1:28をうつ
goboten_udon <- data.frame(X = X_label, Y = mu) %>%
ggplot(mapping = aes(x = X, y = Y, colour = Y)) +
theme_bw() +
geom_point(size = 2.5) +
geom_line(mapping = aes(y = mu), size = 1.5, alpha = 0.5) +
scale_colour_gradient(low = "#00868B", high = "#4169E1") +
guides(colour = "none")
plot(goboten_udon)
「ベイズ・統計モデリング回」は、個人的にやりたかった回なのですが、
実際にやってみると、参加者の皆様が楽しみにしてくださったみたいで、
企画してよかったなぁ、と思いました。
いつもと違う内容にも関わらず、レベル高く発表してくださった
@kilometer00さん、@y__mattuさんには、感謝しかありません。。
本当にありがとうございます。
受付や懇親会準備等で色々手を動かしてくださった運営メンバーの皆様にも
本当に感謝です。
あと三銃士はすごかった。。神やった。。
お忙しい中、素晴らしい発表をしてくださり、本当にありがとうございました。
