AutoEncoderを用いた政令指定都市の区の分析と, 主成分分析との比較
目的
AutoEncoderにより、東京23区および政令都市の区(n=198)のデータの次元削減を行い, 主成分分析と比較する.
gakushukun1.hatenablog.com
前回の記事で主成分分析を行なった各政令指定都市の区のデータについて, 今回はAutoEncoderにより次元削減を行う.
具体的には, 各区のデータが持つ11次元の情報について, 下のような構成のエンコーダー・デコーダーを組み合わせたニューラルネットにより, 2次元の情報に落とし込む.
Kerasにより実装したコードは次の通り. なお, 分析に用いたベクトルは, 主成分分析の時と同様, 地域ごとの人口差による影響を除くために, 各次元を総人口で割っている. また今回は, AutoEncoderの出力層の活性化関数をシグモイド関数としたため, ベクトルの各次元の値が0〜1の間に収まるよう正規化も行なった.
import numpy as np from keras.models import Sequential from keras.layers import Dense, Activation, BatchNormalization from keras.activations import sigmoid from keras.optimizers import Adam from keras.utils import plot_model import matplotlib.pyplot as plt import matplotlib import pandas as pd if __name__ == '__main__': matplotlib.rcParams['font.family'] = 'IPAexGothic' df = pd.read_csv('ku2015.csv', encoding="SHIFT-JIS") df.loc[:, '総人口':] = df.loc[:, '総人口':].astype(float) pop = df.loc[:, '総人口'].values X = df.iloc[:, 2:].apply(lambda x: x / pop, axis=0) X_norm = X.loc[:, :].apply(lambda x: x / x.max(), axis=0) print(len(X_norm.columns)) # エンコーダー encoder = Sequential() encoder.add(Dense(6, input_shape=(11,))) encoder.add(BatchNormalization()) encoder.add(Activation(sigmoid)) encoder.add(Dense(2)) encoder.add(BatchNormalization()) encoder.add(Activation(sigmoid)) # デコーダー decoder = Sequential() decoder.add(Dense(6, input_shape=(2,))) decoder.add(BatchNormalization()) decoder.add(Activation(sigmoid)) decoder.add(Dense(11)) decoder.add(BatchNormalization()) decoder.add(Activation(sigmoid)) model = Sequential([encoder, decoder]) model.compile(optimizer=Adam(), loss='mean_squared_error', metrics=['accuracy']) fit = model.fit(X_norm, X_norm, epochs=2000) y = encoder.predict(X_norm)
上記のコードによりAutoEncoderの学習を行なった後, 元のデータ(コード中のX_norm)をEncoderに入力して得られた出力が, 今回必要な次元削減された情報である.
図の横軸は, Encoderから得られた出力の一次元目を, 縦軸は二次元目を表している. 主成分分析の主成分得点の結果と比較すると, よく似ているが, 千代田区, 中央区, 港区については, 曲がった座標に埋め込まれていることが分かる. この３つの区は, 特別な性質を持つのかもしれない.
上の図の各x=(0〜0.75), y=(0〜0.75)の点について, Decoderに入力して得られた出力を次に示す.
x, y=(0, 0)から復元された情報では, 5:転入, 6:転出, 9:単独世帯,10:婚姻件数,11:離婚件数が少なく, 人の往来の少ない町であると分かる. x, y=(0.75, 0)は, 4:死亡数, 9:単独世帯が多く, 高齢化の進んだ町であると分かる. x, y=(0, 0.75)は, 4:死亡者数が少なく, 8:昼間人口が多いためビジネス街のような場所であると推測される. x,y=(0.75, 0.75)は, 5:転入, 6:転出, 9:単独世帯,10:婚姻件数,11:離婚件数が多いため, 人の往来が激しい町であると分かる.
まとめ
AutoEncoderを用い政令都市の区のデータの分析を行い, 主成分分析と比較した. AutoEncoderでは, 主成分分析とほぼ同じ結果が得られたが, いくつかの町では, 線形変換ではみつからない特徴が取り出されていた.
主成分分析を用いた政令指定都市の区の分析
目的
主成分分析を用いて、東京23区および政令都市の区()の分析を行う.
用いたデータは, e-statの2015年の調査によるもので, 総人口,１５歳未満人口,１５〜６４歳人口, 出生数, 死亡数, 転入者数, 転出者数, 従業・通学していない人口, 昼間人口, 単独世帯, 婚姻件数, 離婚件数の
12次元のベクトルである. 今回は, このベクトルに対して主成分分析を行い, 各因子の寄与率・主成分ベクトル・主成分得点を計算することで, 区の特徴をどういったもので定められているかを分析する.
pythonで実装したコードは次の通りである. なお, 実際に主成分分析を行ったベクトルは, もとの各次元のうち総人口以外を総人口で割った11次元のベクトルである(地域ごとの人口差による影響を除くため). また, 主成分分析を行う前に各ベクトルについて標準化を行なっている.
import pandas as pd import numpy as np from sklearn.decomposition import PCA if __name__ == '__main__': matplotlib.rcParams['font.family'] = 'IPAexGothic' df = pd.read_csv('ku2015.csv', encoding="SHIFT-JIS") print(df) df.loc[:, '総人口':] = df.loc[:, '総人口':].astype(float) pop = df.loc[:, '総人口'].values X = df.iloc[:, 2:].apply(lambda x: x / pop, axis=0) X_norm = X.loc[:, :].apply(lambda x: ((x - x.mean()) / x.std()), axis=0) for i in range(len(X_norm.columns)): print(i + 1, end=' : ') print(X_norm.columns[i]) pca = PCA() pca.fit(X_norm) prop = pca.explained_variance_ / np.sum(pca.explained_variance_)
コード中の変数propに格納されている各主成分の寄与率を図で示す (寄与率は合計で1になるよう正規化した).
図より, 第二主成分までの寄与率の和が0.8であり, 全体の情報のうちの大部分がここに集中していることが分かった.
これに対して, 第一主成分に対応する主成分ベクトルは次のようになった.
このベクトルは, 労働人口(15~64歳人口)・転入出者数・婚姻件数が多いほど+の影響があり, 死亡者数・従業通学していない者の数が多いほど-の影響がある.
このことから, このベクトルの+の方角は区が人の往来が多く賑やかであり, -の方角は, 落ち着いたのどかであることが考えられる.
一方, 第二主成分に対応する主成分ベクトルは次のようになった.
このベクトルは, 死亡者数・単独世帯数が多いほど+の影響があり, 15歳未満人口・出生数が多いほど-の影響がある.
このことから, このベクトルの+の方角に行くほど少子高齢化が進んでいる区であると考えられる.
この２つの因子に対し, 東京23区の各区について, 第一主成分得点をx軸に, 第二主成分得点をy軸にプロットした結果を次に示す.
この図から, 東京23区は, 日本の他の政令指定都市と比べると, 人の往来が多く, また少子高齢化については平均的であることが分かる.
まとめ
主成分分析により, 東京23区と政令都市のデータの分析を行った. これらのデータは主に, ２つの主成分によって特徴づけれらていることが分かった.
ADVI(PyStan)を用いて混合正規分布の推測
目的
混合正規分布の推測をStanのADVI(Automatic Differentiation Variational Inference)を用いて行う.
混合正規分布の事後分布は, StanのNUTS(No-U-Turn-Sampling)ではうまく作れないことが多いので, 代わりに変分ベイズ法のADVI(Automatic Differentiation Variational Inference)を用いて作る.
データは, 真のパラメータが既知である3個の正規分布の混合から発生する(サンプルサイズ). ADVIでは, 推定されるパラメータが揺らいでいるので, 1000回得られたパラメータサンプルのうち, 後半の500個を平均して, モデルにplug-inすることで, 予測分布を構成した.
実験に使ったコードを次に示す.
import pystan import numpy as np from scipy.stats import multivariate_normal, norm import pandas as pd import matplotlib.pyplot as plt if __name__ == '__main__': np.random.seed(seed=123) N = 600 # サンプル数 E = [[1.0, 0.0], [0.0, 1.0]] # 単位行列 t_a = np.array([1/3]*3) # 真の混合比 # t_mu = np.array([[0.0, 2.0], [-np.sqrt(3), -1.0], [np.sqrt(3), -1.0]]) t_mu = np.array([[0.0, 1.5], [-np.sqrt(3) * 2 / 3, -0.75], [np.sqrt(3) *2 / 3, -0.75]]) # 真の平均値 t_sigma = np.array([E, E, E]) # 真の分散共分散(単位行列) t_z = [np.argmax(np.random.multinomial(1, t_a)) for i in range(N)] # 真の潜在変数 # 入力するサンプル Y = [np.random.multivariate_normal(t_mu[t_z[i]], t_sigma[t_z[i]]) for i in range(N)] # Stanに入力 data = {'N': N, 'K': 3, 'M': 2, 'Y': Y} sm = pystan.StanModel(file='mixnorm_.stan') fit = sm.vb(data=data, seed=123) # 結果の取得 sample = pd.read_csv(fit['args']['sample_file'].decode('utf-8'), comment='#') sample = sample.drop([0, 1]) # 予測されたパラメータ predict_a = [np.mean(sample['a.1'][501:1000]), np.mean(sample['a.2'][501:1000]), np.mean(sample['a.3'][501:1000])] predict_mu = [[np.mean(sample['mu.1.1'][501:1000]), np.mean(sample['mu.1.2'][501:1000])], [np.mean(sample['mu.2.1'][501:1000]), np.mean(sample['mu.2.2'][501:1000])], [np.mean(sample['mu.3.1'][501:1000]), np.mean(sample['mu.3.2'][501:1000])]] predict_sigma = np.mean(sample['sigma'][501:1000]) x, y = np.meshgrid(np.arange(-5.0, 5.0, 0.05), np.arange(-5.0, 5.0, 0.05)) pos = np.dstack((x, y)) gauss = [multivariate_normal(mean=t_mu[k], cov=t_sigma[k]) for k in range(3)] z = sum([t_a[k] * gauss[k].pdf(pos) for k in range(3)]) color = ['red', 'blue', 'black'] predict_gauss = [[norm(loc=predict_mu[k][m], scale=predict_sigma) for m in range(2)] for k in range(3)] predict_z = [[sum([predict_a[k] * np.prod([predict_gauss[k][m].pdf(pos[i, j, m]) for m in range(2)]) for k in range(3)]) for j in range(200)] for i in range(200)] # プロット plt.figure() plt.subplot(2, 1, 1) plt.contour(x, y, z) for i in range(N): plt.scatter(Y[i][0], Y[i][1], c=color[t_z[i]], s=10) plt.subplot(2, 1, 2) plt.contour(x, y, predict_z) plt.show()
- Stan
data {
int N;
int K;
int M;
vector[M] Y[N];
}
parameters {
simplex[K] a;
vector[M] mu[K];
real<lower=0> sigma;
}
transformed parameters {
vector[K] lp[N];
for (n in 1:N) {
for (k in 1:K) {
lp[n, k] = log(a[k]);
for (m in 1:M) {
lp[n, k] += normal_lpdf(Y[n, m] | mu[k, m], sigma);
}
}
}
}
model {
a ~ dirichlet(rep_vector(5.0, K));
for (k in 1:K) {
for (m in 1:M) {
mu[k, m] ~ normal(0.0, 5);
}
}
for (n in 1:N) {
target += log_sum_exp(lp[n]);
}
}
generated quantities {
simplex[K] pr_z[N];
for(n in 1:N) {
pr_z[n] = softmax(lp[n]);
}
}
真の分布は, 二次元のユークリッド空間上に定義され, サンプル点は二次元の点からなる.
真の分布と, サンプル点は, 次の図の通りである. 各点の色は, どの正規分布から発生したかを表しているが, 学習モデルにはその情報は与えられない(潜在変数).
このサンプルに対して, ADVIを用いて作った予測分布を次に示す. 良好に推定されていることが分かる.
次に, 3個の山が重なりを持つ場合について実験した. 真の分布とサンプルは次の通りである.
このサンプルに対して, ADVIを用いて作った予測分布を次に示す. 予測分布の右下が他と比べてやや大きくなっているが, これは発生したサンプルの揺らぎが原因であると思われ, 推測自体は良好にできていると考えられる.
実行中に,
WARNING:pystan:Automatic Differentiation Variational Inference (ADVI) is an EXPERIMENTAL ALGORITHM.
という警告が発生するが, 今回は良好に推測できていた. ただし, ここでは全ての正規分布の分散共分散行列を共通の単位行列の定数倍というモデルを用いており, コンポーネントごとに異なる分散共分散を持つモデルでは適切な推測は困難であった.
そもそも, 異なる分散共分散行列を持つ混合正規分布の推測は, 他の方法によっても困難であるので, StanのADVIは有力な方法であると考えられる.
まとめ
StanのADVIを用いて混合正規分布の学習を行った.
StanのADVIは開発中とのことであるが, 一般にNUTSは混合分布には適さないことが多いように思われるので, 今の所はADVIを使うほうがよいかもしれない.
階層ベイズを用いた地域ごとの家賃の変化の推測
目的
RStanによる階層ベイズを用いて, 日本の各地域ごとの家賃の変化を推測する.
今回は, e-statの日本の様々な町ごとの1991年と2000年の家賃のデータを用いて, 各地域の家賃の変化の推測を行う. その際, 日本全体で大まかな変化の傾向があることを仮定して, 階層ベイズモデルを適用する.
具体的には, データに含まれる67箇所の都市を6つの地域(北海道&東北, 関東, 中部, 関西, 中国&四国, 九州)に分類して, 次のようなモデルにより線形回帰を行う.
(ただしはサンプル点
の所属する地域を表す)
このモデルでは, 地域のX(1991年の家賃)とY(2000年の家賃)の関係を示す線形モデルのパラメータ
が, 平均
, 標準偏差
の正規分布から生成されていると考える(
も同様に平均
, 標準偏差
の正規分布から生成されているとする). また, 標準偏差
は,
ごとに定数1.0と
をパラメータとしたガンマ分布から生成しているとする. こうした階層的なモデルを用いることで, それぞれの地域ごとの少ないサンプルで線形回帰した場合に起こる過学習を防げる可能性がある. また. 他の地方の情報が, 事前分布を通じて, 別地方の推定に活用できる.
上記のモデルをStanで定義し, 各パラメータの事後分布のサンプリングを行うコードは次の通りである.
- R
library(rstan) options(mc.cores = parallel::detectCores()) d <- read.csv("yachin.csv") area <- d$X[2:68] town <- d$X.1[2:68] rent1991 <- d$X.2[2:68] / 1000 rent2000 <- d$X.3[2:68] / 1000 df <- data.frame(town, rent1991, rent2000, area, stringsAsFactors=FALSE) data <- list(N=nrow(df), K=6, X=df$rent1991, Y=df$rent2000, KID=df$area) fit <- stan(file="kaisou.stan", data=data, seed=123)
data {
int N;
int K;
real X[N];
real Y[N];
int<lower=1, upper=K> KID[N];
}
parameters {
real a0;
real b0;
real<lower=0> beta0;
real a[K];
real b[K];
real<lower=0> s_a;
real<lower=0> s_b;
real<lower=0> s_Y[K];
}
model {
for (k in 1:K) {
a[k] ~ normal(a0, s_a);
b[k] ~ normal(b0, s_b);
s_Y[k] ~ gamma(1.0, beta0);
}
for (n in 1:N) {
Y[n] ~ normal(a[KID[n]] + b[KID[n]] * X[n], s_Y[KID[n]]);
}
}
generated quantities {
vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = normal_lpdf(Y[n] | a[KID[n]] + b[KID[n]] * X[n], s_Y[KID[n]]);
}
}なお, 元データの家賃は民営借家一ヶ月(あたり)の値であり, 前処理の段階で扱いやすくするために1/1000倍している.
これとは別に, 各地域ごとに別々に線形回帰を行った. 地域ごとに行った回帰により予測された線を黒, 階層ベイズを用いた回帰により予測された線を赤でプロットすると, 次の図ようになった.
(グラフ中の番号は , 1...北海道&東北, 2...関東, 3...中部, 4...関西, 5...中国&四国, 6...九州に該当, また点は各地方に属している都市のデータを表す)
グラフを見ると, 2(関東)のように２つのモデルで予測された線がほぼ同一になるものもあれば, 5(中国&四国)のようにまったく異なる線になる場合があると分かる. 特に5の場合, 家賃が同じ価格帯の都市のデータばかりであったために, 地域単体の回帰では, 日本全体の傾向(傾きが1.00に近い=1991年から2000年にかけて家賃が横ばい)から外れたものになっていると考えられる.
地域ごとのデータで行った回帰と, 階層ベイズによる回帰の汎化性能を比べるため, それぞれの場合での予測分布に対するWAICとLOOの値を次に示す. 値が小さい方が汎化性能が優れているといえる.
| WAIC(個々) | WAIC(階層) | LOO(個々) | LOO(階層) | |
|---|---|---|---|---|
| 1.北海道&東北 | 11.9 | 10.1 | 12.2 | 10.2 |
| 2.関東 | 25.9 | 24.1 | 26.1 | 24.2 |
| 3.中部 | 29.7 | 26.7 | 30.6 | 27.1 |
| 4.関西 | 22.4 | 20.9 | 22.8 | 21.5 |
| 5.中国&四国 | 12.4 | 14.3 | 12.8 | 14.6 |
| 6.九州 | 15.7 | 12.8 | 15.9 | 13.0 |
5.の中国&四国を除いて, 階層ベイズによって予測された回帰の方が優れていると分かる.
RStanで重回帰を試してみる
目的
e-statの統計データから, rstanによる重回帰を用いて年代ごとの趣味の関係を調べる.
今回はRStanを使って重回帰分析を試してみる. 対象としたのは, e-statの中の,
社会生活基本調査 平成28年社会生活基本調査 調査票Ａに基づく結果 生活行動に関する結果 生活行動編（全国） 趣味・娯楽 である. このデータでは, おおまかな世代ごとの趣味の傾向について調査の結果が掲載されている. 例として, 10〜14歳の世代のデータでは, 調査を行ったサンプル9351人のうち, 趣味にスポーツ観覧が含まれると答えた人が1804人, 美術鑑賞が含まれていると答えた人が854人ということが分かる(なお, これらは重複解答を含んでいる).
本エントリでは, 複数ある趣味のカテゴリの項目のうち, スポーツ観覧, 園芸, とテレビ・パソコンゲームの関係を調べてみることとする. 使用するモデルは次の通り.
このモデルを用いて, をスポーツ観覧が趣味の人の割合,
を園芸が趣味の人の割合,
をテレビ・パソコンゲームが趣味の人の割合として回帰を行う. ここでいう割合とは, 調査を行った各世代ごとのサンプルサイズに対して, それぞれの趣味を持っていると回答した人の割合である. なお, 回帰に用いるデータは10〜14歳の世代, 15～19歳の世代, …, 85歳以上の世代の16個である.
重回帰分析に用いたコードを次に示す.
- R
library(rstan) library(maptools) d <- read.csv("FEH_00200533_190529214744.csv", fileEncoding="shift-jis", stringsAsFactors=FALSE) age_str <- strsplit(unlist(d[2:17,6]), "_") age <- c() for(i in 1:length(age_str)) { age <- append(age, age_str[[i]][2]) } sample_size <- as.numeric(unlist(d[2:17,8])) sport <- as.numeric(unlist(d[2:17,14])) / sample_size garden <- as.numeric(unlist(d[2:17,54])) / sample_size game <- as.numeric(unlist(d[2:17,74])) / sample_size df <- data.frame(age, sport, garden, game) data <- list(N=nrow(df), X=df$sport, Y=df$garden, Z=df$game) fit <- stan(file="jukaiki.stan", data=data, seed=123) ms <- extract(fit) x <- seq(0.0, 0.3, length=50) y <- seq(0.0, 0.3, length=50) func <- function(x, y) {mean(ms$a) + mean(ms$b) * x + mean(ms$c) * y} z <- outer(x, y, func) png("hobby.png", width=2000, height=2000, res=300 ) par(family = "HiraKakuProN-W3") contour(x, y, z) par(new=T) plot(df$sport, df$garden, xlab="watching_sports", ylab="gardening", xlim=c(0.0,0.3), ylim=c(0.0,0.3)) pointLabel(x=df$sport, y=df$garden, labels=df$age) dev.off()
- Stan "jukaiki.stan"
data {
int N;
real X[N];
real Y[N];
real Z[N];
}
parameters {
real a;
real b;
real c;
real<lower=0> sigma;
}
model {
for (n in 1:N) {
Z[n] ~ normal(a + b * X[n] + c * Y[n], sigma);
}
}
generated quantities {
vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = normal_lpdf(Z[n] | a + b * X[n] + c * Y[n], sigma);
}
}各パラメータのサンプリング結果の平均は次のようになった.
| a | b | c | sigma |
|---|---|---|---|
| 0.07 | 2.08 | -0.78 | 0.04 |
これらの値を元のモデルに当てはめ, テレビ・パソコンゲームを等高線でプロットした図は次のようになる. なお, 図中の点は分析に用いた元のデータをプロットしたものである. サンプルサイズが小さいものの, 良好に推測できていると考えられる.
これより, スポーツを観戦する人ほど, また園芸をしない人ほどテレビ・パソコンゲームをすることが分かる. また, 図の右下から左上に向かって年齢層が上がっていくことが分かる.
推測により得られた重回帰平面を三次元に図示したものが次である.
まとめ
RStanを用いて重回帰分析を行い, スポーツ観戦と園芸から, テレビ・パソコンゲームへの関係を調べた.
RStanで単回帰を試してみる
目的
estatsの統計データから, rstanによる単回帰を用いて各都道府県の死亡者数と出生児数の関係を調べる. 事後分布からサンプリングされたパラメータから, 情報量規準WAICと交差検証LOOを計算し, 単回帰と二次回帰どちらが真の分布に近いかを推測する.
今回はrstanを使って単回帰分析を試してみる. 今回対象としたのは, estatの中の, 平成29年10月〜平成30年9月間の都道府県別の出生児数・死亡者数のデータである.
このデータについて, 死亡者数を横軸, 出生児数を縦軸で表すと次のようなグラフとなる. それぞれの点が各都道府県のデータを表している.
このデータから, 死亡者数から出生児数への回帰曲線を推測する. 死亡者数を, 出生児数を
とし, 次のような2つのモデルを考える.
それぞれのモデルについてstanでパラメータの事後分布を作り, サンプリングを行った. 元の確率モデルに対して. サンプリングで得られたパラメータを当てはめた数値を保存しておくことで, 擬似的な予測分布を作り, waicとlooの計算を行った. 実際のコードは次の通り.
- 線形回帰(R)
library(rstan) library(loo) d <- read.csv("FEH_00200524_190526235535.csv", stringsAsFactors=FALSE) d[,2] <- d[,2] / 10000 d[,3] <- d[,3] / 10000 x <- d$死亡者数 y <- d$出生児数 data <- list(N=nrow(d), X=x, Y=y) fit <- stan(file="kaiki.stan", data=data, seed=123) log_lik <- extract_log_lik(fit) ms <- extract(fit) a <- mean(ms$a) b <- mean(ms$b) sigma <- mean(ms$sigma) waic01 <- waic(log_lik) loo01 <- loo(log_lik) png("kaiki1.png", width=800, height=800) par(family = "HiraKakuProN-W3") plot(x, y, xlab="死亡者数", ylab="出生児数", xlim=c(0,13), ylim=c(0,13)) abline(a, b) abline(a + sigma * 1.96, b) abline(a - sigma * 1.96, b) dev.off()
- 線形回帰(Stan) "kaiki.stan"
data {
int N;
real X[N];
real Y[N];
}
parameters {
real a;
real b;
real<lower=0> sigma;
}
model {
for (n in 1:N) {
Y[n] ~ normal(a + b * X[n], sigma);
}
}
generated quantities {
vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = normal_lpdf(Y[n] | a + b * X[n], sigma);
}
}
- 二次回帰(R)
library(rstan) library(loo) d <- read.csv("FEH_00200524_190526235535.csv", stringsAsFactors=FALSE) d[,2] <- d[,2] / 10000 d[,3] <- d[,3] / 10000 x <- d$死亡者数 y <- d$出生児数 data <- list(N=nrow(d), X=x, Y=y) fit <- stan(file="kaiki2.stan", data=data, seed=123) log_lik <- extract_log_lik(fit) ms <- extract(fit) a <- mean(ms$a) b <- mean(ms$b) c <- mean(ms$c) sigma <- mean(ms$sigma) fun <- function(x) { a + b * x + c * (x ^ 2) } fun_upr <- function(x) { a + b * x + c * (x ^ 2) + sigma * 1.96 } fun_lwr <- function(x) { a + b * x + c * (x ^ 2) - sigma * 1.96 } waic01 <- waic(log_lik) loo01 <- loo(log_lik) png("kaiki2.png", width=800, height=800) par(family = "HiraKakuProN-W3") plot(x, y, xlab="死亡者数", ylab="出生児数", xlim=c(0, 13), ylim=c(0, 13)) par(new=T) plot(fun, xlab="", ylab="", xlim=c(0, 13), ylim=c(0, 13)) par(new=T) plot(fun_upr, xlab="", ylab="", xlim=c(0, 13), ylim=c(0, 13)) par(new=T) plot(fun_lwr, xlab="", ylab="", xlim=c(0, 13), ylim=c(0, 13)) dev.off()
- 二次回帰(Stan)
data {
int N;
real X[N];
real Y[N];
}
parameters {
real a;
real b;
real c;
real<lower=0> sigma;
}
model {
for (n in 1:N) {
Y[n] ~ normal(a + b * X[n] + c * (X[n] ^ 2), sigma);
}
}
generated quantities {
vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = normal_lpdf(Y[n] | a + b * X[n] + c * (X[n] ^ 2), sigma);
}
}元のデータの数値が大きいため, 線形, 二次回帰いずれの場合も全てのデータを10000分の1倍している. 予測分布のYについての平均を計算することで, 回帰曲線を推測し, サンプリングされたノイズの標準偏差の1.96倍の値を足し引きすることで, 与えられた
に対する
の95%のベイズ予測区間を示した.
プロットした結果を次に示す. グラフ中の真ん中の線が推定された回帰曲線で, その上下の線が95%ベイズ予測区間の上限と下限である.
- 線形回帰
- 二次回帰
95%ベイズ予測区間の外にある都道府県は, が小さい方から, 沖縄県, 北海道, 愛知県であった.
また, 得られたWAICとLOOの値は次のようになった.
| WAIC | LOO | |
|---|---|---|
| 線形回帰 | 62.7 | 64.3 |
| 二次回帰 | 44.4 | 45.0 |
これより, 二次回帰の方が線形回帰よりも真の予測分布に近いと考えられる.
一般に, 死亡者数と出生児数は比例すると思われるが, 日本の都道府県データでは, 死亡者数が多いほど出生児数の割合も大きくなることが分かる. これは, 人口が多い都市圏(東京など)ほど, 高齢者数に比べて若年層が多いということを意味していると考えられる.
まとめ
rstanを用いて, 線形回帰と二次回帰を行い, 死亡者数から出生児数への回帰曲線を推測した. 情報量規準, 交差検証を用いて, 二次回帰の方がより真の分布に近いと推定された.
