Easy to type

個人的な勉強の記録です。データ分析、可視化などをメイントピックとしています。

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

その3です。今回は第2章の話をしていきます。

2章 MCMC正規分布の推測

内容

  • 解析的に事後分布を推定するのは難しい!
    • ベイズの定理を利用しても、分母に位置する正則化定数の積分が複雑で解析的に解けない
  • 事後分布に従う乱数を発生させる、という発想転換を行う(マルコフ連鎖モンテカルロ法; MCMC)

    • 乱数群の結果を比較して、解析する
    • HM法とGS法が具体的なアルゴリズムとしては代表的な存在
      • GS法は特定の事前分布であることが必要(やや恣意的)
      • HM法は自由な事前分布を使うことが出来る
    • 事前分布の設定 : 常識的な知識から決める
    • 乱数の発生 : 初期乱数は使用しない(burn in、あるいはwarm upと呼ばれる)
      • 値が登ったり下りたり(ドリフトした)形状が観測される場合には、好ましくない乱数生成となる
      • 複数chainを発生させて、全部で統一的な傾向が見られてない場合、好ましくない乱数生成が行われていると考えられる
  • 事後分布の解釈をどう行うか?

    • 点推定量
      • 事後期待値(EAP) = 平均値として扱うことが出来る
      • 事後中央値(MED) = 中央値として(ry
      • 事後確率最大値(MAP) = 最頻値(ry
      • 事後分散
      • 事後標準偏差
        • 与えられたデータを元に推定した分布自体の広がりがどれぐらいあるかを示したもの
      • 事後期待値の標準偏差(標準誤差; se)
        • EAPがどれぐらいの広がりを持つか示したもの
    • 確信区間(信用区間)
      • 母数が存在する確率がx % である区間を確信区間とする
    • 予測分布
      • 新たなデータが得られうる分布のこと
      • 事後予測分布
        • f(x^*|x)
        • 推定された母数の事後分布を用いて、未知データx^{*}を予測する統計モデルの平均を計算
        • x^*(t) \sim f(θ(t))
        • 事後分布から母数をサンプリングし、それを統計モデルにあてはめて予測をする
      • 条件付き予測分布
    • ベイズ的推測
      • 何らかの仮説を定めて、データから推定すること
      • 統計的仮説検定では帰無仮説に基づいたが、そんなことはしなくてよいのがベイズ統計
    • 生成量
      • 推定した母数を利用して、各種統計量の分布も計算することが出来る
        • 母数より生成される値なので、生成量という
    • 研究仮説が正しい確率
      • 真偽を判定する二値変数を利用して、研究仮説が正しい確率を評価することが出来る
      • 0 or 1しか使わないので、分散や信頼区間は使えない

条件付き予測分布の取扱について(2017-05-28追記)

ここまでで書いた概要部分でも説明しましたが、議論のされている部分なので補足説明です。 このテキストではデータxが得られた状態で、次に得られるであろう未知データ x^{\ast}についての予測方法を2つ紹介しています。

1つは事後予測分布です。事後予測分布の数式は p_{pred}(x^\ast|x) = \int p(x^\ast|\theta)p(\theta|x)d\thetaと記述されます。これを踏まえたやり方です。 数式だけ見ると混乱しがちなのですが、MCMCで計算する時にはそんなに難しくありません。 MCMCの各サンプリングにおいて得られる\thetap(\theta|x^\ast)を示していますし、この\thetaを構築したモデルに代入して得られる値が分布p_{pred}(x^\ast|x)からサンプリングされたうちの1つということになります。 MCMC積分の代わりをしている、という認識ですね。

もう1つは条件付き予測分布です。これはMCMCの結果として得られる \thetaEAPやMAPといった点推定値をモデルに代入して、予測分布を得るというものです。 これは理解もし易いし、簡単です。分布の中で最もそれっぽいと理解出来た値をモデルに代入して予測をすることができます。 著者の豊田先生としては、

  • 理解が簡単
  • 速度が速い(生成量や予測量として計算しなくて済む)
  • オンライン性(事後予測分布は過去データ xが必要だが、条件付き予測分布は点推定値だけでよい)

という利点を挙げています。

しかしながらリンク先にもあるように、条件付き予測分布には色々問題点もありますし、何より個人的には「もったいない」という点が目立つと考えています。 せっかくMCMCまでして分布を計算したんですから、それを有効に活用しましょう。 それに「仮説が正しい確率」を計算するには、どっちにしろ事後予想分布が必要になります。

以上の理由で、以後の説明では条件付き予測分布は使用せず、事後予想分布で推定を行っています。

章末問題

データの準備

  • 1章のデータをそのまま利用します
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()

Stanによる統計モデルの構築

  • MCMCをStanで行います
  • Stan言語で書かれたファイルを読み込んで、C++コンパイルすることで、高速にサンプリングを行うことが可能となります。
  • 標本データや問題文を見た時に、値の範囲について情報が全く得られなかったので、パラメータ分布には0以上という制限しか利用していません。
import os
import pystan
import pickle

# Stanのモデルを読み込んでコンパイルする
stan_file = os.path.join("stan", "g1_mean.stan")
stan_file_c = os.path.join("stan", "g1_mean.pkl")
model = pystan.StanModel(file=stan_file)
with open(stan_file_c, "wb") as f:
    pickle.dump(model, f)
# Stanのモデル
data {
    int<lower=0> n ;
    real<lower=0> x[n] ;
    int<lower=0> c ;
    real<lower=0, upper=1> percent ;
    real ces ;
}

parameters {
    real<lower=0> mu ;
    real<lower=0> sigma ;
}

model {
    x ~ normal(mu,sigma);
}

generated quantities{
    real xaste ;
    real sigma2 ;
    real cv ;
    real es;
    real percent_point ;
    real<lower=0, upper=1> prob_under_c_cdf ;
    real rate_c ;
    real<lower=0, upper=1> prob_mu_under_c ;
    real<lower=0, upper=1> prob_under_c_bool ;
    real<lower=0, upper=1> prob_es_under_ces ;
    real<lower=0, upper=1> prob_prob_under_c_bool_over_percent ;

    xaste = normal_rng(mu,sigma) ;
    sigma2 = sigma * sigma ;
    cv = sigma / mu ;
    es = (mu - c) / sigma ;
    percent_point = mu + normal_cdf(percent, 0, 1) * sigma ;
    rate_c = xaste / c ;
    prob_under_c_cdf = normal_cdf(c, mu, sigma) ;
    prob_under_c_bool = xaste < c ? 1 : 0 ;
    prob_mu_under_c = mu < c ? 1 : 0 ;
    prob_es_under_ces = es < ces ? 1 : 0 ;
    prob_prob_under_c_bool_over_percent = prob_under_c_bool > percent ? 1 : 0 ;
}

統計モデルによる事後分布のサンプリング

  • コンパイルしたモデルを使って、正規分布からデータが得られた際に、正規分布の母数がどのように分布するか調べます
  • 母数の分布から、各種効果量を推定しています。
  • 仮説が得られる確率については、各サンプリング値において判定を行い、その結果のbool値の頻度で確認します。
    • 付録でついてきたスクリプトにはここらへんの記述がありませんので、何か問題が有った場合には指摘頂けると助かります。
import pandas as pd
import pickle
import pystan
import matplotlib
import matplotlib.pyplot as plt
from IPython.core.display import display
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (10, 40)

# Stanで使うデータの用意
stan_data = {"n": x.size,
                 "x": x, 
                 "c": 41, 
                 "percent": 0.70, 
                 "ces": -0.05}
# 興味のあるパラメータの設定
par = ["mu",
       "sigma",
       "xaste",
       "sigma2",
       "cv",
       "es",
       "percent_point",
       "rate_c",
       "prob_under_c_cdf",
       "prob_under_c_bool",
       "prob_mu_under_c",
       "prob_es_under_ces",
       "prob_prob_under_c_bool_over_percent"]
# 興味のある母数の%点を書く
prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975]

# コンパイルしたモデルの読み込み
stan_file_c = os.path.join("stan", "g1_mean.pkl")
with open(stan_file_c, "rb") as f:
    model = pickle.load(f)

# MCMCでサンプリング
fit = model.sampling(data=stan_data,
                         pars=par,
                         iter=21000,
                         chains=5,
                         warmup=1000,
                         seed=1234,
                         algorithm="HMC")


# 事後分布の表を取得、可視化
summary = fit.summary(pars=par, probs=prob)
summary_df = pd.DataFrame(summary["summary"],
                                       index=summary["summary_rownames"],
                                       columns=summary["summary_colnames"])
display(summary_df)

# 事後分布の可視化
fit.plot(pars=par)
plt.show()
  • 結果(表と分布)
mean se_mean sd 2.5% 5% 25% 50% 75% 95% 97.5% n_eff Rhat
mu 40.660 0.010 0.649 39.383 39.603 40.224 40.659 41.091 41.735 41.947 4120.000 1.004
sigma 6.523 0.014 0.468 5.660 5.797 6.200 6.499 6.828 7.321 7.488 1154.000 1.008
xaste 40.693 0.021 6.560 27.784 29.927 36.314 40.679 45.108 51.424 53.495 96434.000 1.000
sigma2 42.772 0.182 6.178 32.033 33.604 38.445 42.237 46.623 53.602 56.063 1153.000 1.008
cv 0.160 0.000 0.012 0.139 0.142 0.152 0.160 0.168 0.181 0.185 1372.000 1.006
es -0.052 0.002 0.099 -0.248 -0.215 -0.120 -0.053 0.014 0.112 0.144 4065.000 1.004
percent_point 45.605 0.019 0.744 44.204 44.418 45.099 45.583 46.092 46.862 47.132 1509.000 1.009
rate_c 0.993 0.001 0.160 0.678 0.730 0.886 0.992 1.100 1.254 1.305 96434.000 1.000
prob_under_c_cdf 0.521 0.001 0.039 0.443 0.455 0.494 0.521 0.548 0.585 0.598 4055.000 1.004
prob_under_c_bool 0.520 0.002 0.500 0.000 0.000 0.000 1.000 1.000 1.000 1.000 97815.000 1.000
prob_mu_under_c 0.703 0.007 0.457 0.000 0.000 0.000 1.000 1.000 1.000 1.000 4861.000 1.001
prob_es_under_ces 0.510 0.007 0.500 0.000 0.000 0.000 1.000 1.000 1.000 1.000 5481.000 1.002
prob_prob_under_c_bool_over_percent 0.520 0.002 0.500 0.000 0.000 0.000 1.000 1.000 1.000 1.000 97815.000 1.000

f:id:ajhjhaf:20170523020636p:plain

Rhat値が全パラメータで1.1以下になっていることや、Traceplotの様子からMCMCが収束していることが確認できると思います。 というわけで推定値は一応問題がありません。実際に推定値を用いて推論を行うことが出来ます。

仮説と検証

平均値
  • このデータにおける平均値はいくつか
    • 事後分布のEAPが40.660(0.650)であるため、40.660と推定される
  • このデータにおいて、平均値が95%の確率で分布する区間はどこか
    • 事後分布の平均値の2.5%区間と97.5%信頼区間は39.383から41.947である
  • このデータにおいて平均値が片側で95%の信頼区間を持つ区間はどこか
    • 事後分布の平均値の95%信頼区間より、41.735以下になる確率が95%である
    • もしくは39.603以上になる確率が95%である
標準偏差
  • 標準偏差の平均値、両側/片側95%信頼区間はどこか
    • 6.523(0.468)、5.660から7.488、5.797以上あるいは7.321以下である
予測分布
  • 新しい人からデータをとったときに、95%の確率で何点を取るでしょうか
    • 27.769から53.520点の間のいずれかの点数をとる
生成量
  • 分散の点推定と95%信頼区間推定
    • 42.772(6.178)[32.033, 56.063]
  • 基準値を41とした際の、変動係数の点推定と95%信頼区間推定
    • 0.160(0.012)[0.139, 0.185]
  • 基準値を41とした際の、効果量の点推定
    • -0.524(0.099)
  • 基準値を41とした際の、効果量の95%区間推定(両側/片側)
    • 両側[-0.248, 0.144]、片側[-0.215, ] or [, 0.112]
  • 70%点の点推定と区間推定
    • 45.605(0.744)[44.204, 47.132]
  • 41以下の値が観測される確率と95%信頼区間
    • 分布関数を利用 : 0.521(0.039)[0.443, 0.598]
    • 確率的命題を利用 : 0520
  • 41と測定値の比の、点推定値と95%信頼区間推定値
    • 0.993
研究仮説が正しい確率
  • 平均値が41より小さい確率
    • 0.703
  • 基準値を41とした際の効果量が、-0.05より小さい確率
    • 0.510
  • 41未満の測定値が観測される確率が、70%より大きい確率
    • 0.512

はじめての 統計データ分析 ―ベイズ的〈ポスト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値時代〉の統計学―

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

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

目次

本の紹介

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

d.hatena.ne.jp

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

良い点

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

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

Research Questionベースの解析

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

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

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

答えが用意されている

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

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

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

悪い点

入門書ではない

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

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

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

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

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

なにをするのか

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

ブログ始めました。

自己紹介

ブログ始めました。

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

ブログの目的

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

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

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