見出し画像

「ベイズ統計モデリングによるデータ分析入門」をPythonとPyMCで写経 ～ Vol.10 デザイン行列を用いた単回帰モデルの推定

ネイピア DS

書籍の著者　馬場真哉 先生

この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第３部第４章「デザイン行列を用いた一般化線形モデルの推定」の Python 写経活動記録です。

書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
今回も引き続き「単回帰モデル」を堪能します。
書籍によると「formula構文」と「デザイン行列による表現」に慣れる章です。

では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀


はじめに

このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

■ 記事の範囲
この記事はテキスト第３部第４章の以下の節を取り扱います。

4.2 分析の準備
4.4 formula構文を用いたデザイン行列の作成
4.5 デザイン行列を使うためのStanファイルの修正
4.6 MCMCの実行

ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部、第３部第１章をご覧いただき、本記事との繋がりをご確認下さいませ。

■ コード記述法
Jupyter Notebook 形式でコードを記述します。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np
import pandas as pd

# ベイズ統計モデリング
import pymc as pm                      # pymc
import arviz as az                     # 分析・可視化

# デザイン行列
from patsy import dmatrices

# 統計処理
import scipy.stats as stats            # 正規分布乱数の生成

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()                        # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'

第４章 デザイン行列を用いた一般化線形モデルの推定

データの読み込み

テキスト 4.2 節に相当します。
テキストの仮想のビールの売り上げデータを引用いたします。
前々回記事から利用しているデータです。
csv ファイルを pandas のデータフレーム形式で変数 file_beer_sales_2 に読み込みます。

# p.180 分析対象データの読み込み

# ファイルの読み込み
file_beer_sales_2 = pd.read_csv('./data/3-2-1-beer-sales-2.csv')

# 結果の表示
print('file_beer_sales_2.shape: ', file_beer_sales_2.shape)
file_beer_sales_2.head(3)

【実行結果】
標本サイズ 100、変数 sales は 売り上げ（単位：万円）、temperature は気温（おそらく摂氏℃）です。
気温とビール売り上げの関係を単回帰モデルで分析します。

🚀🚀🚀

デザイン行列と係数ベクトル

前回記事までの単回帰モデルでは、線形予測子 $${\mu}$$ の切片と傾きに別々の係数を当てて表現していました。

$$
\mu_i = \underbrace{\beta_0}_{切片} + \underbrace{\beta_1}_{傾き} \underbrace{x_i}_{説明変数}
$$

テキスト p.181 式(3.41) の数式を一部改変して引用

上式は識別子 $${i}$$ 単位であり、標本サイズ $${n}$$ 個の数式が並びます。
これをベクトル・行列でひとまとめにして表現してみます。

$$
\underbrace{\begin{bmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_n \end{bmatrix}}_{\bm \mu}
= \begin{bmatrix} \beta_0 + \beta_1 x_1 \\ \beta_0 + \beta_1 x_2 \\ \vdots \\ \beta_0 + \beta_1 x_n\end{bmatrix}
= \underbrace{\begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots \\ 1 & x_n \end{bmatrix}}_{\bm X} 
\underbrace{\begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix} }_{\bm \beta}
= \bm X \bm \beta = \bm \mu \\
$$

テキスト p.182 式(3.46) を参考にして構成

$${\bm X}$$ は説明変数の行列であり、書籍は「デザイン行列」と呼んでいます。

  • １列目の $${1}$$ のみの列ベクトルは定数項です。
    切片 $${\beta_0}$$ と掛け算する対象です。

  • ２列目は 説明変数ベクトル $${\bm x}$$ です。
    単回帰モデルでは説明変数が１個なので１列＝ベクトルですが、重回帰モデルなどの説明変数が複数存在するモデルでは説明変数が複数列＝行列になって $${\bm X}$$ に合流します。

$${\bm \beta}$$ は係数ベクトルです。
単回帰モデルでは切片 $${\beta_0}$$ と傾き $${\beta_1}$$ の２要素で構成されますが、重回帰モデルなどの説明変数が複数存在するモデルでは説明変数の数 $${k}$$ 個の係数要素が合流します。

🚀🚀🚀

デザイン行列とformula構文

テキスト 4.4 節に相当します。
テキストにならい、formula 構文を活用してデザイン行列を作成してみます。

formula 構文のオリジナルは R 言語（だと思うの）です。
Python には R 言語ライクに formula 構文を使ってデザイン行列を作成するライブラリがいくつかあります。
ここでは patsy ライブラリを利用します。
pandas データフレームの列名（＝変数名）を用いて、目的変数と説明変数等の関係を formula 構文 で表現し、デザイン行列等を作成できます。

patsy 公式サイトの Quick Startページのリンクはこちら

今回の例題データは 目的変数：sales、説明変数：temperature であり、formula 構文で単回帰モデルを記述すると：

$$
\mathtt{sales} \sim \mathtt{temperature}
$$

です。
チルダ $${\sim}$$ の左辺に目的変数、右辺に説明変数を書きます。
上式の書き方の場合、定数項を含む formula と同じ意味合いになります。

$$
\underbrace{\mathtt{sales}}_{目的変数} \sim \underbrace{\mathtt{1}}_{定数項} + \underbrace{\mathtt{temperature}}_{説明変数}
$$

実はこの formula 構文は、前々回・前回記事の単回帰分析の際にこっそり使っています。
statsmodels は裏側で patsy ライブラリを用いて formula 構文 ⇒ デザイン行列作成の機能を実現しています。

formula を定義して、デザイン行列を作成しましょう。
patsy の dmatrices 関数を利用します。
目的変数とデザイン行列を生成して、先頭５行を表示します。

# p.182 デザイン行列の作成

# formula構文の設定
formula_lm = 'sales ~ temperature'

# 目的変数Y, デザイン行列（説明変数）Xの作成
Y_dm, X_dm = dmatrices(formula_lm, file_beer_sales_2, return_type='dataframe')

# デザイン行列の先頭5行の表示
X_dm.head()

【実行結果】
こちらはデザイン行列 $${\bm X}$$ を体現したデータです。
定数項の変数名が「Intercept」（切片）になっています。

# 目的変数の先頭5行の表示
Y_dm.head()

【実行結果】
こちらは目的変数のデータです。

【コード補足】
こちらは formula の定義です。変数 formula_lm に設定しています。

# formula構文の設定
formula_lm = 'sales ~ temperature'

こちらは dmatrices 関数に formula を与えて、目的変数 Y_dm、デザイン行列 X_dm を作成しています。
３つの引数は、①formula、②データフレーム名、③戻り値をデータフレームで返す設定、です。

# 目的変数Y, デザイン行列（説明変数）Xの作成
Y_dm, X_dm = dmatrices(formula_lm, file_beer_sales_2, return_type='dataframe')

🚀🚀🚀

ベイズモデリング by PyMC

テキスト 4.5 節に相当します。

① モデルの概要
デザイン行列と係数ベクトルを用いて、次のモデルを実装します。

$$
\begin{align*}
Y &\sim \text{Normal}\ (\bm \mu,\ \sigma^2) \\
\bm \mu &= \bm X \bm \beta \\
\bm \beta &\sim \text{Normal}\ (0,\ (1e5)^2,\ \text{dims=2}) \\
\sigma &\sim \text{HalfNormal}\ ((1e5)^2) \\
\end{align*}
$$

目的変数 $${Y}$$（sales）は正規分布に従うと仮定しています。
正規分布の平均パラメータ $${\bm \mu}$$ は、デザイン行列 $${\bm X}$$ と係数ベクトル $${\bm \beta}$$ の積です。
係数ベクトル $${\bm \beta}$$ は平均０、標準偏差 100000 の正規分布（2次元）に従うとします。
標準偏差パラメータ $${\sigma}$$ は標準偏差 100000 の半正規分布に従うとします。

② PyMC のモデル設定
上式のモデルを PyMC で記述します。

# モデリング

# coordsの設定
coords = {'id': file_beer_sales_2.index.values,   # 観測データの識別子
          'coefs': ['intercept', 'temperature']}  # 係数ベクトルbの識別子

# モデルの定義
with pm.Model(coords=coords) as model:
    
    ## dataの設定
    # 目的変数: 売上データ
    Y = pm.Data('Y', value=Y_dm.values.flatten(), dims='id')
    # 説明変数: デザイン行列
    X = pm.Data('X', value=X_dm.values, dims=('id', 'coefs'))

    ## 事前分布: 無情報事前分布的な分布
    # 係数ベクトル
    b = pm.Normal('b', mu=0, sigma=1e5, dims='coefs')
    # 標準偏差
    sigma = pm.HalfNormal('sigma', sigma=1e5)

    ## 線形予測子：デザイン行列と係数ベクトルの積
    mu = pm.Deterministic('mu', X @ b, dims='id')

    ## 尤度関数: 正規分布を仮定
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='id')

【実行結果】なし

【コードの補足】
説明変数（デザイン行列）は２次元 $${\text{dims}=(\mathtt{id},\ \mathtt{coefs})}$$ です。
coords「$${\mathtt{coefs}}$$」は係数ベクトルの２つの要素の名前 $${\mathtt{['intercept', 'temperature']}}$$ も定義しています。

# 説明変数: デザイン行列
X = pm.Data('X', value=X_dm.values, dims=('id', 'coefs'))

線形予測子 $${\mathtt{mu}}$$ は デザイン行列と係数ベクトルの積です。

# 線形予測子：デザイン行列と係数ベクトルの積
mu = pm.Deterministic('mu', X @ b, dims='id')

🚀

モデルの数式ライク表示とグラフィカルモデル描画をします。

# モデルの表示
model

【実行結果】
最下行の obs が尤度、上の２行が事前分布です。
数式表現では次元や線形予測子の中身が分からないですね（汗）

# モデルの可視化
pm.model_to_graphviz(model)

【実行結果】
説明変数（デザイン行列）の X の次元は $${100 \times 2}$$、係数ベクトル b は $${2}$$、線形予測子 mu・尤度 obs・目的変数（Data）の Y は、軸 id を設定した 100 個のデータを持っています。
説明変数 X と係数ベクトル b は線形予測子の形式で mu のパラメータになっています。
mu と sigma は尤度 obs のパラメータになっています。

🚀🚀🚀

MCMC の実行

テキスト 4.6 節に相当します。
MCMCを実行しましょう。
NUTS サンプラーに nutpie を利用します。

%%time
# p.183 MCMCの実行 ※高速化したいときは nuts_sampler='nutpie'
with model:
    idata = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1,
                      nuts_sampler='nutpie')

【実行結果】
Divergences（ダイバージェンス）は０件です。

収束の確認をします。
MCMC サンプルの要約表を表示します。

# p.184 要約統計量の表示
var_names = ['b', 'sigma']
az.summary(idata, var_names=var_names, hdi_prob=0.95)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。

係数ベクトル b には係数の名前 intercept、temperature が付記されていて分かりやすいです！
これが、coords と dims の設定の力なのです💪💪💪

前回記事の単回帰モデルの要約表と比べてみます。
（参考：前回記事の単回帰モデルの要約表）

事後平均（mean）、標準誤差（sd）、95% HDI（hdi_2.5%、hdi_97.5%）の値は、若干異なる結果になっています。

結果が変わる原因を Gemini に訊いてみました。
長文ですがぜひお読み下さい！

PyMCにおける単回帰モデルで、係数（$${b_{0},b_{1}}$$）を個別に定義する場合とデザイン行列 $${\bm X}$$ を用いて一括定義する場合で要約統計量に差が生じる主な理由は、モデルの数学的な「構造」自体は同じでも、MCMCサンプリングにおける挙動や数値計算上のプロセスが異なるためです。
主な理由は以下の通りです。

  1. サンプリングのランダム性と収束性
    MCMC（マルコフ連鎖モンテカルロ法）は乱数を用いたシミュレーションであるため、計算のたびに結果はわずかに変動します。

    • 初期値の影響：定義方法が異なると、サンプラーが探索を開始する初期値や、パラメータ空間を探索する「歩幅」が微妙に変化します。

    • 連鎖（Chain）の独立性：PyMCがデフォルトで使用するNUTS（No-U-Turn Sampler）などのアルゴリズムでは、パラメータがベクトルとして一括定義されているか個別かによって、勾配計算の順序や効率が変わることがあり、それが最終的なサンプルの分布（事後分布）に微細な差異をもたらします。

  2. 事前分布の暗黙的な相関

    • 独立性の仮定：$${b_{0}}$$ と $${b_{1}}$$ を個別に定義した場合、サンプラーはこれらを完全に独立した変数として扱います。

    • ベクトル定義：デザイン行列を用いる形式（例：pm.Normal('beta', ..., shape=2)）で一括定義すると、数学的には独立であっても、計算上は一つのベクトルオブジェクトとして処理されます。多変量正規分布などを事前分布に用いる場合、パラメータ間の共分散を考慮しやすくなり、探索の効率（ハミルトニアンの軌道など）が向上し、結果として推定値の精度（標準誤差など）に影響を与えることがあります。

  3. 多重共線性への耐性と幾何学的構造
    回帰モデルにおいて、切片 $${b_{0}}$$ と傾き $${b_{1}}$$ はしばしば強い負の相関を持ちます（一方が増えれば他方が減ることでデータに適合しようとするため）。 

    • 幾何学的探索：デザイン行列を用いる形式は、高次元のパラメータ空間をベクトルとして効率よく探索するのに適しています。個別に定義すると、サンプラーがこれらの相関関係を解明するのに時間がかかり、収束が不十分な場合に要約統計量（特に標準誤差や有効サンプルサイズ ESS）に差が出やすくなります。

  4. 数値計算上の誤差
    デザイン行列を用いた行列演算（例：pm.math.dot(X, beta)）と、個別のスカラー演算（例：b0 + b1 * x）では、浮動小数点演算の順序や内部的な最適化が異なります。これにより、極めて微小ではありますが、対数尤度の計算結果に差が生じ、それがサンプリングの採択判定に影響を及ぼす可能性があります。

結論として：十分なサンプル数（例：数千〜数万）を確保し、モデルが完全に収束（R-hatが1.01未満など）していれば、これらの差は実用上無視できるほど小さくなります。もし大きな差がある場合は、サンプリングの収束不足や事前分布の設定ミス（片方にだけ異なるスケールを適用している等）を確認することをお勧めします。 

差が生じるのは仕方がないことのようです。
大きな差が無いことを確認したので、次に進みましょう🚀

トレースプロットを描画します。

# トレースプロットの描画
az.plot_trace(idata, compact=False, var_names=var_names,
              backend_kwargs={'tight_layout': True});

【実行結果】
右側のチャートがゲジゲジしています。

以上のチェックに基づいて、収束していると考えましょう。
テキスト 第４章 はここまでとなります。

🚀🚀🚀

回帰直線の可視化

せっかく単回帰モデルを構築したので…
第５章を先取りして、回帰直線＋95%信用区間、回帰直線＋95%予測区間のチャートを描画しましょう！
テキスト p.199 ～ p.200 に相当します。

描画に利用するデータは MCMC サンプルを PyMC の外で加工して作成します。

① 回帰直線＋95%信用区間
単回帰モデルの予測平均（＝線形予測子）$${\mu}$$ をMCMCサンプルから算出して、回帰直線と95%信用区間を描画します。

# p.199 図3.5.3 回帰直線95%信用区間付き

## データの作成
# MCMCサンプルからβ0, β1, σを取り出し
beta0s, beta1s = az.extract(idata.posterior).b.data
sigmas = az.extract(idata.posterior).sigma.data
# x軸の値の設定
x_val = np.linspace(X_dm.iloc[:, 1].min(), X_dm.iloc[:, 1].max(), 1001)
# μの算出
mu_preds = (beta0s + np.outer(x_val, beta1s)).T
# μの95%信用区間の算出
mu_95ci = np.quantile(mu_preds, q=[0.025, 0.975], axis=0).T

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
plt.scatter(X_dm['temperature'], Y_dm['sales'], color='navy', s=10,
            label='観測値')
# サンプルの平均値の描画
plt.plot(x_val, mu_preds.mean(axis=0), label='回帰直線')
# サンプルの95%信用区間の塗りつぶし
plt.fill_between(x_val, mu_95ci[:, 0], mu_95ci[:, 1], alpha=0.2,
                 label='95%信用区間')
# 修飾
plt.xlabel('温度 [℃]', fontsize=12)
plt.ylabel('売上', fontsize=12)
plt.xticks(range(10, 31, 2))
plt.legend();

【実行結果】
予測平均 $${\mu}$$ の 95%信用区間は、予測値の平均（期待値）が 95% の確率で収まる区間です。
かなり狭い範囲なんだと感じます。

【コードの補足】
◆ コード１：パラメータの事後分布からのサンプルを取得

# MCMCサンプルからβ0, β1, σを取り出し
beta0s, beta1s = az.extract(idata.posterior).b.data
sigmas = az.extract(idata.posterior).sigma.data

こちらのコードは MCMC サンプルを格納した idata に対して arviz の extract　関数を適用して、平坦化（１次元化）した上で、data 属性で各パラメータの MCMC サンプルを numpy 配列化しています。

◆ コード２：予測平均 $${\mu}$$（＝線形予測子）の算出

# x軸の値の設定
x_val = np.linspace(X_dm.iloc[:, 1].min(), X_dm.iloc[:, 1].max(), 1001)
# μの算出
mu_preds = (beta0s + np.outer(x_val, beta1s)).T

こちらのコードは $${\mu = \beta_0 + \beta_1 x}$$ を計算しています。

◆ コード３：$${\mu}$$ の 95% 信用区間の算出

# μの95%信用区間の算出
mu_95ci = np.quantile(mu_preds, q=[0.025, 0.975], axis=0).T

コード２で計算した $${\mu}$$（shape=(4000, 1001)）のサンプルから 2.5% 点と 95.5% 点を算出します。

🚀

② 回帰直線＋95%予測区間
$${\mu}$$ の（MCMC的）サンプルに観測値の「ばらつき」（正規分布乱数）を加味して「予測分布からのサンプルに相当するデータ」を作成し、回帰直線と95%予測区間を描画します。

# p.200 図3.5.4 回帰直線95%予測区間付き

## MCMCサンプルデータからyの予測値を算出
# yの予測値の算出 ※scipy.statsで正規分布乱数を付加
y_preds = (mu_preds.T + stats.norm.rvs(loc=0, scale=sigmas, random_state=123)).T
# yの95%予測区間の算出
y_95ci = np.quantile(y_preds, q=[0.025, 0.975], axis=0).T

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
plt.scatter(X_dm['temperature'], Y_dm['sales'], color='navy', s=10,
            label='観測値')
# yの予測値の平均値の描画
plt.plot(x_val, y_preds.mean(axis=0), color='tab:red', label='回帰直線')
# yの予測値の95%予測区間の塗りつぶし
plt.fill_between(x_val, y_95ci[:, 0], y_95ci[:, 1], color='lightpink', alpha=0.3,
                 label='95%予測区間')
# 修飾
plt.xlabel('温度 [℃]', fontsize=12)
plt.ylabel('売上', fontsize=12)
plt.xticks(range(10, 31, 2))
plt.legend();

【実行結果】
ばらつきが加わるので 95% 予測区間の幅はかなり広くなります。

前回記事の事後予測 95% HDI とよく似ています（当たり前）。
今回記事の方が縁が滑らかな理由は、単に x 軸の値の粒度が細かいからです。

（参考：前回記事の事後予測 95% HDI）

【コードの補足】
◆ コード：y の予測値（観測値ベース）の算出

# yの予測値の算出 ※scipy.statsで正規分布乱数を付加
y_preds = (mu_preds.T + stats.norm.rvs(loc=0, scale=sigmas, random_state=123)).T

$${\mathtt{scipy.stats}}$$ の $${\mathtt{norm.rvs()}}$$ で正規分布乱数（平均：０、標準偏差：$${\sigma}$$ のMCMCサンプル）を予測平均 $${\mu}$$ に付加して、観測値ベース（ばらつき込み）の予測値を算出します。

以上で デザイン行列を用いた単回帰モデルの学びを終わりにします。
面白かったですね。

今回の記事は以上です。
楽しかったですね！

シリーズの記事

次の記事（Bambi版）

前の記事

Stan版

目次

ブログの紹介

note で８つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．統計・データ分析とつながる

シリーズ「統計・データ分析とつながる」は、統計・データ分析との「つながり」を発掘して、コラム風に仕立てたブログシリーズです。
生成 AI の力を借りながら、統計・データ分析の入り口をイメージして、自由気ままに書きました。
たとえば…
・日常生活と統計のつながり
・統計検定２級からその先へのつながり
気楽にお読みいただけたら嬉しいです🍀

３．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

４．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

５．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

６．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

７．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

８．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

いいなと思ったら応援しよう！

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します！

この記事が参加している募集

#

最近の学び

31万件