片方しかない靴下の数をゼロ過剰ポアソン分布でモデリングしてみた


ShoesZIPhist.png ZIPFullResult.png
author: Mr.Unadon (見習い飯炊き兵)
動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1



結論


【今回、僕が仮定した考え方】

「世の中には、片方の靴下をなくすタイプのひとと、なくさないタイプがいる」

「靴下をなくす人の中には、のこった片方を捨ててしまうなどで、片方だけの靴下を一足も持っていない人がいる」

「靴下をなくす人が持っている片方の靴下の数には個人差があり、ポアソン分布する」

【結果】

「靴下をなくす人は世の中の82%くらいだろう」

「靴下をなくした人の中では、 もう片方を捨ててしまうなどで持ってないひともおり、平均的に1.2足くらい、片方だけの靴下を所持している」

「持っている靴下の数は個人で違っており、ひとそれぞれ」


はじめに


本稿は『StanとRでベイズ統計モデリング』の勉強会での発表の一部です。

11章で出てきましたゼロ過剰ポアソン分布について、実データで楽しく学んでいただけたら嬉しいなと思い、整理したものとなります。

なお、詳細は上記書籍を御覧ください。この記事よりも間違いなく理解が進みます。


ところで、みなさま。

靴下って、なぜか片方だけなくなりませんか?

酔っ払って上着や財布を失うことはあっても、靴下だけは両方揃って履いて帰ってくる。

なのに、片方だけない靴下を、僕は6足持っていました。

ほんと不思議です。僕の靴下たち…どこいったんやろか…

不思議で不思議でたまらなかったので、友人に「靴下って片方なくなるよね」と話しました。

そこで、「え?」という反応が帰ってきたとき、私は生まれて初めて知ったのです。

世の中には、片方の靴下をなくさない人もいるんだ!

そして同時に、

片方しか無い靴下の所持数はゼロ過剰ポアソン分布に違いない!

とも、感づきました。

ゼロ過剰ポアソン分布とはなんぞや。

靴下片方をなくさない人がいて、なくす人もいる。

「ある人が、靴下なくす人である確率」はどれくらいだろう?

片方なくすタイプの人は、「平均何足」くらい、なくしてるんだろう?

そんな疑問に答えることのできる、統計学の一つのモデルです。

一般の方 + 統計にお詳しい方々の両方がご覧になられる可能性がありますので、 少し難しく、あるいは少し物足りなく感じられる内容になっていくこと、どうかご容赦くださいませ。


1. 「片方しか無い靴下の所持数」データの概要


先日、Twitterでアンケートを取らせて頂きました。

「あなたは、片方しかない靴下を何足もっていますか?」

ShoesZIPTweet.png

ご回答くださった方々は合計1304名(!!!)。

リツイートにより拡散してくださった方々は、(重複あり)125名様でした。

みなさま、本当にありがとうございます。

なお、今回のアンケートの最大値は7足でした。7足以上持っている方は、7に投票していただいたと思います。

実際に、「7足以上持ってるよ」というコメントもありました。こういうデータを、「打ち切りデータ」といいます。

ここも後々考慮していきますが、

まずはデータを見てみましょう。

データ解析に用いるR言語というプログラミング言語で、

片方しか無い靴下の数ごとに、 その人数をグラフ化してみます。

invisible({rm(list=ls());gc();gc()})
###########################################
# 片方しか無い靴下の所持数を 
# Zero-Inflated Poisson Distributionで
#   Author: @MrUnadon
###########################################
library(rstan)
library(tidyverse)
library(formattable)

#stan parallel
rstan_options(auto_write=TRUE)
options(mc.cores = parallel::detectCores())

#Data
zero<-rep(0,746)
one<-rep(1,142)
two<-rep(2,142)
thr<-rep(3,154)
fou<-rep(4,23)
fiv<-rep(5,22)
six<-rep(6,9)
sev<-rep(7,66)
dat<-data.frame(Y=c(zero,one,two,thr,fou,fiv,six,sev))
    

#Histogram
ggplot()+theme_set(theme_bw(base_size = 18,base_family = "HiraKakuProN-W3"))
gg1<-ggplot(dat,aes(x=Y))+
    geom_histogram(position="identity",fill="black",colour="black",alpha=0.5,bins=8)+
    theme(plot.subtitle = element_text(vjust = 1), 
    plot.caption = element_text(vjust = 1), 
    axis.title = element_text(family = "serif", 
        size = 18), axis.text = element_text(family = "serif", 
        size = 14), plot.title = element_text(family = "serif", 
        size = 22)) +labs(title = "片方しか無い靴下の数の分布", 
    x = "片方しか無い靴下の個数", y = "人数")

#描画
plot(gg1)

いぃぃよいしょっと。

ShoesZIPhist.png
Figure1 【度数分布】片方しか無い靴下の所持数(回答者数1304名)

このグラフですが、縦軸が人数、横軸が片方しか無い靴下の所持数です。

0足という方々が多いという特徴を持っていますね。

また、7足以上の方も一定数いらっしゃったようで、その方々のデータが7に含まれて、少し増えていることがわかります。


2. ゼロ過剰ポアソン分布について


ゼロ過剰ポアソン分布は、「靴下をなくす確率」と「なくした靴下の平均的な数」を推定するのに便利な統計学の道具です。

靴下をなくすという事象の生起確率と、無くした数という数えられる値の平均値を、 統計理論に基づいて考えることができます。

専門的に言えば、ゼロ過剰ポアソン分布は、ベルヌイ分布とポアソン分布の混合分布。

下記の図で、およそのイメージを掴んでいただけたらと思います。

ShoesZIPExplain1.png

今回は、単純に無くした靴下の数の平均値や、無くした靴下の数が1足以上の方の割合を求めるだけにとどめません。

もう一歩踏み込んで、「靴下をなくすタイプな人がどれくらいいるか」、「平均的に何足くらいなくしているのか」を推定していきます。


3. ゼロ過剰ポアソン分布で片方しか無いくつ下をモデリング


ゼロ過剰ポアソン分布を使ったモデリングでは、以下の2つの仮定をおきます。

  1. 靴下を片方なくすタイプの人と、なくさないタイプの人がいる(その割合(確率)はなくす確率thetaのベルヌイ分布に従う)
  2. 靴下を片方なくすタイプの人が紛失した靴下の数には個人差がある(平均=分散=λのポアソン分布に従う)

こういった仮定のもの、そういう眼差しのもとで、データを見ていきます。

イメージで書くと、次のようになります。 このとき、「ポアソン分布を仮定する」という仮定の中に、 「靴下なくすタイプだけど、今は片方だけを1足も持っていない」という場合も含まれています。 捨てちゃったとか、そういうかんじです。

ShoesZIPZIPExplain2.png

以下にベイズ統計の分析コードを紹介します。

所持数0のときは、ベルヌイ分布で0が出る場合と、ベルヌイ分布で1が出て、かつ、ポアソン分布で0が出る場合の尤度の積(対数尤度の和)を計算します。 所持数1以上のときは、ベルヌイ分布で1が出て、かつ、ポアソン分布で対象の値が出る場合の尤度の積(対数尤度の和)を計算します。

一気に結果の出力まで行きます。収束は事前に確認しております。

#stanに渡すデータ
datastan=list(Y=dat$Y, N=length(dat$Y))

#モデル
model<-"
data{
    int N;
    int Y[N];
}

parameters{
    real<lower=0,upper=1>theta;
    real<lower=0>lambda;
}

model{
    for (i in 1 :N){
        if (Y[i]==0){
            target += log_sum_exp(
                bernoulli_lpmf(0|theta),
                bernoulli_lpmf(1|theta)+poisson_lpmf(0|lambda)
            );
        } else {
            target += bernoulli_lpmf(1|theta) + poisson_lpmf(Y[i]|lambda);
        }
    }
}
"

#分析の実行
fitZIP<-stan(model_code = model,data=datastan)

#結果の出力
theta<-rstan::extract(fitZIP)$theta%>%as.numeric()
lambda<-rstan::extract(fitZIP)$lambda%>%as.numeric()
res<-data.frame(sample=c(theta,lambda),pars=rep(c("theta","lambda"),each=length(lambda)))

#表で結果の整理
summary(fitZIP)$summary[,c(1,4,6,8,9,10)]%>%data.frame()%>%
    formattable()

#グラフ化
Color=c("#5CACEE", "#00868B")
ggplot(res,aes(x=sample,y=..density..,colour=pars,fill=pars))+
    geom_histogram(position="identity",alpha=0.7)+
    geom_density(position="identity",alpha=0.4)+
    scale_color_manual(values=Color)+
    scale_fill_manual(values=Color)+
    facet_wrap(~pars,scales = "free")
ふむふむ。このモデルで考えると、
靴下をなくすタイプな確率は46.1%くらい。平均的な片方靴下の数は2.7足くらいですね。
Table1. ゼロ過剰ポアソン分布による推定結果

図にも出しておきましょう。事後分布はこちら。

ShoesZIPpars.png
Figure ゼロ過剰ポアソンモデルで推定した靴下なくす確率と、靴下数の平均(事後分布)

この結果をそのまま鵜呑みにすれば、

約半数の人たちが靴下をなくすタイプであるということになります。で、平均的に2足くらい持っているわけですね。


4. 打ち切りデータのゼロ過剰ポアソンモデリング


しかし、上の分析方法は実は適切とは考えにくい側面があります。

なぜなら、「7足以上」とお答えくださった方々が、全部7に吸収されてしまっているからです。

そのことを考慮していない分析は、結果が少しゆがんでしまっていると考えられます。

「7足以上の人が全部7に回答している。」

こういうデータを打ち切りデータと呼びますが、

以下のような方法で対処することができます。

詳細はやはり、『StanとRでベイズ統計モデリング』の第7章を御覧くださいませ。

以下は.stanファイルのStanモデルです。

data{
    int N_obs;
    int Y[N_obs];
    int N_cen;
    int L;
}

parameters{
    real<lower=0,upper=1>theta;
    real<lower=0>lambda;
}

model{
    for (i in 1 :N_obs){
        if (Y[i]==0){
            target += log_sum_exp(
                bernoulli_lpmf(0|theta),
                bernoulli_lpmf(1|theta)+poisson_lpmf(0|lambda)
            );
        } else {
            target += bernoulli_lpmf(1|theta) + poisson_lpmf(Y[i]|lambda);
            target +=bernoulli_lpmf(1|theta) + N_cen * poisson_lcdf(L|lambda);
        }
    }
}

上記のStanモデルを実行するRコードはこちら。

一気に結果の出力まで行います。

#
# censored data 
#
N_cen<-dat%>%dplyr::filter(Y==7)%>%summarize(count=n())%>%as.integer()
N_obs<-dat%>%dplyr::filter(Y<7)%>%nrow(.)
Y<-dat%>%dplyr::filter(Y<7)
censoredData=list(N_obs=N_obs,N_cen=N_cen,L=6,Y=Y$Y)

#stanfit
fitZIPcens<-stan(file = "ShoesZIP.stan",data=censoredData)

#結果の値を表で出力
summary(fitZIPcens)$summary[,c(1,4,6,8,9,10)]%>%data.frame()%>%
    formattable()

#HMC
thetaC<-rstan::extract(fitZIPcens)$theta%>%as.numeric()
lambdaC<-rstan::extract(fitZIPcens)$lambda%>%as.numeric()
resC<-data.frame(sample=c(thetaC,lambdaC),pars=rep(c("theta","lambda"),each=length(lambdaC)))

Color=c("#104E8B", "#006400")
ggplot(resC,aes(x=sample,y=..density..,colour=pars,fill=pars))+
    geom_histogram(position="identity",alpha=0.7)+
    geom_density(position="identity",alpha=0.4)+
    scale_color_manual(values=Color)+
    scale_fill_manual(values=Color)+
    facet_wrap(~pars,scales = "free")
打ち切りデータを考慮したゼロ過剰ポアソンモデルで考えると、
靴下をなくすタイプな確率は81.9%くらい。平均的な片方靴下の数は1.2足くらいですね。
Table2. 打ち切りデータを考えた場合の、ゼロ過剰ポアソン分布による推定結果

図にも出しておきましょう。事後分布はこちら。

ShoesZIPresult2.png
Figure ゼロ過剰ポアソンモデルで推定した靴下なくす確率と靴下数の平均(事後分布)

「無くした靴下の数はポアソン分布に従う」という仮定が効いてますね。

というのは、 「靴下をなくすタイプの人は8割くらいなんだけど、今は所持数ゼロ。平均は1.2足」という、 単純にデータを見ているだけでは分からない情報が得られています。

片方なくしたからもう片方を捨ててしまった。 とか、そういう方々を上手く反映できているのかなぁ。

真実はわかりません。しかし、そういう眼差しで考えた時の数値を求めることはできました。


5. 結果のまとめと総括


統計学の得意技は、「要約」です。

「平均値」などという一つの値にまとめて、理解しやすくする機能をもっています。

これは便利なのですが、

ついつい “いろんな人がいるんだよ” “ひとによって、個性があるんぞよ”ということを忘れてしまいます。

平均値も大事ですが、それにとらわれ過ぎると、個性を見落とします。バランスが大切ですね。

また、

統計学が教えてくれるのは、真実ではありません。

「こういうふうに考えてデータを見たいんだ」と思った時の、その答えを教えてくれます。

「片方しかない靴下の話」であっても、本当に色々な見方ができます。 最後に、「いろんなモノの見方」で捉えた、結果を表にまとめます。

どれが真実かはわかりません。

しかし僕は、

「靴下をなくすタイプの人が一定数いて」、

「なくす数は個人によって違う」、

「なくすタイプであっても、捨ててしまって今は片方だけの靴下を一足も持っていない人もいる」

という見方が現状一番説得力があると考えておりますので、

下記の表の一番右、

「打ち切りデータを考慮したゼロ過剰ポアソン分布の推定結果」を採用したいと思います。

ZIPFullResult.png

いやーしかし、ホントどこいったんやろ僕の片方の靴下…

Enjoy!!

Written on June 16, 2017