異なるN人のグループ分けの総数が知りたい

はじめに

識別可能なN人を任意のグループに分けるとき、その分け方の総数を知りたい。

はじめにPythonを使って力業で数え上げて計算してみた。次に、ベル数を使った計算を行った。

例:4人のグループ分けの総数

a, b, c, dの4人がいたとしよう。分け方は以下の15通り(のはず)。

グループの分け方 メンバーの分け方
oooo abcd
ooo | o abc | d
abd | c
acd | b
bcd | a
oo | oo ab | cd
ac | bd
ad | bc
oo|o|o ab|c|d
ac|b|d
ad|b|c
bc|a|d
bd|a|c
cd|a|b
o|o|o|o a|b|c|d

Pythonで計算してみる

n=4の場合

以下の手順でグループ分けの総数を求めた。

  1. 人数分のグループ(空リストn個からなるリスト)を用意する。
  2. 一人ずついずれかのグループに配属させる(リスト内リストにappendする)。
  3. 全員がいずれかのグループに配属したら、グル-プ分けの状態を整頓する(sorted関数でリスト内要素を並び替える)。
  4. 得られたグループ分けを記録する(group_listに追加する)。
  5. これ(1~4)をすべての配属パターンについて行う。
  6. group_listから重複を取り除く(一度set関数で集合にする)。
  7. group_listの要素数を数える(len関数)。

簡単に計算するために標準ライブラリitertoolsに含まれる関数を使う。
docs.python.org

import string
from itertools import product


def main(n):
    assert n <= len(string.ascii_letters), ValueError
    menbers = [s for s in string.ascii_letters[0:n]]
    group_list = []
    assignment_list = product(range(n), repeat=n)
    for assignment in assignment_list:
        groups = [[] for _ in range(n)]
        for e, i in enumerate(assignment):
            groups[i].append(menbers[e])
        groups = '|'.join(sorted([''.join(agroup) for agroup in groups]))
        group_list.append(groups)
    group_list = sorted(list(set(group_list)))
    print(len(group_list))
    print(group_list)
    return


if __name__ == '__main__':
    main(n=4)

実行すると次のようになる。

15
['a|b|c|d', '|ab|c|d', '|ac|b|d', '|ad|b|c', '|a|bc|d', '|a|bd|c', '|a|b|cd', '||abc|d', '||abd|c', '||ab|cd', '||acd|b', '||ac|bd', 
'||ad|bc', '||a|bcd', '|||abcd']

全部で15通りなので、手計算の結果と一致する。

n=8までの場合

N人の場合を計算してみる。値が大きくなると計算時間が非常に大きくなるので、n=8まででテスト。

import string
from itertools import product
import pandas as pd


def calc(n):
    assert n <= len(string.ascii_letters), ValueError
    menbers = [s for s in string.ascii_letters[0:n]]
    group_list = []
    assignment_list = product(range(n), repeat=n)
    for assignment in assignment_list:
        groups = [[] for _ in range(n)]
        for e, i in enumerate(assignment):
            groups[i].append(menbers[e])
        groups = '|'.join(sorted([''.join(agroup) for agroup in groups]))
        group_list.append(groups)
    group_list = sorted(list(set(group_list)))
    print(len(group_list))
    print(group_list)
    return len(group_list)


def main(n):
    df = []
    for k in range(n):
        df.append([k, calc(k)])
    df = pd.DataFrame(df, columns=['n_persons', 'n_groups'])
    df.to_csv('result.csv', index=False)
    return


if __name__ == '__main__':
    main(n=9)

力業でやっているのでまあまあな時間がかかる。出力されるCSVファイルを結果として示す。

n_persons n_groups
0 1
1 1
2 2
3 5
4 15
5 52
6 203
7 877
8 4140

人数が増えるとグループ分けの総数は爆発的に増える。

f:id:cyanatlas:20201021113213j:plain

計算方法の難点

上で示したPythonコードでは、 n^n 通りのグループ分けを計算した後で、重複を取り除くという方法をとっている。

そのため、最初に行うグループ分けの数が、 n = 9 のときでさえ3.8億通りであり、n が大きくなると計算負荷は爆発的に増加してしまう。

手元のノートパソコンでは、 n \geq 9 のときは応答が返ってこなくなってしまったので、計算できなかった。使用するメモリ量を減らすような書き方をしないといけないかもしれない。

ベル数を使って計算する

異なるn個のものをグループに分けるやり方の総数は、ベル数と呼ばれている。

実際、wikiに記載されているベル数を見ると

1, 1, 2, 5, 15, 52, 203, 877, 4140, ...

となっていて、先にPythonで数え上げたときの結果と一致している。

ja.wikipedia.org


グループ分けの一覧ではなく、グループ分けの総数だけが知りたいのであれば、数え上げまでせずにベル数を計算してやれば十分だ。

from scipy.special import comb
import pandas as pd
import matplotlib.pyplot as plt


def bell(n):
    if n == 0:
        return 1
    elif n == 1:
        return 1
    else:
        res = 0
        for k in range(n):
            res += comb(n - 1, k, exact=True) * bell(k)
        return res


def main(n):
    df = []
    for i in range(n + 1):
        df.append([i, bell(i)])
    df = pd.DataFrame(df, columns=['n', 'bell'])
    df.to_csv('bell.csv')

    fig, ax = plt.subplots()
    ax.plot(df['n'], df['bell'], marker='o')
    ax.set_xlabel('n')
    ax.set_ylabel('Bell number')
    fig.tight_layout()
    fig.savefig('bell.jpg', dpi=300)
    return


if __name__ == '__main__':
    main(n=10)

f:id:cyanatlas:20201021122511j:plain

数え上げる方法よりも早く計算できたため、より大きなnに対しても計算が可能になった。

[統計] 統計でおすすめの本(主に回帰)

はじめに

統計学を勉強するために読んだ本の中から、オススメの本をいくつか紹介しつつまとめてみる。

統計学と一口に言っても、トピックは多岐にわたるので、この記事では回帰に関する記述が豊富な書籍に絞った。

諸注意

  • 本稿の筆者は統計学の専門家ではありません。
  • 基本的には独立なデータに対する回帰の話題に絞っています。例えば時系列解析の書籍などは本稿では扱っていません。

対象者について

本の良し悪しは絶対的なものでなく、往々にして目的依存であると考えている。そこで、本を紹介する前に、私がどのような目的で統計学の本を読んでいるか説明しておく。

  • 統計手法がどのような目的を持つか知りたい(どんなタイプのデータから何を導けるか)
  • 統計手法がどのような数学的理論に基づいているのか知りたい(「このコマンドを打てばP値が出ます」では不満。とは言え、数理統計学の本を読んで定理と証明を厳密に知りたいというほどではない)
  • 得られた解析結果が統計学的にどう解釈されるか知りたい(有意差がある・ない、AIC最小ってつまりどういう"意味"?)

私と似たような目的を持って統計学の本を探している人には、以下で紹介する本が役に立つかもしれない。

統計の基本~単回帰

統計学入門

統計学入門 (基礎統計学Ⅰ)

統計学入門 (基礎統計学Ⅰ)

  • 発売日: 1991/07/09
  • メディア: 単行本

オススメ度 ★★★☆☆

良いところ

  • 検定や区間推定の細かいところまで詳しく書かれている。例えば、t統計量を用いた区間推定をどのように行うのか、その計算の詳細まで書かれており、付録の表を使えば手計算でも区間推定が可能なほど詳しい。
  • 統計学での考え方をはっきり書いている。例えば、区間推定で得られた区間がどういう意味なのか、詳しく説明されている。
  • 個人的には、第1章で紹介されている統計学の歴史がとても興味深かった。統計学が一見するとさまざまな手法の寄せ集めにも見える理由に納得がいった。

良くないところ

  • 書名に「入門」とあるものの、書き方が固く読みづらい。もちろん、固い本に慣れている人ならあまり支障はないかもしれない。
  • 話の行く末が見えづらい。本の前半は用語の紹介や古典的手法の説明が続くので(非常に大事なことが書かれているものの)、退屈な印象を受ける。
  • 基本的かつ古典的な手法しか載っていない。最終章まで読んでも、一変数の直線回帰(単回帰)までしかできるようにならない。もちろん、新しい手法を理解するために古典的な手法を抑えておくべきではある。

総評
統計学の考え方をちゃんと学びたい人には、自信を持ってオススメできる。この本が読める人には1冊目としてオススメしたいが、固すぎて読めないという人も多いと思う。

一般化線形モデル(GLM)

一般化線形モデルは正規線形モデル(直線回帰)の拡張。正規線形モデルを勉強したいという場合でも、その拡張である一般化線形モデルを勉強しておき、正規線形モデルをその特殊な例とみなす方が、見通しが良いように感じる。また、分散分析や共分散分析は、正規線形モデルに基づく分析手法であるため、これらの手法を学ぶ場合にも一般化線形モデルに触れておいて損はないと思う。本によっては分散分析も扱っている。

データ解析のための統計モデリング入門

オススメ度 ★★★★☆

統計を扱ったインターネットの記事ではよく見かける本で、緑本と呼ばれて親しまれているようだ。

良いところ

  • 前半でGLMとGLMM、後半でベイズ階層モデルが扱われており、統計モデルを概観することができる。とくに前半のGLMの話は非常に分かりやすい。
  • 実際に動くRのコードも載っており、自分で手を動かしながら学ぶことができる。
  • AICについてのまとまった記述がある。AICについては、非専門家向けに日本語で書かれた解説は少ないように思うので、これは貴重だと思う。

良くないところ

  • 数学的な話・理論的な話は省略されている。その辺りは無理にこの本で理解するよりも、他の本を読んだ方が良い。例えば下で紹介するDobson『一般化線形モデル入門』など。
  • 統計モデルを使った回帰の話に特化している。例えば相関に関する話題はない。まあこの点は良くないところというわけではないが・・・

総評
とても読みやすい本なので、GLMを勉強する際の1冊目としてオススメできる。

一般化線形モデル入門

一般化線形モデル入門 原著第2版

一般化線形モデル入門 原著第2版

オススメ度 ★★★★★

数式を使いながら一般化線形モデルを解説している本。

良いところ

  • 一般化線形モデルを数式を使ってちゃんと定式化し、最尤推定なども数式で示している。式展開は丁寧で分かりやすい。本文の説明が丁寧なので、途中の式を読み飛ばしながらでも十分内容を理解できる。
  • 指数型分布族、デザイン行列(計画行列)、連結関数(リンク関数)、スコア統計量、情報量といった用語は三章までで分かりやすく説明されている。
  • 最終章には、相関のあるデータに対する回帰手法が紹介されている。

良くないところ

  • 他の本に比べると数式が多いので(とは言え数理統計の本ほどではないが)、とっつきにくさは感じる。ただ、言葉の説明が曖昧で分かりにくいところは、数式を見に行けば厳密な意味が理解できるという利点はある。

総評
数式が多くてとっつきにくいかもしれないが、とても良い本なので是非読んで欲しい。とは言え、いきなり読むのは難しいので、RでGLMを使えるくらいにはGLMを理解してから読むことをオススメしたい。

自然科学の統計学

自然科学の統計学 (基礎統計学)

自然科学の統計学 (基礎統計学)

  • 発売日: 1992/08/01
  • メディア: 単行本

オススメ度 ★☆☆☆☆

上で紹介した『統計学入門』の続編。本全体としては必ずしも統計モデルに焦点を当てているわけではないが、第2章で正規線形モデル、第3章で分散分析、第4章で一般化線形モデルが扱われている。

良いところ

  • 第4章ではカウントデータに直線回帰をしてはいけない理由が丁寧に説明されていたり、最尤法がなぜ良いのか(最尤推定量は推定量としてどういった点で優秀なのか)が説明されていたりと、他の本にはあまり書かれていないような細かい話まで抑えている。

良くないところ

  • 前の『統計学入門』よりも数学的な話が多くなり、さらに固く読みづらい(実のところ第2章は難しくてまだ理解できてない・・・)。

総評
定量の話などは一読の価値があると思うが、万人にオススメできる本ではないと感じる。数理統計学の話題を数理統計学の本を勉強することなしに理解したいという人には良いかもしれない。

R

プログラミングが得意な人はRをわざわざ勉強する必要はないかもしれない(ネット上にも情報は豊富にある)。しかし実際には、R特有の落とし穴や頻出バグがあるので、一度体系的に勉強してみても良いと思う。自分の場合は、「Rでの正しいコーディングルールが知りたかった(スペースの入れ方や慣習的な命名規則)」というのと「グローバル変数の慣習的な扱いが知りたかった(ネットで探すとグローバル変数だらけのコードが多くてちょっと気持ち悪かった)」というのがモチベーションとしてあった。

アート・オブ・Rプログラミング

オススメ度 ★★★☆☆

アート・オブ・Rプログラミング

アート・オブ・Rプログラミング

いくつかRの本を開いてみたが、この本が一番良さそうだと感じた。信頼と実績のオライリー

ただし、これはRの文法書であって統計の本ではないことに注意。統計処理を行うRのコード集が欲しいという人には向かない。

[統計] GLMの最尤推定を自分でやってみる

はじめに

GLMで対数尤度を自分で計算して自分で最尤推定をやってみる、という回。

ロジスティック回帰

今回はロジスティック回帰を行ってみる。

ガス濃度とカブトムシの死亡数のデータを使う(Dobson, 2008, p.144)。

x n y
1.6907 59 6
1.7242 60 13
1.7552 62 18
1.7842 56 28
1.8113 63 52
1.8369 59 53
1.8610 62 61
1.8839 60 60

統計モデルは次のような式で表される。

 y_i \sim Binomial(n_i, \pi_i)

 \pi_i \sim N(\beta_1 + \beta_2 x_i, \sigma^2)

対数尤度

2パラメータのロジスティック回帰の場合の対数尤度関数を示す (Dobson, 2008, p.145)。

 l = \sum_{i=1}^N \left[y_i (\beta_1 + \beta_2 x_i) - n_i \log [1 + exp(\beta_1 + \beta_2 x_i) ] + \log\begin{pmatrix} n_i \\ y_i \end{pmatrix} \right]

対数尤度の計算

パラメータの値を変えて対数尤度を計算してみる。

from itertools import product

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
from scipy.special import comb


def loglikelihood_function(Y, N, X, beta1, beta2):
    res = 0
    for y, n, x in zip(Y, N, X):
        res += y * (beta1 + beta2 * x) - n * np.log(1 + np.exp(beta1 + beta2 * x)) + np.log(comb(n, y))
    return res


def calc_loglikelihood():
    df = pd.read_csv('beetle.csv')
    Y = df['y'].values
    X = np.ones(shape=[len(df), 2]) * np.nan
    X = df['x'].values
    N = df['n'].values
    RESOLUTION = 100
    beta1_arr = np.linspace(-80, -40, RESOLUTION)
    beta2_arr = np.linspace(15, 55, RESOLUTION)
    res = []
    for beta1, beta2 in tqdm(list(product(beta1_arr, beta2_arr))):
        loglikelihood = loglikelihood_function(Y=Y, N=N, X=X, beta1=beta1, beta2=beta2)
        res.append([beta1, beta2, loglikelihood])
    res = pd.DataFrame(res, columns=['beta1', 'beta2', 'loglikelihood'])
    res.to_csv('loglikelihood.csv')
    print(res[res.index == res['loglikelihood'].idxmax()])
    return


def plot_heatmap():
    df = pd.read_csv('loglikelihood.csv')
    df = df.pivot(index='beta2', columns='beta1', values='loglikelihood')
    df.to_csv('heatmap.csv')
    fig, ax = plt.subplots()
    sns.heatmap(data=df, ax=ax, cmap='YlGnBu')
    ax.set_xticklabels(['{:.2f}'.format(float(x.get_text())) for x in ax.get_xticklabels()])
    ax.set_yticklabels(['{:.2f}'.format(float(x.get_text())) for x in ax.get_yticklabels()])
    fig.tight_layout()
    fig.savefig('heatmap_loglikelohood.jpg', dpi=300)
    return


if __name__ == '__main__':
    calc_loglikelihood()
    plot_heatmap()
         beta1      beta2  loglikelihood
4947 -60.20202  33.989899      -18.72801

f:id:cyanatlas:20201018141256j:plain

図では分かりにくいが、対数尤度が最大になるのは \beta_1=-60, \beta_2=34 あたり。

statsmodelsでGLM

statsmodelsライブラリでロジスティック回帰を行ってみる。

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf


def glm_statsmodels():
    df = pd.read_csv('beetle.csv')
    df['z'] = df['n'] - df['y']
    glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
    fit = glm.fit()
    print(fit.summary())
    return


if __name__ == '__main__':
    glm_statsmodels()
                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:             ['y', 'z']   No. Observations:                    8
Model:                            GLM   Df Residuals:                        6
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -18.715
Date:                Sun, 18 Oct 2020   Deviance:                       11.232
Time:                        13:56:22   Pearson chi2:                     10.0
No. Iterations:                     6
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    -60.7175      5.181    -11.720      0.000     -70.871     -50.563
x             34.2703      2.912     11.768      0.000      28.563      39.978
==============================================================================

となって、対数尤度は \beta_1=-60, \beta_2=34 あたりで最大となることが分かる。この結果は、先に計算した結果と一致する。

実際の最尤推定の手順

ここでは、パラメータに対して総当たり的に対数尤度を計算し、対数尤度が最大になる場所を探して最尤推定量を求めた。しかし、この方法は非効率的である。

そこで、山登り的に計算を進めて最尤推定量を求める方法が考えられる。実際には、ニュートン法ニュートンラプソン法)をGLM向けに改良したスコア法が用いられ、対数尤度の一階微分(スコア統計量)が0になる場所が探索される。

メモ

このデータと最尤推定の例はDobson (2008)で紹介されているが、p.146の表7.3に示されている対数尤度は桁が1つ間違っていると思われる。具体的には、-186.235となっているところは-18.6235だろう。

[統計] GLMでの区間推定

はじめに

前回の記事では、正規分布の平均を区間推定する方法を説明した。

cyanatlas.hatenablog.com

同様に、一般化線形モデル(GLM)でもパラメータを区間推定することがあるが、この区間推定をどうやっているのか気になったので調べてみた。

GLM

パラメータを p 個持つGLMを考える。

 Y_i \sim \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p

正確に書けば次のようになる

 g(E[Y_i]) = \boldsymbol{x}^{T} \boldsymbol{\beta}
 \boldsymbol{x}^{T} = (x_{i1}, x_{i2}, \dots, x_{ip})
 \boldsymbol{\beta}^{T} = (\beta_1, \beta_2, \dots, \beta_p)

ここで g はリンク関数、\boldsymbol{x} は説明変数の行列、\boldsymbol{\beta} はパラメータのベクトルである。i 番目の応答変数  Y_i の期待値は、リンク関数を介して、説明変数とパラメータの線形結合によって説明される。

GLMのパラメータの最尤推定量が従う分布

GLMのパラメータ \boldsymbol{\beta}最尤推定\boldsymbol{b} は、近似的に正規分布に従うことが分かっている (Dobson, 2008, p.87)。

 \boldsymbol{b} \sim N(\boldsymbol{\beta}, \mathcal{I}^{-1})

ここで  \mathcal{I}^{-1} は、フィッシャー情報行列 \mathcal{I}逆行列である。フィッシャー情報行列は、対数尤度から計算される量であり、対数尤度関数が最尤推定量の周辺でどれだけ凸型になっているかを表している(詳しくはフィッシャー情報量 - Wikipedia)。

したがって、最尤推定量は、真の値を平均とし、フィッシャー情報行列の逆行列を分散共分散行列として持つ正規分布に従うということである。上の式を推定量と真の値の差で書き直せば、次のようになる。

 \boldsymbol{b} - \boldsymbol{\beta} \sim N(0, \mathcal{I}^{-1})

この式から、最尤推定量が真の値とどれだけ離れている可能性があるかを評価することができるので、信頼区間を求めることができる。

例:パラメータ2つのロジスティック回帰

例として次のデータのようなデータ(  N=40 )を考える。

x y
2.614672698355135 1
0.2311712788458009 0
4.130713626870605 1
3.4980589040126513 0
1.2242786762488267 1
0.5373907435482667 0
1.8970854531195753 1
0.3183052275680498 0
3.9466709572816496 0
3.5309785781654632 1
3.3761561081893414 0
0.720001443910317 0
2.064750554734893 0
2.009122195532624 0
2.2860103892990615 0
4.7927043238872375 1
1.1305070181779313 0
4.027287985243035 0
2.3250037833200765 0
4.554878548309989 0
0.6777507705415364 0
1.0177827317855659 0
4.042219354035978 1
1.5374167164460428 1
0.419902823072304 0
3.2341776629722996 0
1.2807428825075295 0
1.1185882637563382 0
0.37086337116461865 1
0.6935115140150583 1
4.775067677568351 1
3.5673200749785567 1
3.4595515158432804 1
3.195628563216748 0
0.3674178782304305 1
0.4163850367959582 0
0.23154185542532957 0
0.11134517325856519 0
4.287761228825091 1
2.2963424298511153 0

f:id:cyanatlas:20201018014744j:plain

パラメータを2つ持つロジスティック回帰を考える。

 logit(q_i) = \beta_0 + \beta_1 x_i

 i 番目のサンプル  Y_i が1となる確率  q_i は、切片  \beta_0 とパラメータ  x_i の係数  \beta_1 によって決まるというモデルである。

ちなみにデータを生成するときに使ったパラメータ(真のパラメータ)は \beta_0 = -2.0, \beta_1 = 0.8 である。

ではPythonのstatsmodelsライブラリを使って、ロジスティック回帰を行ってみよう。0の数もカウントする必要があるので、0の数をカウントした列を追加して、ロジスティック回帰を行う。

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv('test.csv')
df['z'] = 1 - df['y']
glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
fit = glm.fit()
print(fit.summary())

先に示したデータに対してロジスティック回帰を行い、print(fit.summary())で出力される結果を以下に示す。

                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:             ['y', 'z']   No. Observations:                   40
Model:                            GLM   Df Residuals:                       38
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -24.577
Date:                Sun, 18 Oct 2020   Deviance:                       49.153
Time:                        01:57:24   Pearson chi2:                     40.2
No. Iterations:                     4
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -1.4937      0.648     -2.306      0.021      -2.763      -0.224
x              0.4327      0.231      1.871      0.061      -0.021       0.886
==============================================================================

サマリーの右下に2つのパラメータの信頼区間が表示されている。今回はせっかくなので、この値をフィッシャー情報行列から求めてみよう。

statsmodelsでは、GLMオブジェクトのinformationメソッドを使うことで、フィッシャー情報行列を得ることができる(おそらく符号が逆になっていると思われたので、ソースコードの中では符号を反転させている)。

以下は  \beta_1区間推定を行う例。フィッシャー情報行列の逆行列を求め、その  (1, 1) 成分(  \beta_1 - b_1 の分散)から標準偏差を求める。±標準偏差×1.96の範囲に95%が含まれるので、ここから信頼区間を得ることができる。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv('GLM_test/test.csv')
df['z'] = 1 - df['y']
glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
fit = glm.fit()
x_est = fit.params['x']  # Maximum likelihood estimation of beta1
fim = - glm.information(params=fit.params)  # Fisher Information Matrix (FIM)
inv_fim = np.linalg.inv(fim)  # Inverse of FIM
std1 = np.sqrt(inv_fim[1, 1])  # Standard error of beta1
lower = x_est - 1.96 * std1
upper = x_est + 1.96 * std1
print('[{:.3f}, {:.3f}]'.format(lower, upper))

出力結果は

[-0.021, 0.886]

となって、先にサマリーで確認した信頼区間と一致する。

信頼区間の検証

GLMのパラメータの信頼区間の求め方が分かったところで、95%信頼区間が本当に95%信頼区間になっているか、データの生成と信頼区間の推定を繰り返して検証してみる。

先程と同じく  N=40 のデータにロジスティック回帰を行って \beta_1 の信頼区間を求める。これを M=1000 回繰り返して、信頼区間の的中率を計算する。

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
from tqdm import tqdm

# Parameters of a population
beta0 = -2.0
beta1 = 0.8

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 1000

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def logistic(x):
    return 1 / (1 + np.exp(-x))


def make_data(n):
    x_arr = np.random.uniform(low=0, high=5, size=n)
    err_arr = np.random.normal(loc=0, scale=0.3, size=n)
    probs = logistic(beta0 + beta1 * x_arr + err_arr)
    y_arr = np.array([np.random.binomial(n=1, p=p) for p in probs])
    df = pd.DataFrame({'x': x_arr, 'y': y_arr})
    return df


def calc_confint(df):
    df['z'] = 1 - df['y']
    glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
    fit = glm.fit()
    x_est = fit.params['x']  # Maximum likelihood estimation of beta1
    fim = - glm.information(params=fit.params)  # Fisher Information Matrix (FIM)
    inv_fim = np.linalg.inv(fim)  # Inverse of FIM
    std1 = np.sqrt(inv_fim[1, 1])  # Standard error of beta1
    lower = x_est - 1.96 * std1
    upper = x_est + 1.96 * std1
    print(fit.summary())
    return (lower, upper)


def estimation():
    df = []
    for i in tqdm(list(range(M))):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower <= beta1) and (beta1 <= upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(beta1, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を示す。

f:id:cyanatlas:20201018022108j:plain

この例では、的中率は96.5%であり、おおよそ95%になっている。

正規分布近似とWald統計量

以上の方法では、GLMのパラメータの最尤推定量は真の値を平均とする正規分布に従うという法則を用いていた。しかし、これは近似であるということに注意が必要である。データの数が少ないときには近似の精度が悪くなるので、正規分布に近似する代わりに、パラメトリックブートストラップを行う方が良いかもしれない。

最尤推定量と真の値の差が正規分布に従うのであれば、その二乗和はカイ二乗分布に従う(パラメータは独立でないので、分散共分散行列の逆行列を挟む必要がある。詳しくはDobson(2008, p.87)など)。この考え方で求まるのが、Wald統計量であり、このWald統計量を用いてWald検定が行われる。したがって、Wald検定は上で説明した信頼区間と類似した考え方だろうと思うのだが、詳細を知らないので確実なことは言えない。信頼区間の場合は分布の両側から2.5%を削って求めたが、Wald検定の場合はパラメータ  \beta が0から離れているかがポイントになるので、0に近い方から分布を5%削って有意性を調べるのだろうか。この辺りは久保(2012, p.51)に簡単な解説がある。とは言え、多くの場合はWald検定ではなく尤度比検定を使うだろうと思うので、この辺りの詳細はあまり気にする必要はないかもしれない。

参考文献

  • Annette J, Dobson (2008)『一般化線形モデル入門』(訳:田中豊・森川敏彦・山中竹春・冨田誠、原著第二版)共立出版
  • 久保拓弥 (2012)『データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC』(シリーズ 確率と情報の科学)岩波書店

[統計] 信頼区間とは何か

はじめに

与えられたデータの性質を知るために、回帰分析を行って係数(パラメータ)を推定するということが行われる。

このとき、得られた推定結果をp値や信頼区間AICといったさまざまな指標で評価する。

どの評価指標を使うべきかについては、これまでに数多くの議論がされてきた。ここでその詳細は追わないが、その議論の中で、p値よりも解釈の分かりやすい信頼区間を使うべきだとの指摘もあった。

では信頼区間とは何だったか。どうやって求めるのか。これが今回のテーマである。

点推定と区間推定

母集団のパラメータをどう推定するか。まず思いつくのは、最も可能性の高そうなある1つの値を示す方法である。これは点推定と呼ばれる。

点推定 point estimation:母集団のパラメータをある1つの値で指定する方法。

しかし、推定には常に誤差が伴うものである。点推定の結果が正しいとは限らないので、推定誤差の大きさも評価したい。

そこで、推定誤差の大きさも同時に評価できるような、区間推定という方法が考えられる。

区間推定 interval estimation:真のパラメータの値が入る確率がある値 1-\alpha 以上と保証される区間 [L, U] を示す方法。ここで 1-\alpha は信頼係数と呼ばれる。

区間推定に関する上記の説明は、直観的には理解しやすい。しかし、立ち止まってよくよく考えてみると、分からないところがある。

区間推定では「真のパラメータがある区間に入る確率」なるものを考えている。真のパラメータの値は、我々が調査する前からただ1つに決まっているはずなのに、である。

例えば、目の前に密閉されて中の見えない箱が1つあったとしよう。ある調査の結果、「その箱の中に猫が入っている確率は95%以上であると保証されます」と報告された。これは一体どう理解するべきだろうか。箱の中には最初から猫がいるかいないか決まっているはずであり、箱を開ける瞬間に猫がいるかどうか確率的に決まるとは考えにくい(少なくとも我々が通常認識している時空間的スケールでは)。では、同じ箱が100個あったら、95個には猫が入っているということだろうか?しかし、目の前には箱がたった1つあるだけである。

区間推定の定義

区間推定の正確な定義を見よう。『統計学入門』のp.225から引用する。

 同一の母集団から抽出した標本でも、標本ごとに信頼区間の推定値は変化する。\theta は未知ではあるが決まった定数である。したがって、一つの標本から信頼区間を具体的な数値として推定してやれば、これは信頼区間に含まれるか含まれないかのいずれかしかない。すなわち、具体的に数値として計算した現実の信頼区間に対して、”1-\alpha の確率で \theta を含む”ということはない。信頼区間の意味は、繰り返し多くの異なった標本について信頼区間をここで述べた方法によって何回も計算した場合、\theta区間内に含むものの割合が 1-\alpha となるというころである。

母集団から取り出された n 個の標本に基づいて、母集団のパラメータ \theta区間推定するということを、m 回繰り返すという状況を考えてみよう。

1回目の区間推定の結果を [L_1, U_1]i 回目の区間推定の結果を [L_i,U_i] と書くことにしよう。

すると m 回の区間推定のうち、区間内に真のパラメータ \theta が含まれる回数(\theta \in [L_i, U_i] である回数)は、m が十分大きいとき、おおよそ (1 - \alpha)m 回になるということである。

区間推定の妥当性

区間推定の正確な定義を参照した上で、区間推定の妥当性について考えてみる。

ここでは信頼係数 1-\alpha を0.95としておこう。

区間推定では「95%の確率で真のパラメータはだいたいこの区間にあるよ」ということが示される。

では「95%の確率で」というのがどういう意味かというと、区間推定の定義に立ち戻れば、「同じ母集団から今回と同じ数だけデータを取ってきて同じように区間推定した場合、95%の区間推定は真の値を含んでいるよ」ということである。

したがって、今回の区間推定が的中しているのか失敗しているのかは分からないけれども、区間推定が失敗する確率(5%)が無視できるとすれば、真のパラメータは今回の区間推定で得られた区間に入っているはずである。

以上が定義に基づく区間推定の解釈だが、正直言ってなかなかややこしい。

以下では具体的に区間推定を行う方法を説明する。

正規分布に従う母集団の平均を区間推定する(t 統計量を用いる方法)

正規分布  N(\mu, \sigma^2) に従う母集団の平均を、母集団から取り出した n 個の標本 (X_1, X_2, \dots, X_n) から推定する場合を考える。

標本平均  \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_it 分布に従うことが知られている。この結果を応用すると、母平均 \mu の信頼係数 1-\alpha の信頼区間

 \displaystyle [\bar{X} - t_{\alpha/2} (n - 1) s / \sqrt{n}, \bar{X} + t_{\alpha/2} (n - 1) s / \sqrt{n}]

となる。ここで s は標本分散の平方根であり、

 \displaystyle s^2 = \frac{1}{n-1} \sum_{i=1}^{n}(X_i-\bar{X})^2

から計算される。また、t_{\alpha} (n - 1) は自由度 n-1t 分布の上側確率 100\alpha %のパーセント点である。

では実際に計算してみよう。以下のPythonコードでは、平均 \mu=2.456 かつ標準偏差 \sigma=0.327正規分布に従う集団から、N=40 個の標本を取り出し、信頼係数 1-\alpha=0.95 で母平均の区間推定を行うという操作を M=100 回繰り返している。

import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt


# Parameters of a population
mu = 2.456
sigma = 0.327

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 100

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def make_data(n):
    return pd.DataFrame({'x': np.random.normal(loc=mu, scale=sigma, size=n)})


def calc_confint(df):
    n = len(df)
    x_bar = df['x'].mean()
    t_alpha_half = stats.t.ppf(q=1 - alpha / 2, df=n - 1)
    s2 = np.sum((df['x'].values - x_bar * np.ones(n)) ** 2) / (n - 1)
    s = np.sqrt(s2)
    lower = x_bar - t_alpha_half * s / np.sqrt(n)
    upper = x_bar + t_alpha_half * s / np.sqrt(n)
    return (lower, upper)


def estimation():
    df = []
    for i in range(M):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower < mu) and (mu < upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(mu, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を図にして示す。

f:id:cyanatlas:20201014000003j:plain

垂直な青線で各回の区間推定の結果を表示している。水平な赤線は真の母平均を示す。青線が赤線と交点を持つ場合、その回の区間推定は真のパラメータを含んでいた(的中した)ということになる。逆に、青線が赤線と交点を持たない場合は、その回の区間推定は真のパラメータを含んでいなかった(失敗した)ということになる。

この例では、100回の区間推定のうち、95回では区間推定で得られた区間に真の値が含まれていた。この結果は、信頼係数 1 - \alphaと対応している(乱数によってばらつきがあるので、必ず100回中95回というわけではない)。

ではもっと回数を増やしてみよう。M=1000 回の例を示す。

f:id:cyanatlas:20201014002510j:plain

この場合は、1000回中967回で区間推定が的中していた。

もっと増やしてみよう。最後は10000回。

f:id:cyanatlas:20201014005904j:plain

10000回中9537回で的中していた。おおよそ95%である。

正規分布に従う母集団の平均を区間推定する(ブートストラップ法)

母集団が正規分布に従う場合、標本平均は t 分布に従うことがあらかじめ分かっていた。上の例では、それを利用して、比較的簡単な計算で区間推定を行うことができた。

しかし、標本から得られた統計量の従う分布が分かっていない場合も多い。また、標本数が多い場合は漸近的に正規分布に従うが、標本数が少ないときは正規分布に近似すると精度が悪いという場合もある。

そこで、統計量の従う分布が数学的に分かっていなくても、コンピュータで乱数を生成することで注目する統計量の分布を生成してしまおう、というのがブートストラップ法である。

ここでは、先程と同じように、正規分布  N(\mu, \sigma^2) に従う母集団の平均を、母集団から取り出した n個の標本 (X_1, X_2, \dots, X_n) から推定する事例を、ブートストラップ法を使って考えたい(もちろん、本来であれば、この事例には上で説明したような t 分布を使う方法の方が良い)。

以下のPythonコードでは、平均 \mu=2.456 かつ標準偏差 \sigma=0.327正規分布に従う集団から、N=40 個の標本を取り出し、信頼係数 1-\alpha=0.95 で母平均の区間推定を行うという操作を M=100 回繰り返している。

ただし、区間推定の方法は先程とは異なる。得られた N=40 個の標本の平均と(標本)分散を持つ正規分布から N=40 個の乱数を生成してその平均を計算する、という操作を N_{bootstrap} = 2000 回繰り返すことによって、N=40 個の標本の平均が従う確率分布を計算し(この分布は t 分布と同じになる、はず)、この分布に基づいて区間推定を行っている。

実行には時間がかかるので注意。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm


# Parameters of a population
mu = 2.456
sigma = 0.327

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 100

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def make_data(n):
    return pd.DataFrame({'x': np.random.normal(loc=mu, scale=sigma, size=n)})


N_bootstrap = 2000


def calc_confint(df):
    samples = [np.mean(np.random.normal(loc=df['x'].mean(), scale=df['x'].std(), size=len(df))) for _ in range(N_bootstrap)]
    lower = np.quantile(samples, q=alpha / 2)
    upper = np.quantile(samples, q=1 - alpha / 2)
    return (lower, upper)


def estimation():
    df = []
    for i in tqdm(list(range(M))):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower < mu) and (mu < upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(mu, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を示す。おおよそ95%で的中している。

f:id:cyanatlas:20201014005200j:plain

M=1000 回の場合。こちらもおおよそ95%で的中している。

f:id:cyanatlas:20201014013759j:plain

M=10000 回で計算してみたところ、93.99%で的中していた。上の二つの例も踏まえると、95%よりも少ない方にバイアスがかかっているかもしれない。原因が分かったら追記したい。

まとめ

  • 区間推定で得られた区間の厳密な解釈には注意が必要(95%の”確率”で区間は真のパラメータを含んでいる)
  • 区間推定には、あらかじめ統計量の従う確率分布が分かっている場合(t 分布など)には、その分布を利用して区間推定が行われる。
  • 確率分布が分かっていない場合でも、ブートストラップ法によって区間推定を行うことができる。ただし、計算時間は長くなる。

参考資料・参考文献

[Python][Seaborn] Swarmplotが便利

Swarmplotの利点

Seabornライブラリのswarmplotがとても便利、というだけの記事。

import matplotlib.pyplot as plt
import seaborn as sns

df = sns.load_dataset('iris')
fig, ax = plt.subplots()
sns.swarmplot(x='species', y='sepal_length', data=df, ax=ax)
fig.savefig('sepal_length.jpg', dpi=300)
plt.close(fig)

f:id:cyanatlas:20201010174339j:plain

swarmplotを使うと、同じ値が複数あるときに横にずらして並べてくれるので、全体的な傾向が分かりやすい。

この図の場合、sepal_length(花のがくの長さ)の平均は、setosa種<versicolor種<virginica種、という傾向にありそうだが、同時に分散も大きくなっていそうだというのが見てとれる。

Swarmplotの限界

ただしswarmplotにも限界はあって、同じ値を持つデータ数が多くなりすぎると、横にずらすスペースがなくなって、最終的に重なってしまう。

下の図は、先程のデータセットを使ってpetal_length(花弁の長さ)をプロットしたもの。

f:id:cyanatlas:20201010175148j:plain

左のsetosa種では似たような値が多いために、横の方で複数の点が重なってしまっている。そうなると、頻度が分かりにくくなり、swarmplotを使う利点は小さくなる。

もしswarmplotで綺麗に描けないくらい値の似たデータが多いなら、カテゴリごとにヒストグラムを描いて並べるのが良いだろう。ただその場合は、ヒストグラムの切り方(bins)に注意を払う必要がある。

[Python][Pandas] いずれかの文字列と合致するorを含む行を抽出する

はじめに

文字列のリストが用意されていて、PandasのDataFrameからいずれかの文字列と合致するor含む行を取り出したいときがある。

Forループを使えば書けるのだが、Numpyを使えばもっとシンプルに書けることに気づいたのでまとめておく。

テスト用のDataFrameを読み込む

seabornライブラリからirisデータセットを読み込む。

import seaborn as sns

df = sns.load_dataset('iris')
df.to_csv('iris.csv')

irisデータセットは、文字列の種名を含んでいる。

sepal_length sepal_width petal_length petal_width species
0 5.1 3.5 1.4 0.2 setosa
1 4.9 3.0 1.4 0.2 setosa
2 4.7 3.2 1.3 0.2 setosa
3 4.6 3.1 1.5 0.2 setosa
4 5.0 3.6 1.4 0.2 setosa
5 5.4 3.9 1.7 0.4 setosa
6 4.6 3.4 1.4 0.3 setosa
7 5.0 3.4 1.5 0.2 setosa
8 4.4 2.9 1.4 0.2 setosa
9 4.9 3.1 1.5 0.1 setosa
10 5.4 3.7 1.5 0.2 setosa
11 4.8 3.4 1.6 0.2 setosa
12 4.8 3.0 1.4 0.1 setosa
13 4.3 3.0 1.1 0.1 setosa
14 5.8 4.0 1.2 0.2 setosa
15 5.7 4.4 1.5 0.4 setosa
16 5.4 3.9 1.3 0.4 setosa
17 5.1 3.5 1.4 0.3 setosa
18 5.7 3.8 1.7 0.3 setosa
19 5.1 3.8 1.5 0.3 setosa
20 5.4 3.4 1.7 0.2 setosa
21 5.1 3.7 1.5 0.4 setosa
22 4.6 3.6 1.0 0.2 setosa
23 5.1 3.3 1.7 0.5 setosa
24 4.8 3.4 1.9 0.2 setosa
25 5.0 3.0 1.6 0.2 setosa
26 5.0 3.4 1.6 0.4 setosa
27 5.2 3.5 1.5 0.2 setosa
28 5.2 3.4 1.4 0.2 setosa
29 4.7 3.2 1.6 0.2 setosa
30 4.8 3.1 1.6 0.2 setosa
31 5.4 3.4 1.5 0.4 setosa
32 5.2 4.1 1.5 0.1 setosa
33 5.5 4.2 1.4 0.2 setosa
34 4.9 3.1 1.5 0.2 setosa
35 5.0 3.2 1.2 0.2 setosa
36 5.5 3.5 1.3 0.2 setosa
37 4.9 3.6 1.4 0.1 setosa
38 4.4 3.0 1.3 0.2 setosa
39 5.1 3.4 1.5 0.2 setosa
40 5.0 3.5 1.3 0.3 setosa
41 4.5 2.3 1.3 0.3 setosa
42 4.4 3.2 1.3 0.2 setosa
43 5.0 3.5 1.6 0.6 setosa
44 5.1 3.8 1.9 0.4 setosa
45 4.8 3.0 1.4 0.3 setosa
46 5.1 3.8 1.6 0.2 setosa
47 4.6 3.2 1.4 0.2 setosa
48 5.3 3.7 1.5 0.2 setosa
49 5.0 3.3 1.4 0.2 setosa
50 7.0 3.2 4.7 1.4 versicolor
51 6.4 3.2 4.5 1.5 versicolor
52 6.9 3.1 4.9 1.5 versicolor
53 5.5 2.3 4.0 1.3 versicolor
54 6.5 2.8 4.6 1.5 versicolor
55 5.7 2.8 4.5 1.3 versicolor
56 6.3 3.3 4.7 1.6 versicolor
57 4.9 2.4 3.3 1.0 versicolor
58 6.6 2.9 4.6 1.3 versicolor
59 5.2 2.7 3.9 1.4 versicolor
60 5.0 2.0 3.5 1.0 versicolor
61 5.9 3.0 4.2 1.5 versicolor
62 6.0 2.2 4.0 1.0 versicolor
63 6.1 2.9 4.7 1.4 versicolor
64 5.6 2.9 3.6 1.3 versicolor
65 6.7 3.1 4.4 1.4 versicolor
66 5.6 3.0 4.5 1.5 versicolor
67 5.8 2.7 4.1 1.0 versicolor
68 6.2 2.2 4.5 1.5 versicolor
69 5.6 2.5 3.9 1.1 versicolor
70 5.9 3.2 4.8 1.8 versicolor
71 6.1 2.8 4.0 1.3 versicolor
72 6.3 2.5 4.9 1.5 versicolor
73 6.1 2.8 4.7 1.2 versicolor
74 6.4 2.9 4.3 1.3 versicolor
75 6.6 3.0 4.4 1.4 versicolor
76 6.8 2.8 4.8 1.4 versicolor
77 6.7 3.0 5.0 1.7 versicolor
78 6.0 2.9 4.5 1.5 versicolor
79 5.7 2.6 3.5 1.0 versicolor
80 5.5 2.4 3.8 1.1 versicolor
81 5.5 2.4 3.7 1.0 versicolor
82 5.8 2.7 3.9 1.2 versicolor
83 6.0 2.7 5.1 1.6 versicolor
84 5.4 3.0 4.5 1.5 versicolor
85 6.0 3.4 4.5 1.6 versicolor
86 6.7 3.1 4.7 1.5 versicolor
87 6.3 2.3 4.4 1.3 versicolor
88 5.6 3.0 4.1 1.3 versicolor
89 5.5 2.5 4.0 1.3 versicolor
90 5.5 2.6 4.4 1.2 versicolor
91 6.1 3.0 4.6 1.4 versicolor
92 5.8 2.6 4.0 1.2 versicolor
93 5.0 2.3 3.3 1.0 versicolor
94 5.6 2.7 4.2 1.3 versicolor
95 5.7 3.0 4.2 1.2 versicolor
96 5.7 2.9 4.2 1.3 versicolor
97 6.2 2.9 4.3 1.3 versicolor
98 5.1 2.5 3.0 1.1 versicolor
99 5.7 2.8 4.1 1.3 versicolor
100 6.3 3.3 6.0 2.5 virginica
101 5.8 2.7 5.1 1.9 virginica
102 7.1 3.0 5.9 2.1 virginica
103 6.3 2.9 5.6 1.8 virginica
104 6.5 3.0 5.8 2.2 virginica
105 7.6 3.0 6.6 2.1 virginica
106 4.9 2.5 4.5 1.7 virginica
107 7.3 2.9 6.3 1.8 virginica
108 6.7 2.5 5.8 1.8 virginica
109 7.2 3.6 6.1 2.5 virginica
110 6.5 3.2 5.1 2.0 virginica
111 6.4 2.7 5.3 1.9 virginica
112 6.8 3.0 5.5 2.1 virginica
113 5.7 2.5 5.0 2.0 virginica
114 5.8 2.8 5.1 2.4 virginica
115 6.4 3.2 5.3 2.3 virginica
116 6.5 3.0 5.5 1.8 virginica
117 7.7 3.8 6.7 2.2 virginica
118 7.7 2.6 6.9 2.3 virginica
119 6.0 2.2 5.0 1.5 virginica
120 6.9 3.2 5.7 2.3 virginica
121 5.6 2.8 4.9 2.0 virginica
122 7.7 2.8 6.7 2.0 virginica
123 6.3 2.7 4.9 1.8 virginica
124 6.7 3.3 5.7 2.1 virginica
125 7.2 3.2 6.0 1.8 virginica
126 6.2 2.8 4.8 1.8 virginica
127 6.1 3.0 4.9 1.8 virginica
128 6.4 2.8 5.6 2.1 virginica
129 7.2 3.0 5.8 1.6 virginica
130 7.4 2.8 6.1 1.9 virginica
131 7.9 3.8 6.4 2.0 virginica
132 6.4 2.8 5.6 2.2 virginica
133 6.3 2.8 5.1 1.5 virginica
134 6.1 2.6 5.6 1.4 virginica
135 7.7 3.0 6.1 2.3 virginica
136 6.3 3.4 5.6 2.4 virginica
137 6.4 3.1 5.5 1.8 virginica
138 6.0 3.0 4.8 1.8 virginica
139 6.9 3.1 5.4 2.1 virginica
140 6.7 3.1 5.6 2.4 virginica
141 6.9 3.1 5.1 2.3 virginica
142 5.8 2.7 5.1 1.9 virginica
143 6.8 3.2 5.9 2.3 virginica
144 6.7 3.3 5.7 2.5 virginica
145 6.7 3.0 5.2 2.3 virginica
146 6.3 2.5 5.0 1.9 virginica
147 6.5 3.0 5.2 2.0 virginica
148 6.2 3.4 5.4 2.3 virginica
149 5.9 3.0 5.1 1.8 virginica

DataFrameからある一つの文字列と合致する行を取り出す

DataFrameからある一つの文字列と合致する行を取り出すには、次のように書けば良い。

import seaborn as sns

df = sns.load_dataset('iris')
df_setosa = df[df['species'] == 'setosa'].copy()
df_setosa.to_csv('setosa.csv')

これで種名が'setosa'の行だけを取り出すことができる。

DataFrameから複数の文字列のいずれかと合致する行を取り出す

irisデータセットから、種名が'setosa'あるいは'virginica'である行だけを取り出したい。
この場合、Numpyのanyメソッドを使うとシンプルに書くことができる。

import numpy as np
import seaborn as sns

df = sns.load_dataset('iris')
species_list = ['setosa', 'virginica']
df_se_vi = df[
    np.array([df['species'] == aspecies for aspecies in species_list]).any(axis=0)
].copy()
df_se_vi.to_csv('setosa_virginica.csv')

この書き方の場合、species_listの要素数が増減しても変更の必要がないというメリットがある。

DetaFrameから複数の文字列のいずれかを含む行を取り出す

条件式の==str.containsメソッドに変えれば良い。

import numpy as np
import seaborn as sns

df = sns.load_dataset('iris')
part_list = ['set', 'vir']
df_set_vir = df[
    np.array([df['species'].str.contains(part) for part in part_list]).any(axis=0)
]
df_set_vir.to_csv('set_vir.csv')