pythonのprastクラスによる集計(4) 分布のチェック

2018/10/04

 今回はデータの分布を確認する方法です。

 最小値、最大値、平均値、あるいは最頻値などを取得します。

 データが正規分布に即しているといえるかどうかの検証も行います。

 処理するcsvファイルは第1回〜第3回と同じ data01.csv です。

 当Webページで紹介するスクリプトや素材データ一式は、
prast04.zip という圧縮ファイルに同梱しておきます。

    


《このページの目次》


    

1. pandasライブラリのdescribeメソッド

 pythonのpandasライブラリには describe() という便利なメソッドがあります。

 身長などの数値データに適用した場合は
最小値、最大値、平均値などの分布の要約を示してくれます。

 数値データ以外(性別とか意見などの種類を扱いデータ)のときは
種類数や最頻値(最も頻度の多い種類)を示します。

 Prastクラスにも describe() を設けましたが、基本的に同じ情報を返します。

 そこで、まずは pandasの describe() について概要を記します。

    

(1) 数値データへの適用

 describe() は、DataFrame または Series に適用します。

 DataFrameに適用した場合は、その各列(Series)に適用した結果を
取りまとめて出力する形になります。

 data01.csv の height に適用する場合は次のようにします。

1# (coding: cp932)
2import pandas as pd
3dtf = pd.read_csv("data01.csv")
4ser = dtf["height"].describe()
5print(ser)

 出力結果は下のとおり。

count    392.000000
mean     163.461990
std        7.090396
min      144.800000
25%      158.675000
50%      163.600000
75%      168.025000
max      185.200000
Name: height, dtype: float64

    

 出力の意味は次のとおりです。

 標準偏差(Standard Deviation:SD)は数値のばらつきの度合いを表すもので、
この値が大きいほどばらつきが大きいことになります。

 標本標準偏差と不偏標準偏差がありますが、
describe() は不偏標準偏差の方を返します。

 不偏標準偏差は、標本そのものの標準偏差ではなく
その母集団に係る標準偏差の推定値です。

 25%, 50%, 75% はパーセンタイルの値で、
最小値から最大値までを小さい順で並べたとき、
25%の値(158.675)以下の人が25%いて、
50%の値(163.6)以下の人が半数いるという意味になります。

 50%の値は、集団をちょうど半々に分けるので中央値といいます。

 中央値は、極端な値に左右されにくいというのが特徴です。

 たとえば、10人のうち9人の貯蓄額が100万円、
1人だけが1千万円だったとすると
平均値は 190万円になりますが、中央値は 100万円です。

 最高額が1千万円でなく1億円だとしても中央値は変化しません。

    

(2) 数値以外のデータへの適用

 data01.csv でいえば、gender の要素は m:男性, f:女性です。
数値ではありません。

 このようなデータにも describe() を適用できます。

 4個の要素からなるSeries ['a', 'a', 'b', 'c'] があるとします。

 aが2個、b, c がそれぞれ1個です。

 これに describe を適用すると下の結果を得ます。

count     4
unique    3
top       a
freq      2

 最頻値が複数ある場合(要素個数が同じ種類が他にもある場合)は
どれか1種類だけが表示されます。

    

(3) describe() のオプション指定

 describe() では percentiles, include, exclude のオプションを指定できます。

 includeは DataFrame の中のどの列を対象にするかをDataTypeで指定するもので
describe(include=['object']) とすれば object型が対象になります。
つまり、int, float といった数値型は対象外になります。

 excludeは、対象外にするDataTypeを指定します。

 include, exclude を省略すると、数値型の列が対象になります。

    

 percentilesは、どのパーセンタイルを得るかを指定するものです。

 デフォルト値が [.25, .5, .75] なので 25%, 50%, 75% が出力されます。

 describe(percentiles=[]) のように空のリストを指定しても
50%(中央値)は出力されるようです(中央値は必ず出力される)。

 describe(percentiles=[.8]) の場合、50%, 80% が出力されます。

目次に戻る


2. Prastクラスのdescribeメソッド

 Prastクラスのdescribeメソッドも同じ分布情報を返します。

 異なる点は、指定により分類処理を行う点です。

Prast.describe(clist, gkey=None, ddof=True, percentiles=[.25, .5, .75])
dod = psx.describe("height")
dod = psx.describe("height", "gender")

 戻り値のdodは OrderedDict です。

 dod["身長"] に分布の要約がDataFrameの形で入ります。

 heightの別名(身長)がキーになります。

 第1引数だけを指定したときは(heightだけを指定)
身長の分布情報を返します。

 describe("height", "gender") のように第2引数も指定した場合は
身長を男女別に分類した上で、それぞれの分布を1つのDataFrameにして返します。

  度数 平均 std 最小値 25% 50% 75% 最大値
男性 198 166.7 6.22 151.6 162.7 166.8 170.8 185.2
女性 192 160.1 6.37 144.8 156.1 160.1 164.4 179.8
全体 390 163.5 7.11 144.8 158.6 163.6 168.1 185.2

 第1引数を ["height", "weight"] のように
リストの形で複数 指定することもできます。

 その場合、戻り値は dod["身長"], dod["体重"] のように参照できます。

 オプションとして `ddof=False’ を指定すると、
標準偏差が不偏標準偏差でなく標本標準偏差になります。

 percentilesオプションの指定方法は、pandas.describe() と同じです。

    

 身長と体重の分布を男女別に表示するスクリプトを掲げておきます。

 1# distrib01.py (coding: cp932)
 2import pandas as pd
 3from prast import Prast
 4
 5dtf = pd.read_csv("data01.csv")
 6psx = Prast(dtf, "data01_c.txt", "data01_i.txt", "cp932")
 7dod = psx.describe(["height", "weight"], "gender")
 8for key in dod.keys():
 9    print(key)
10    print(dod[key])
11    print('')

 なお、「度数」とか「最小値」などの日本語名でなく
count, min などのままがよければ describe() を呼び出す前に
psx.etoj1 = None という1行を置いてください。

 数値データ以外の場合は psx.etoj2 = None です。

 etoj1, etoj2 は辞書データとして定義されています。

目次に戻る


3. 小数点以下の桁数の調整

 Prastクラスのメンバー変数に dround というのがあります。辞書型です。

 describe() が返す値の小数点以下の桁数を調整するものです。

 デフォルト値は次のとおり。

dround = {'count': 0, 'mean': 1, 'std': 2, 'min': 1,
    '%': 1, 'max': 1}

 countが0となっているので、これは整数値になります。

 stdが2、それ以外は1です。

 標準偏差の小数点以下が2桁、それ以外の最小値などは1桁となるわけです。

 パーセンタイルの桁数は dround['%'] で一括して設定します。

psx.dround['mean'] = 2  # 平均値
psx.dround['std'] = 3  # 標準偏差
dod = psx.describe("height", "gender")

 上のようにすると、平均値が2桁、標準偏差が3桁になります。

 zip圧縮ファイルに入っている distrib02.py がこのサンプルになっています。

 なお、psx.dround = {} とか psx.dround = None とすれば
桁数の調整は行われません。

目次に戻る


4. 正規分布か否かの確認

 分布が正規分布か否かを確認する方法の1つに shapiro test があります。

 pythonのscipyライブライにある shapiro() でその検定が可能です。

import scipy.stats as stats
w, p = stats.shapiro(x)

 戻り値は統計量w, 有意確率pの2つです。

 p<0.05 だと「正規分布である」という仮説が棄却されます。

 Prastクラスにも shapiro() を設けました。

Prast.shapiro(clist, gkey=None)
dod = psx.shapiro(["height", "weight"])
dod = psx.shapiro(["height", "weight"], "gender")

 第1引数 clist で指定された列ラベルの分布をチェックします。

 第2引数は、あってもなくてもかまいません。

 指定すると、その列を分類の手がかりにします。

 gender を指定した場合は男女別に分類した上で分布をチェック。

 戻り値 dod は OrderedDict ですが、第2引数がない場合は
DataFrameが dod['shapiro'] にセットされます。内容は次のとおり。

shapiro ||w|p_value|正規分布 |身長|0.997316|0.777182|yes |体重|0.980778|4.35878e-05|no

 第2引数として "gender" を指定した場合は
dod["身長"], dod["体重"] の2つのDataFrameが得られます。

 それぞれ男女別に分類されて正規性が検証されます。

 身長のDataFrameは下のとおり。

  w p_value 正規分布
男性 0.994998 0.755965 yes
女性 0.993902 0.617443 yes
全体 0.997338 0.785926 yes

 標本サイズがごく少数であるなどの理由で shapiro test が行えない場合は
w, p_value, 正規分布の3つとも NaN になります。

 shapiro test のサンプルスクリプトを掲げておきます。

 身長と体重について、男女別にその分布の正規性を確認。

 1# distrib03.py (coding: cp932)
 2import pandas as pd
 3from prast import Prast
 4
 5dtf = pd.read_csv("data01.csv")
 6psx = Prast(dtf, "data01_c.txt", "data01_i.txt", "cp932")
 7dod = psx.shapiro(["height", "weight"], "gender")
 8for key in dod.keys():
 9    print(key)
10    print(dod[key])
11    print('')

目次に戻る


5. 日時データのdescribe

 日時データ(datetime)を describe() にかけた場合、
数値としてではなく object型に準ずる扱いになるようです。

 平均値、標準偏差、パーセンタイルなどは算出されず
度数、種類数、最頻値、最頻度数が出力されます。

 このほかに first, last の2つも示されます。

 firstが最も初期の日時、lastが最後の日時です。

 日時データは、特殊とはいえ、内部では数値として保持されているはず。

 2つの日時の差を計算できたりするのは数値だからです。

 ということで、妙な意地を出して datetime の平均値や中央値などを
describeで表示させるよう試みました。

 あまり意味がないような気がしますが、番外編ということで……

 方法は簡単で、describe() にかける前に
UNIX時間(基準時刻からの秒数)に変換し、
得られた結果をdatetimeに復元するだけです。

 UNIX時間は、世界標準時の 1970-01-01 00:00:00 を始点として
そこからの経過秒数を意味するようです。

 日本時間は9時間先をいくので同日の 09:00 が始点ということになるでしょうか。

    

 サンプルのスクリプトを掲げます。

 ID, birthday の2列×100人分の data02.csv を処理します。

 birthday は日付だけで時刻はありません。

 1# distrib04.py (encoding: cp932)
 2import pandas as pd
 3from datetime import datetime
 4from prast import Prast
 5
 6ckey = "birthday"  # 処理対象(datetime型)の列名
 7dtf = pd.read_csv("data02.csv", parse_dates=[ckey])
 8    # ↓ datetimeを基準時点からの経過秒数に変換
 9dtf[ckey] = dtf[ckey].map(pd.Timestamp.timestamp).astype(int).apply(
10    lambda x: x - 60*60*9)
11psx = Prast(dtf)
12for key in psx.dround.keys():  # mean, std などをすべて整数値に
13    psx.dround[key] = 0
14dod = psx.describe(ckey, percentiles=[.5])
15for key in dod.keys():
16    dtfx = dod[key]
17    for cname in dtfx.columns:
18        if cname in ["count", u"度数"]:
19            continue
20        dtfx[cname] = dtfx[cname].map(datetime.fromtimestamp)
21    dtfx = dtfx.T  # 行と列を入れ替え
22    print(dtfx)

    

 csvファイルを単純に読み込むと birthday が datetime でなく object型になります。

 そこで、read_csv(……, parse_dates=["birthday"]) として
2列目のbirthdayを datetime として読み込んでいます。

 parse_dates=[1] のように番号を書いてもいいようです。

 datetime→UNIX時間への変換は pd.Timestamp.timestamp(x) で行います。

 UNIX時間→datetimeの変換は datetime.fromtimestamp(x) です。

 スクリプトが出力する結果は下のとおり。

                      birthday
度数                         100
平均         1992-07-21 21:07:12
std        1972-06-10 00:23:19
最小値        1988-08-25 00:00:00
50%        1992-08-09 12:00:00
最大値        1996-12-18 00:00:00

 標準偏差は 1972念6月ですから、1970念1月からの経過期間をみて
約2年半のインターバルであることがわかります。

 中央値が 1992-08-09 なので、これより前に生まれた人が半数、
後に生まれた人も半数ということになります。

 スクリプトでは1秒未満のミリ秒は切り捨てるようにしました。

 また、日本時間と世界標準時間のずれ(9時間)を調整するため
UNIX時間から 60*60*9 を引き算していますが、
これは環境によっては余計なお世話ですね。

〜 以上 〜

Copyright (C) T. Yoshiizumi, 2018 All rights reserved.