からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【R】1次元ガウス分布の作図

はじめに

 機械学習で登場する確率分布について色々な角度から理解したいシリーズです。

 1次元ガウス分布の計算とグラフの作成をR言語で行います。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

1次元ガウス分布の作図

 1次元ガウス分布(Gaussian Distribution)の計算と作図を行います。

 利用するパッケージを読み込みます。

# 利用するパッケージ
library(tidyverse)
library(gganimate)

 分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。

定義式の確認

 まずは、1次元ガウス分布の定義式を確認します。

 1次元ガウス分布は、次の式で定義されます。

$$ \mathcal{N}(x | \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right) $$

 ここで、$\mu$は平均、$\sigma$は標準偏差、$\sigma^2$は分散です。
 確率変数の値$x$は実数となります。平均パラメータ$\mu$は実数、標準偏差パラメータ$\sigma$は0より大きい値$\sigma > 0$を満たす必要があります。

 この式の対数をとると、次の式になります。

$$ \log \mathcal{N}(x | \mu, \sigma^2) = - \frac{1}{2} \left\{ \log (2 \pi) + \log \sigma^2 + \frac{(x - \mu)^2}{\sigma^2} \right\} $$

 1次元ガウス分布の平均と分散は、次の式で計算できます。

$$ \begin{aligned} \mathbb{E}[x] &= \mu \\ \mathbb{V}[x] &= \sigma^2 \end{aligned} $$


 1次元ガウス分布は、実数$x$を生成することから、1次元ガウス分布の平均パラメータのパラメータの事前分布として利用されます。

確率密度の計算

 1次元ガウス分布に従う確率密度を計算する方法をいくつか確認します。

 パラメータを設定します。

# 平均を指定
mu <- 1

# 標準偏差を指定
sigma <- 2.5

# 確率変数の値を指定
x <- 1

 1次元ガウス分布の平均パラメータ$\mu$、標準偏差パラメータ$\sigma > 0$、確率変数がとり得る値$x$を指定します。設定した値に従う確率密度を計算します。

 まずは、定義式から確率密度を計算します。

# 定義式により確率密度を計算
C <- 1 / sqrt(2 * pi * sigma^2)
dens <- C * exp(-0.5 * (x - mu)^2 / sigma^2)
dens
## [1] 0.1595769

 1次元ガウス分布の定義式

$$ \begin{aligned} C_{\mathcal{N}} &= \frac{1}{\sqrt{2 \pi \sigma^2}} \\ \mathcal{N}(x | \mu, \sigma^2) &= C_{\mathcal{N}} \exp \left( - \frac{1}{2 \sigma^2} (x - \mu)^2 \right) \end{aligned} $$

で計算します。$C_{\mathcal{N}}$は1次元ガウス分布の正規化係数、$\pi$は円周率です。

 対数をとった定義式から計算します。

# 対数をとった定義式により確率密度を計算
log_C <- -0.5 * log(2 * pi) - log(sigma)
log_dens <- log_C - 0.5 * (x - mu)^2 / sigma^2
dens <- exp(log_dens)
dens; log_dens
## [1] 0.1595769
## [1] -1.835229

 対数をとった定義式

$$ \begin{aligned} \log C_{\mathcal{N}} &= - \frac{1}{2} \log (2 \pi) - \log \sigma \\ \log \mathcal{N}(x | \mu, \sigma^2) &= \log C_{\mathcal{N}} - \frac{1}{2 \sigma^2} (x - \mu)^2 \end{aligned} $$

を計算します。計算結果の指数をとると確率が得られます。

$$ \mathcal{N}(x | \mu, \sigma^2) = \exp \Bigr( \log \mathcal{N}(x | \mu, \sigma^2) \Bigr) $$

 指数と対数の性質より$\exp(\log x) = x$です。

 次は、関数を使って確率を計算します。
 1次元ガウス分布の確率密度関数dnorm()を使って計算します。

# ガウス分布の関数により確率密度を計算
dens <- dnorm(x = x, mean = mu, sd = sigma)
dens
## [1] 0.1595769

 変数の引数xx、平均の引数meanmu、標準偏差の引数sdsigmaを指定します。

 log = TRUEを指定すると対数をとった確率密度を返します。

# ガウス分布の対数をとった関数により確率密度を計算
log_dens <- dnorm(x = x, mean = mu, sd = sigma, log = TRUE)
dens <- exp(log_dens)
dens; log_dens
## [1] 0.1595769
## [1] -1.835229

 計算結果の指数をとると確率密度が得られます。

統計量の計算

 1次元ガウス分布の平均と分散を計算します。

 平均を計算します。

# 平均を計算
E_x <- mu
E_x
## [1] 1

 1次元ガウス分布の平均は、定義より$\mu$です。

$$ \mathbb{E}[x] = \mu $$

 分散を計算します。

# 分散を計算
V_x <- sigma^2
V_x
## [1] 6.25

 1次元ガウス分布の分散は、定義より$\sigma$の2乗です。

$$ \mathbb{V}[x] = \sigma^2 $$


分布の可視化

 ggplot2パッケージを利用して1次元ガウス分布のグラフを作成します。

 1次元ガウス分布の確率変数がとり得る値$x$ごとの確率を計算します。

# 平均を指定
mu <- 0

# 標準偏差を指定
sigma <- 1

# 作図用のxの点を作成
x_vals <- seq(from = mu - sigma*4, to = mu + sigma*4, length.out = 250)

# 1次元ガウス分布を計算
dens_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  density = dnorm(x = x_vals, mean = mu, sd = sigma) # 確率密度
)
head(dens_df)
## # A tibble: 6 x 2
##       x  density
##   <dbl>    <dbl>
## 1 -4    0.000134
## 2 -3.97 0.000152
## 3 -3.94 0.000173
## 4 -3.90 0.000196
## 5 -3.87 0.000222
## 6 -3.84 0.000251

 $x$がとり得る値を作成してx_valsとします。この例では、平均muを中心に標準偏差sigma4倍を範囲とします。
 x_valsの各要素に対応する確率密度を求めてデータフレームに格納します。

  1次元ガウス分布のグラフを作成します。

# 1次元ガウス分布を作図
ggplot(data = dens_df, mapping = aes(x = x, y = density)) + # データ
  geom_line(color = "#00A968") + # 折れ線グラフ
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("mu=", mu, ", sigma=", sigma)) # ラベル

f:id:anemptyarchive:20220130175143p:plain
1次元ガウス分布のグラフ


 この分布に平均と最頻値、標準偏差の情報を重ねて表示します。

# 統計量を重ねた1次元ガウス分布を作図
ggplot(data = dens_df, mapping = aes(x = x, y = density)) + # データ
  geom_line(color = "#00A968") + # 分布
  geom_vline(xintercept = mu, color = "orange", size = 1, linetype = "dashed") + # 平均
  geom_vline(xintercept = mu - sigma, color = "orange", size = 1, linetype = "dotted") + # 平均 - 標準偏差
  geom_vline(xintercept = mu + sigma, color = "orange", size = 1, linetype = "dotted") + # 平均 + 標準偏差
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("mu=", mu, ", sigma=", sigma)) # ラベル

f:id:anemptyarchive:20220130175156p:plain
1次元ガウス分布のグラフ

 1次元ガウス分布のグラフを描画できました。

パラメータと分布の形状の関係

 パラメータ$\mu, \sigma$が及ぼす分布への影響をアニメーション(gif画像)で可視化します。

平均の影響

 平均パラメータ$\mu$の値を少しずつ変更して、分布の変化をアニメーションで確認します。

・作図コード(クリックで展開)

# muとして利用する値を指定
mu_vals <- seq(from = -5, to = 5, by = 0.1)
#length(mu_vals) # フレーム数

# 標準偏差を指定
sigma <- 1

# 作図用のxの点を作成
x_vals <- seq(from = median(mu_vals) - sigma*4, to = median(mu_vals) + sigma*4, length.out = 250)

# muの値ごとに分布を計算
anime_dens_df <- tidyr::tibble()
for(mu in mu_vals) {
  # 1次元ガウス分布を計算
  tmp_dens_df <- tidyr::tibble(
    x = x_vals, # 確率変数
    density = dnorm(x = x_vals, mean = mu, sd = sigma), # 確率密度
    parameter = paste0("mu=", round(mu, 1), ", sigma=", sigma) %>% 
      as.factor() # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_dens_df <- rbind(anime_dens_df, tmp_dens_df)
}
head(anime_dens_df)
## # A tibble: 6 x 3
##       x density parameter     
##   <dbl>   <dbl> <fct>         
## 1 -4      0.242 mu=-5, sigma=1
## 2 -3.97   0.234 mu=-5, sigma=1
## 3 -3.94   0.226 mu=-5, sigma=1
## 4 -3.90   0.219 mu=-5, sigma=1
## 5 -3.87   0.211 mu=-5, sigma=1
## 6 -3.84   0.203 mu=-5, sigma=1

 $\mu$がとり得る値を作成してmu_valsとします。
 for()ループを使ってmu_valsの値ごとに確率密度を計算して、anime_dens_dfに追加していきます。

 gif画像を作成します。

# アニメーション用の1次元ガウス分布を作図
anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = x, y = density)) + # データ
  geom_line(color = "#00A968") + # 折れ線グラフ
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gaussian Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_dens_graph, nframes = length(mu_vals), fps = 100)


f:id:anemptyarchive:20220130175216g:plain
平均と1次元ガウス分布の形状の関係

 $\mu$が大きくなるに従って、$x$が大きいほど確率密度が大きくなる(山が右に移動する)のを確認できます。

標準偏差の影響

 標準偏差パラメータ$\sigma$の値を少しずつ変更して、分布の変化をアニメーションで確認します。

・作図コード(クリックで展開)

# sigmaとして利用する値を指定
sigma_vals <- seq(from = 1, to = 10, by = 0.1)
#length(sigma_vals) # フレーム数

# 平均を指定
mu <- 0

# 作図用のxの点を作成
x_vals <- seq(from = mu - max(sigma_vals)*2, to = mu + max(sigma_vals)*2, length.out = 250)

# sigmaの値ごとに分布を計算
anime_dens_df <- tidyr::tibble()
for(sigma in sigma_vals) {
  # 1次元ガウス分布を計算
  tmp_dens_df <- tidyr::tibble(
    x = x_vals, # 確率変数
    density = dnorm(x = x_vals, mean = mu, sd = sigma), # 確率密度
    parameter = paste0("mu=", mu, ", sigma=", round(sigma, 1)) %>% 
      as.factor() # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_dens_df <- rbind(anime_dens_df, tmp_dens_df)
}
head(anime_dens_df)
## # A tibble: 6 x 3
##       x  density parameter    
##   <dbl>    <dbl> <fct>        
## 1 -20   5.52e-88 mu=0, sigma=1
## 2 -19.8 1.35e-86 mu=0, sigma=1
## 3 -19.7 3.24e-85 mu=0, sigma=1
## 4 -19.5 7.54e-84 mu=0, sigma=1
## 5 -19.4 1.71e-82 mu=0, sigma=1
## 6 -19.2 3.79e-81 mu=0, sigma=1

 $\sigma$がとり得る値を作成してsigma_valsとします。
 for()ループを使ってsigma_valsの値ごとに確率密度を計算して、anime_dens_dfに追加していきます。

 gif画像を作成します。

# アニメーション用の1次元ガウス分布を作図
anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = x, y = density)) + # データ
  geom_line(color = "#00A968") + # 折れ線グラフ
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gaussian Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_dens_graph, nframes = length(sigma_vals), fps = 100)


f:id:anemptyarchive:20220130175254g:plain
標準偏差と1次元ガウス分布の形状の関係

 $\sigma$が大きくなるに従って、裾が広く(山が低く)なるのを確認できます。

乱数の生成

 1次元ガウス分布の乱数を生成してヒストグラムを確認します。

 パラメータを指定して、1次元ガウス分布に従う乱数を生成します。

# 平均を指定
mu <- 1

# 標準偏差を指定
sigma <- 2.5

# データ数を指定
N <- 1000

# 1次元ガウス分布に従う乱数を生成
x_n <- rnorm(n = N, mean = mu, sd = sigma)
x_n[1:5]
## [1]  3.80204201  0.03885278  4.61619643 -1.69259087 -0.01516701

 1次元ガウス分布の乱数生成関数rnorm()のデータ数(サンプルサイズ)の引数nN、平均の引数meanmu、標準偏差の引数sdsigmaを指定します。

 作図に利用するため、サンプルをデータフレームに格納し、また$x$の値と分布を作成しておきます。

# サンプルを格納
data_df <- tidyr::tibble(x = x_n)

# 作図用のxの点を作成
x_vals <- seq(from = mu - sigma*4, to = mu + sigma*4, length.out = 250)

# 1次元ガウス分布を計算
dens_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  density = dnorm(x = x_vals, mean = mu, sd = sigma) # 確率密度
)


 ヒストグラムを作成します。

# サンプルの頻度を作図
ggplot(data = data_df, mapping = aes(x = x)) + # 
  geom_histogram(breaks = seq(from = min(x_vals), to = max(x_vals), length.out = 50), fill = "#00A968") + # 
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("mu=", mu, ", sigma=", sigma, ", N=", N)) # ラベル

f:id:anemptyarchive:20220130175339p:plain
1次元ガウス分布の乱数のヒストグラム:頻度

 geom_histogram()でヒストグラムを作成します。breaks引数に区切り位置を指定します。この例では、seq()を使ってx_valsの最小値から最大値の範囲をlength.out引数に指定した数に区切ります。

 サンプルの密度を分布と重ねて描画します。

# サンプルの密度を格納
ggplot(data = data_df, mapping = aes(x = x, y = ..density..)) + # 
  geom_histogram(breaks = seq(from = min(x_vals), to = max(x_vals), length.out = 50), 
                 fill = "#00A968") + # 
  geom_line(data = dens_df, mapping = aes(x = x, y = density), 
            color = "darkgreen", linetype = "dashed") + # 
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("mu=", mu, ", sigma=", sigma, ", N=", N)) # ラベル

f:id:anemptyarchive:20220130175359p:plain
1次元ガウス分布の乱数のヒストグラム:密度

 geom_histogram()y = ..density..を指定すると、頻度を密度に変換して描画します。

 データ数が十分に増えると元の分布(破線のグラフ)に形が近付きます。

 サンプルサイズとヒストグラムの変化をアニメーションで確認します。

・作図コード(クリックで展開)

# フレーム数を指定
N_frame <- 300

# 乱数を1つずつ生成
lambda_n <- rep(NA, times = N_frame)
anime_freq_df <- tidyr::tibble()
anime_data_df <- tidyr::tibble()
anime_dens_df <- tidyr::tibble()
for(n in 1:N_frame) {
  # 1次元ガウス分布に従う乱数を生成
  x_n[n] <- rnorm(n = 1, mean = mu, sd = sigma)
  
  # ラベル用のテキストを作成
  label_text <- paste0("mu=", mu, ", sigma=", sigma, ", N=", n)
  
  # n個の乱数を格納
  tmp_freq_df <- tidyr::tibble(
    x = x_n[1:n], # サンプル
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # n番目の乱数を格納
  tmp_data_df <- tidyr::tibble(
    x = x_n[n], # サンプル
    parameter = as.factor(label_text) # フレーム切替用のラベル
  )
  
  # n回目のラベルを付与
  tmp_dens_df <- dens_df %>% 
    dplyr::mutate(parameter = as.factor(label_text))
  
  # 結果を結合
  anime_freq_df <- rbind(anime_freq_df, tmp_freq_df)
  anime_data_df <- rbind(anime_data_df, tmp_data_df)
  anime_dens_df <- rbind(anime_dens_df, tmp_dens_df)
}
head(anime_freq_df)
## # A tibble: 6 x 2
##        x parameter           
##    <dbl> <fct>               
## 1 -0.295 mu=1, sigma=2.5, N=1
## 2 -0.295 mu=1, sigma=2.5, N=2
## 3  1.56  mu=1, sigma=2.5, N=2
## 4 -0.295 mu=1, sigma=2.5, N=3
## 5  1.56  mu=1, sigma=2.5, N=3
## 6  0.771 mu=1, sigma=2.5, N=3

 乱数を1つずつ生成して、結果をanime_***_dfに追加していきます。

 ヒストグラムのアニメーションを作成します。

# アニメーション用のサンプルの頻度を作図
anime_freq_graph <- ggplot() + # 
  geom_histogram(data = anime_freq_df, mapping = aes(x = x), 
                 breaks = seq(from = min(x_vals), to = max(x_vals), length.out = 50), 
                 fill = "#00A968") + # 頻度
  geom_point(data = anime_data_df, mapping = aes(x = x, y = 0), 
             color = "orange", size = 5) + # サンプル
  gganimate::transition_manual(parameter) + # フレーム
  labs(title = "Gaussian Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_freq_graph, nframes = N_frame, fps = 100)


 密度のアニメーションを作成します。

# アニメーション用のサンプルの密度を作図
anime_prop_graph <- ggplot() + # 
  geom_histogram(data = anime_freq_df, mapping = aes(x = x, y = ..density..), 
                 breaks = seq(from = min(x_vals), to = max(x_vals), length.out = 50), 
                 fill = "#00A968") + # 密度
  geom_line(data = anime_dens_df, mapping = aes(x = x, y = density), 
            color = "darkgreen", linetype = "dashed") + # 元の分布
  geom_point(data = anime_data_df, mapping = aes(x = x, y = 0), 
             color = "orange", size = 5) + # サンプル
  gganimate::transition_manual(parameter) + # フレーム
  ylim(c(0, max(dens_df[["density"]]) + 0.3)) + # y軸の表示範囲
  labs(title = "Gaussian Distribution", 
       subtitle = "{current_frame}") # ラベル

# gif画像を作成
gganimate::animate(anime_prop_graph, nframes = N_frame, fps = 100)


f:id:anemptyarchive:20220130182240g:plainf:id:anemptyarchive:20220130182027g:plain
1次元ガウス分布の乱数の推移

 サンプルが増えるに従って、元の分布に近付くのを確認できます。

分布の生成

 1次元ガウス分布が共役事前分布となる1次元ガウス分布の平均パラメータを生成して分布を作図します。

 平均$\mu_{\mathrm{prior}}$・標準偏差$\sigma_{\mathrm{prior}}$・分散$\sigma_{\mathrm{prior}}^2$の1次元ガウス分布

$$ \mathcal{N}(\mu | \mu_{\mathrm{prior}}, \sigma_{\mathrm{prior}}^2) = \frac{1}{\sqrt{2 \pi \sigma_{\mathrm{prior}}^2}} \exp \left( - \frac{(x - \mu_{\mathrm{prior}})^2}{2 \sigma_{\mathrm{prior}}^2} \right) $$

を事前分布として、$\mu$を生成します。
 生成した$\mu$を平均として用いて、標準偏差$\sigma$・分散$\sigma^2$の1次元ガウス分布

$$ \mathcal{N}(x | \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right) $$

を描画します。

 平均パラメータ$\mu$を生成します。

# 超パラメータを指定
mu_prior <- 1
sigma_prior <- 2.5

# サンプルサイズを指定
N <- 10

# 1次元ガウス分布の平均パラメータを生成
mu_n <- rnorm(n = N, mean = mu_prior, sd = sigma_prior)
mu_n[1:5]
## [1]  2.0383114  4.3423351  0.4718459 -0.8101986  0.5675290

 事前ガウス分布のパラメータ(超パラメータ)$\mu_{\mathrm{prior}}, \sigma_{\mathrm{prior}}$を指定して、1次元ガウス分布に従う乱数を生成し、平均パラメータ$\mu$として利用します。

 まずは、目安となるように$\mu$の期待値$\mathbb{E}[\mu]$による分布を求めます。

# 平均パラメータの期待値を計算
E_mu <- mu_prior

# 標準偏差パラメータを指定
sigma <- 1

# 作図用のxの点を作成
x_vals <- seq(from = E_mu - sigma*5, to = E_mu + sigma*5, length.out = 250)

# 平均パラメータの期待値による1次元ガウス分布を計算
res_dens_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  density = dnorm(x = x_vals, mean = E_mu, sd = sigma), # 確率密度
  parameter = paste0("E[mu]=", E_mu) %>% 
    as.factor() # 色分け用のラベル
)
head(res_dens_df)
## # A tibble: 6 x 3
##       x    density parameter
##   <dbl>      <dbl> <fct>    
## 1 -4    0.00000149 E[mu]=1  
## 2 -3.96 0.00000182 E[mu]=1  
## 3 -3.92 0.00000221 E[mu]=1  
## 4 -3.88 0.00000270 E[mu]=1  
## 5 -3.84 0.00000328 E[mu]=1  
## 6 -3.80 0.00000398 E[mu]=1

 平均パラメータの期待値$\mathbb{E}[\mu] = \mu_{\mathrm{prior}}$をE_lambdaとします。
 生成された1次元ガウス分布の標準偏差パラメータ$\sigma$を(既知の値として)指定します。

 $x$として利用する値を作成してx_valsとします。この例では、平均E_muを中心に標準偏差のsigma5倍を範囲とします。
 dnorm()で1次元ガウス分布の確率密度を計算します。平均の引数meanE_mu、標準偏差の引数sdsigmaを指定します。

 次に、パラメータのサンプルmu_nの値ごとに分布を計算します。

# サンプルごとに分布を計算
for(mu in mu_n) {
  # 1次元ガウス分布を計算
  tmp_dens_df <- tidyr::tibble(
    x = x_vals, # 確率変数
    density = dnorm(x = x_vals, mean = mu, sd = sigma), # 確率密度
    parameter = paste0("mu=", round(mu, 2)) %>% 
      as.factor() # 色分け用のラベル
  )
  
  # 結果を結合
  res_dens_df <- rbind(res_dens_df, tmp_dens_df)
}

 各分布の情報をres_dens_dfに追加します。

 $N + 1$個の1次元ガウス分布を作図します。

# サンプルによる1次元ガウス分布を作図
ggplot() + 
  geom_line(data = res_dens_df, 
            mapping = aes(x = x, y = density, color = parameter, 
                          alpha = parameter, linetype = parameter), size = 0.8) + # 分布
  scale_linetype_manual(values = c("dashed", rep("solid", times = N))) + # 線の種類
  scale_alpha_manual(values = c(1, rep(0.5, times = N))) + # 透過度
  labs(title = "Gaussian Distribution", 
       subtitle = paste0("mu_pri=", mu_prior, ", sigma_pri=", sigma_prior, ", sigma=", sigma), 
       x = expression(lambda)) # ラベル

f:id:anemptyarchive:20220130175556p:plain
サンプルしたパラメータによる1次元ガウス分布のグラフ

 期待値による分布(破線)を中心に分布しています。

参考文献

  • 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.

おわりに

 ヒストグラムを描くときの頻度から密度に変換する計算をよく分かってない。

 2022年1月30日は、Juice=Juiceの江端妃咲さんの15歳のお誕生日です!

 若い、眩しい、キラキラ粒子が舞ってるのを感じる。

【次の内容】

www.anarchive-beta.com