顧客１人あたりの売上やら何やら、本当に様々なデータが裾の長い度数分布の形状を持ちます。
そういったデータに対して何か説明変数を設定してモデリングをする必要に迫られるシーンを想定して、下図のごときデータにフィットするモデリングを考えてみます。本記事では目的変数をyとおき、説明変数(単変量)をxとおきます。

f:id:sakatoken:20190611000403p:plain

ダミーデータ生成など

ダミーデータ生成とライブラリ読み込みをしておきます。

library(tidyverse)
library(GGally)
library(rstan)
library(tidybayes)

d <- tibble(
    x = abs(rnorm(2000, mean =1.5))
) %>%
    mutate(x = x / max(x)) %>%
    mutate(
        mu = -4 * x,
        sigma = exp(-0.5*x -0.5)
    ) %>%
    mutate( 
        y = as.numeric(map2(mu, sigma, function(m, s){rlnorm(1, meanlog = m, sdlog = s)}))
    )

分布を仮定する

この問題を考えるには「裾の長い非負の度数分布を見たときにどのような分布を仮定するべきか？」に明快に答えを出したいところですが、今回は大変安直ながら指数分布と対数正規分布を見比べて

xに応じて平均と分散の両方をモデリングできる
xに応じて分布の峰が0から離れる分布のほうが、ビジネスの実問題にはフィッティングさせやすいと思われる

という理由で対数正規分布を選びました。

※ダミーデータでは対数正規分布が真のモデルになっています。

モデルの記述

xと対数正規分布のパラメータの間には次の仮定を設定しました。

対数正規分布のパラメータには平均と分散があります。今回はそれらの両方に対し、xによる線形予測子を仮定します。
問題によりますが、ここではxは実数全体を取るため、分散（非負）と結びつけるにあたって指数を使います。
線形予測子に含む係数部分には無情報事前分布を仮定します。

stanコードは次の通りです。単変量であるはずのxにmatrixを指定しているのは、Rコード側で切片項を追加するためです。

data {
    int N;
    int D;
    vector<lower=0>[N] Y;
    matrix<lower=0, upper=1>[N, D] X;
}

parameters {
    vector[D] b;
    vector[D] a;
}

transformed parameters {
    vector<lower=0>[N] s;
    vector[N] mu;
    
    s = exp(X * b);
    mu = X * a;
}

model {
    Y ~ lognormal(mu, s);
}

前処理とMCMC実行

以上でモデルを定義できたので、前処理してMCMCを実行する部分のRコードを記述していきます。

# 外れ値 除去
d.f <- d %>%
    filter(
        x < quantile(x, .98), 
        y < quantile(y, .98)
    )

# 最大値を使ってスケーリング
d.f.max <- d.f %>%
    summarize_all(max) %>%
    select(x,y)
d.f <- d.f %>%
    mutate(
        x = x / d.f.max$x,
        y = y / d.f.max$y
    )

# stanにわたす用のlist型データを作成
d.stan <- list(
    N = nrow(d.f),
    Y = d.f$y,
    X = d.f %>%
        select(x) %>%
        mutate(intercept = 1) %>% # 切片項を追加
        as.matrix(),
    D = 2
)

# MCMCを実行
fit <- stan(file = './modeling.stan', data = d.stan, seed = 0,
            iter = 700, warmup = 200, chains = 3, pars=c('b', 'a'))

# 収束判定にはtraceplot とRhat値のチェックをするけどコード略。

モデルの実データへのフィット具合をグラフで確認

以上で推定したモデルを使ってyの予測値を新たに生成し、これが実測のyを言い当てるかどうかを確認します。ただし冒頭の図で確認した通り、今回は特定のxの値に対してyは広くばらつく分布を持っているため、予測値と実測値の誤差の小ささという基準でフィットを評価しようとすると「まるでデタラメなモデルである」という結果になります。

そこで「モデルは個別のyを言い当てることはできないが、xがある値に固定されたときのyの分布を言い当てることが出来る」という立場を取り、グラフを作成することとします。

具体的には以下の手順を取ります。

推定されたパラメータを使って、xからμとσを算出する。
μとσにしたがって、新しいy_predを生成する。
μの値でデータを等間隔の階級に区切る。
階級別にyとy_predのヒストグラムを重ねてプロットする。
グラフでy_predが十分にyに追従できているのかどうかを確認する。

以下、コードなどです。

まずは先程のMCMCによる推定結果から、パラメータの値を控えておきます。

Inference for Stan model: modeling.
3 chains, each with iter=700; warmup=200; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=1500.

       mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
b[1]  -0.30    0.00 0.07  -0.43  -0.34  -0.30  -0.25  -0.15   930    1
b[2]  -0.53    0.00 0.03  -0.59  -0.55  -0.53  -0.50  -0.46   906    1
a[1]  -2.71    0.00 0.05  -2.81  -2.75  -2.71  -2.68  -2.61   592    1
a[2]  -0.65    0.00 0.03  -0.70  -0.66  -0.65  -0.63  -0.59   657    1
lp__ 298.37    0.06 1.40 294.74 297.71 298.69 299.39 300.09   619    1

Samples were drawn using NUTS(diag_e) at Thu Jun  6 18:01:56 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

beta = c(-0.3, -0.53)
alpha = c(-2.71, -0.65)

続いて、グラフをプロットするまでのRコードです。

# μとσにしたがって、新しいy_predを生成する。
d.f.pred <- d.f %>%
    mutate(
        sigma = exp(x * beta[1] + beta[2]),
        mu = x*alpha[1] + alpha[2]
    ) %>%
    mutate(
        y_pred = rlnorm(n(), meanlog = mu, sdlog = sigma)
    ) %>%
    filter(y_pred < quantile(y_pred, .98)) # 外れ値 除去

# 階級別にyとy_predのヒストグラムを重ねてプロットする。
d.f.pred %>%
    mutate(class = cut_interval(mu, n = 8)) %>% # μの値でデータを等間隔の階級に区切る。
    select(y, y_pred, class) %>%
    gather("name", "value", y, y_pred) %>% 
    arrange(class) %>%
    mutate(name = if_else(name == "y", "実測", "予測")) %>%
    ggplot(
        aes(value, fill = name)
    ) +
    geom_histogram(bins = 15, position = "identity", alpha = 0.5) +
    facet_wrap(~class, ncol = 2, scales="free") +
    labs(
        x = "推定μ(対数正規分布)",
        fill = "予測と実測の色分け",
        title="推定μの階級別ヒストグラム"
    )