はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その2

前置き

環境構築

まずPythonを実行する環境を作ります。 一応メモとして書きますが、他にもPythonの導入の記事はたくさんありますので、そちらも参考ください。 以下のものはOSXを想定しています。

自分は基本的にはpyenv + pyenv-virtualenvを利用しています。pyenvで切り替えつつ、学習する本の内容に合わせてpyenv-virtualenvで環境を作っています。 導入は各リポジトリのページに書いてあるとおりです。

github.com github.com

pythonは3系を使います。今回利用するライブラリはそんなに変なものではないので、2系である必要はありません。 解析はスクリプトでやるよりも、jupyterを使う方が僕は好きです。Stan言語をいじれない問題点はありますが、パラメータを変えるのが楽です。

pyenv install 3.6.1
pyenv virtualenv 3.6.1 3.6.1_postp
mkdir -p ~/Desktop/postp # 作業ディレクトリは適当に変えてください。
cd ~/Desktop postp
pyenv local 3.6.1_postp
pip install numpy scipy pandas matplotlib cython seaborn jupyter
pip install pystan

1章 データの整理とベイズの定理

内容

一応箇条書きで書いていますが、メモです。 ちゃんと勉強したいときは、本文を参考にしてくださいね。

• 分布の紹介
  • 正規分布
  • 一様分布
• データの解釈
  • なぜ正規分布にあてはめる?
    • 観測データの実際の分布は、より複雑な謎の分布に従うと予測される
    • 正規分布はパラメータが少なく、扱いが容易

       * データが分布にあてはめられることで、分布の数式を元にしたマニピュレーションが可能
      

  • どのようにあてはめる?
    • 頻度主義 : 標本平均、標本偏差を推定して、あてはめる
    • ベイズ統計 : データからパラメータの分布を推定し、その曖昧さを踏まえて分布を推定
      • ベイズの定理より、データが得られた時のパラメータ分布が得られる

              * ベイズの定理のために必要なもの
        

        • 尤度
        • 事前分布
          • 無情報事前分布の利用
        • 正則化定数(扱いが難しい!)
        • そこでMCMCですよ。

章末問題

• pythonにて解析をしていきます。

データの準備

僕の好みですが、printではなくloggingモジュールを利用して途中経過を書いていきます。計算時間なども把握できるのが便利です。

import numpy as np
from logging import getLogger, Formatter, StreamHandler, DEBUG

# printではなくloggerを使う
def get_logger():
    logger = getLogger(__name__)
    logger.setLevel(DEBUG)
    log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s'
    formatter = Formatter(log_fmt)
    stream_handler = StreamHandler()
    stream_handler.setLevel(DEBUG)
    stream_handler.setFormatter(formatter)
    logger.addHandler(stream_handler)
    return logger

# データの準備
x = np.array([36,38,51,40,41,52,43,31,35,37,49,43,43,41,36,53,43,26,45,37,33,38,33,35,36,28,46,41,32,49,43,38,46,46,46,45,44,40,38,37,35,39,31,55,48,32,37,37,45,39,42,40,40,50,38,51,29,44,41,42,43,36,38,33,32,42,43,40,46,54,37,24,47,35,35,47,38,31,41,39,40,43,37,45,38,42,48,43,38,48,47,44,42,36,50,36,55,51,38,33])
logger = get_logger()

度数分布表

いきなりですが、pandasでは階級幅を決めるメソッドが無いので超めんどくさい… 手動でデータを入れるbinsを作成した後、pd.cutを使って仕分けます。

import pandas as pd
df = pd.DataFrame(x, columns=["weight"])

width = 5
minim_range = int((df.min() / width).apply(np.floor) * width)
maxim_range = int((df.max() / width).apply(np.ceil) * width)
bins = range(minim_range, maxim_range + 1, width)
c = pd.cut(df["weight"], bins)

d_df = df.groupby(c).count()
d_df.columns = ["count"]
d_df.columns.name = "range"
d_df.index.name = None 
d_df

ヒストグラムの作成

seabornを使って可視化します。matplotlibではないのは、やはり好みです pandas.plotと違ってseaborn / matplotlibのhistメソッドはbinsに配列を受け取ってくれます。

import seaborn as sns
import pandas as pd
%matplotlib inline
df = pd.DataFrame(x, columns=["weight"])

width = 5 
minim_range = int((df.min() / width).apply(np.floor) * width)
maxim_range = int((df.max() / width).apply(np.ceil) * width)
bins = range(minim_range, maxim_range + 1, width)
sns.distplot(df["weight"], kde=False, bins=bins)

統計量

最頻値だけはメソッドが無いので、countをとって最大の値を持つものを返しています。 pandasのvar, stdメソッドはデフォルトパラメータだと母分散、母標準偏差を返す事に注意してください。 numpyは標本分散、標本標準偏差を返す。謎の挙動。

import pandas as pd
df = pd.DataFrame(x, columns=["weight"])

# 標本平均
df_mean = df["weight"].mean()
logger.info("平均値は{0}".format(df_mean))

# 標本分散
df_var = df["weight"].var(ddof=0)
logger.info("標本分散は{0}".format(df_var))

# 標本標準偏差
df_std = df["weight"].std(ddof=1)
logger.info("標本標準偏差は{0}".format(df_std))

# データのソート
df["weight"].sort_values()

# 最大値
df_max = df["weight"].max()
logger.info("最大値は{0}".format(df_max))

# 最小値
df_min = df["weight"].min()
logger.info("最小値は{0}".format(df_min))

# 中央値
df_medium = df["weight"].median()
logger.info("中央値は{0}".format(df_min))

# 最頻値
df_most_freq = df["weight"].value_counts().iloc[0]
logger.info("最頻値は{0}".format(df_most_freq))

正規分部の扱い

scipy.statsのnormオブジェクトを利用します。 片側だけから特定の確率になる値を求めるメソッドが無いですが、intervalメソッドで両側の場合の値を求めて、2で割って考えれば良いです。 これは正規分布の対称性を利用したものですね。

import pandas as pd
from scipy.stats import norm 

df = pd.DataFrame(x, columns=["weight"])

# パラメータを設定
df_mean = df["weight"].mean()
df_std = df["weight"].std(ddof=1)
logger.info("正規分布のパラメータ: 平均は{0}、標準偏差は{1}".format(df_mean, df_std))

# 分布の平均が40.64なので、40付近の値が観測されやすい

# 45以上の値が観測される確率
# 1から45までの値が観測される確率を引く(1-cdf)
# scipyでは、sfを使えば良い
from scipy.stats import norm 
prob_upper_45 = norm.sf([45], loc=df_mean, scale=df_std)[0]
logger.info("45以上の値が観測される確率は{0}".format(prob_upper_45))

# 35以上45未満の値が観測される確率
# 45以下が観測される確率から、35以下が観測される確率を引く
prob_upper_35 = norm.cdf([35], loc=df_mean, scale=df_std)[0]
prob_upper_45 = norm.cdf([45], loc=df_mean, scale=df_std)[0]
prob_between_35_45 = prob_upper_45 - prob_upper_35
logger.info("35以上45未満の値が観測される確率は{0}".format(prob_between_35_45))

# 95%信頼区間
# 正規分布において、平均 ± 1.96 * 標準偏差の範囲が該当することを利用する
#min_edge = df_mean - df_std * 1.96
#max_edge = df_mean + df_std * 1.96
#prob_min = norm.cdf([min_edge], loc=df_mean, scale=df_std)[0]
#prob_max = norm.cdf([max_edge], loc=df_mean, scale=df_std)[0]
# prob_sig_range = prob_max - prob_min
# が、scipyだとintervalを利用すれば良い
min_edge, max_edge = norm.interval(0.95, loc=df_mean, scale=df_std)
logger.info("95%信頼区間は{0}から{1}".format(min_edge, max_edge))

# p(x>a) = 0.05であるような点a
# F(平均 + 1.64 * 標準偏差)が大体0.95になる
#prob_less_095 = norm.sf([df_mean + df_std * 1.64], loc=df_mean, scale=df_std)[0]
# 何も情報が無いとしたら、正規分布性を元に、90%信頼区間を出して、上端を利用すれば良い
_, less_095 = norm.interval(0.9, loc=df_mean, scale=df_std)
prob_less_095 = norm.cdf([less_095], loc=df_mean, scale=df_std)[0]
logger.info("a = {0} であれば、p(x>a)となる確率は{1}である".format(less_095, prob_less_095))

# 3つの4分位点
# 50%信頼区間と中央値を出せば良い
quad_first, quad_third = norm.interval(0.5, loc=df_mean, scale=df_std)
quad_second = df_medium
logger.info("第1四分位点は{0}, 第2四分位点は{1}, 第3四分位点は{2}".format(quad_first, quad_second, quad_third))

Rで出来ることがPythonだとめんどくさかったりしますね。 特定の幅長での度数表や、最頻値について簡単に求めることが出来る方法を知っていらっしゃる方が居たら、是非コメントを頂けると助かります！

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その1

このブログの開設目的の、

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

についての記事インデックスです。

目次

• その1 ←このページ
• その2 環境構築と第1章 データ整理とベイズ概要
• その3 第2章 : MCMCによる正規分布のパラメータ推定
• その4 第3章 : 独立した2群の差の推定
• その4.1 第3章おまけ : 変分ベイズによる推定
• その5 第4章 : 対応ある2群の差の推定
• その6 第5章 : 要因毎の効果推定
• その7 第6章 : カウントデータと比率データの推定

本の紹介

この本は、昨今はびこる統計データ分析についてベイジアンの知見から解説した本です。 他の方の感想では、 statmodeling.hatenablog.com

d.hatena.ne.jp

のようなものがあります。個人的に一冊まるまるよんだ感想としては、次のような印象を持ちました。

良い点

MCMCについての細かい知識は要らない

ベイズの定理やMCMCを利用した多くの解析の本では、数式がガッツリ出てきて頭を抱える方も多いと思います。 しかし、この本ではStan言語を利用して解析を行っていきます。そのため細かい数式の変形や、推論の部分についてはStanがなんとかしてくれますので、本にも書いてありませんし、やる必要もありません。

Research Questionベースの解析

実際にデータ分析をする際には、様々なデータ提供者からの要望に答える形で解析をすることが多いと思います。 この本では、そのような問いをResearch Questionと呼び、各章で該当するResearch Questionを最初に挙げて、解決する形で進めていきます。 結果として、どのような形で今やっている解析が活きるものなのか、具体的にイメージしながら進める事ができます。

ベイズベースでのデータ解析への「入門」

豊田先生が後書きにも記していますが、普通データ解析や統計学を勉強するときには、頻度主義に基づく統計的仮説検定を勉強します(僕もそうしました)。 一方でベイズの定理を基づくベイジアン統計学は、最近計算機処理が高速化したことやデータ量が多くなっていることが原因で、頻繁に使われています。 しかしベイジアン統計学は、頻度主義とはぜんぜん違うため、入門が難しいと僕も思います。その中で本書は、統計的仮説検定の手法でやるような内容をベイジアン統計学で行った場合にどうなるかを解説しています。 その点で、最初に取り組む際に何が利点となるのか分かりやすいと思います。

答えが用意されている

最近は多いですが、解析に使っているスクリプトがR + Stanで用意されています。

朝倉書店｜はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

更に、各章末問題は巻末に答えがかなり丁寧に書いてあります。 大学の教科書もそうですが、答えがない場合自分が間違っていても気が付かないことがあります。答えがあるのは非常に助かりますね。

悪い点

入門書ではない

幾つもの指摘がありましたが、入門書ではありません。 たとえば、

• Stanの文法や使い方について丸々省略されている

• 統計モデリングの必要性が良くわからない

• ベイズの定理の導出なども駆け足気味

といった点が目立ちました。1つ目、2つ目については久保先生の通称緑本、ベイズの定理について入門レベルから丁寧に勉強したい場合には超入門が良かったと思います。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

図解・ベイズ統計「超」入門 あいまいなデータから未来を予測する技術 (サイエンス・アイ新書)

なにをするのか

以降からは、各章の簡単な解説をした上で、Pythonで章末問題を解いていきたいと思います。 Stan言語については、PythonのインターフェイスであるPyStanを利用します。 といっても、Stan言語で書かれている部分についてはRと共通して使うことが可能なので、PyStanを使うためにどのような小手先のテクニックがあるかを書いていくつもりです。

ブログ始めました。

自己紹介

ブログ始めました。

普段は次世代シーケンサー(NGS)を用いた遺伝統計学や、ゲノム配列から得られる知見の解析、それらのツール・DB開発などを行っています。

ブログの目的

自身で勉強したことを内部のWikiに書いていたのですが、もっと外部に公開してマサカリを投げてもらったり、何かしらの役に立つことを期待してブログを開設します。 基本的には、

• Pythonによるモデリング、機械学習

• Javascriptを用いたデータ可視化(D3.js)、アプリケーション開発

• 次世代シーケンサーから得られるデータの解析について

を主に書いていきたいと思います。特にPythonについては日本語知識があまり転がっていないので、Rや海外サイトから輸入することが多いかもしれません。 3番目のNGS解析については、需要がなさそうなのであまり書かない気がします笑

それでは宜しくお願いします。Amazon.co.jpアソシエイトです。