データ分析に取り掛かるときに、データセットの概要を数値で表現できると便利なことがあります。
集計軸などは対象とするデータセットによって異なりますが、基本事項として要約統計量を計算しておくと、共通用語の下でデータの傾向を把握でき、作業の手戻りが少なくなることを期待できます。
Wikipedia 日本語版の要約統計量のページでは３つの統計量について説明があります。
実際のデータに対して pandas (Python Data Analysis Library) を使ってこれらの値を計算してみます。
- モーメントから求められる要約統計量
- 順序から求められる要約統計量
- 度数から求められる要約統計量
データの用意
今回は世界銀行の API (data.worldbank.org) から各国の年ごとの GDP のデータを使います。
世界銀行は Indicators、Projects、World Bank financial data という３種類の API を提供しており、For Developers のページに説明があります。Indicators は時系列データで、8,000を超えるデータセットがあります。その中に NY.GDP.MKTP.CD というコードで表現されるデータがあります。API の使い方は Indicator Queries のページに説明がありますが、実際にその API を使って構築されているページが GDP at market prices (current US$) です。
なお、似たものとしては NY.GDP.PCAP.KD があり、こちらは一人当たり GDP を表します。 - GDP per capita (constant 2005 US$)
API のパラメータを調整
API レスポンスのデータ形式は JSON, XML の両方が提供されていますので、結果を処理するプログラムを書けば良いのですが、 pandas においては pandas_datareader がこの機能を提供してくれます。 pandas_datareader は元々は pandas.io としてコア機能に含まれていたものですが、外部モジュールに切り出されたものです。既存のスクリプトで pandas をバージョンアップして警告文が表示される場合には、依存ライブラリを追加し、 import 文を修正しましょう。
それでは、実際にデータを読み込んでみます。
from pandas_datareader import wb
df = wb.download(indicator='NY.GDP.MKTP.CD', start=1960, end=2014)
df.count()
NY.GDP.MKTP.CD 165
dtype: int64
165件あるようです。国の数としては少ないため、 head() で先頭の数件を表示してみます。
df.head()
| NY.GDP.MKTP.CD | ||
|---|---|---|
| country | year | |
| Canada | 2014 | 1.785387e+12 |
| 2013 | 1.838964e+12 | |
| 2012 | 1.832716e+12 | |
| 2011 | 1.788796e+12 | |
| 2010 | 1.614014e+12 |
マルチインデクスになっており、上位が国、下位が年度を表しています。
インデクスの中身は .index 属性で確認できます。
print(df.index.names)
print(df.index.levels[0])
['country', 'year']
Index(['Canada', 'Mexico', 'United States'], dtype='object', name='country')
国はカナダ、メキシコ、アメリカの３か国になっています。これは、ソースコードの実装で pandas_datareader.wb のデフォルトがこの３つになっているためです。各国の GDP を比較したいので、リクエスト条件を調整することにします。
wb.download は country 引数を受け付けますので、ここに任意のコードを与えます。 all で全てのコードを指定できるのですが、集約コード (Aggregates - Regions and Income Levels | Data) も含まれてしまいますので、Country Queries を使ってフィルタします。
get_countries() 関数で一覧を取得できますので、その結果を見てフィルタ条件を決めます。
df_countries = wb.get_countries()
df_countries.head()
| | adminregion | capitalCity | iso3c | incomeLevel | iso2c | latitude | lendingType | longitude | name | region |
|---|---|------------|-----|----------------------|----|----------|--------|--------|--------|--------|----------|--------|
| 0 | | Oranjestad | ABW | High income: nonOECD | AW | 12.51670 | Not classified | -70.0167 | Aruba | Latin America & Caribbean (all income levels) |
| 1 | South Asia | Kabul | AFG | Low income | AF | 34.52280 | IDA | 69.1761 | Afghanistan | South Asia |
| 2 | | | AFR | Aggregates | A9 | NaN | Aggregates | NaN | Africa | Aggregates |
| 3 | Sub-Saharan Africa (developing only) | Luanda | AGO | Upper middle income | AO | -8.81155 | IBRD | 13.2420 | Angola | Sub-Saharan Africa (all income levels) |
| 4 | Europe & Central Asia (developing only) | Tirane | ALB | Upper middle income | AL | 41.33170 | IBRD | 19.8172 | Albania | Europe & Central Asia (all income levels) |
インデクス番号２のレコードで incomeLevel、lendingType, region の３つが Aggregates になっています。これはアフリカ全体を集約した結果のレコードです。
レコードの出現頻度を確認するために region で集計してみます。
df_countries.groupby('region')['name',].count().reset_index()
| region | name | |
|---|---|---|
| 0 | Aggregates | 49 |
| 1 | East Asia & Pacific (all income levels) | 37 |
| 2 | Europe & Central Asia (all income levels) | 57 |
| 3 | Latin America & Caribbean (all income levels) | 41 |
| 4 | Middle East & North Africa (all income levels) | 21 |
| 5 | North America | 3 |
| 6 | South Asia | 8 |
| 7 | Sub-Saharan Africa (all income levels) | 48 |
集約レコード (Aggregates) が49件あることが分かりましたので、これを除外して国コードを構築することにします。結果を確認すると、インデクス番号２のレコードが除外されていることが分かります。
df_countries[(df_countries['region'] != 'Aggregates')]['iso2c'].head()
0 AW
1 AF
3 AO
4 AL
5 AD
Name: iso2c, dtype: object
対象とするデータセットを取得
対象とする国コードの一覧ができましたので、 country 引数にリストとして与えることで API を呼び出します。また、データ期間は2014年の一年間に限定します。
なお、ソースコード内にハードコードされた国コードに含まれないコードも存在するため、警告メッセージが表示されます。ハードコードされた一覧と最新の一覧の不一致により発生する事象ですので、処理はそのまま続けて問題ありません。
df = wb.download(indicator='NY.GDP.MKTP.CD',
country=df_countries[(df_countries['region'] != 'Aggregates')]['iso2c'].tolist(),
start=2014, end=2014)
df.head()
| NY.GDP.MKTP.CD | ||
|---|---|---|
| country | year | |
| Aruba | 2014 | NaN |
| Andorra | 2014 | NaN |
| Afghanistan | 2014 | 2.003822e+10 |
| Angola | 2014 | NaN |
| Albania | 2014 | 1.321151e+10 |
このままだと少し扱い難いため、データフレームを少し変形します。
d = df.unstack()['NY.GDP.MKTP.CD']
dd = pd.DataFrame(data={'gdp_2014': d['2014']}, index=d.index).reset_index()
dd.head()
| country | gdp_2014 | |
|---|---|---|
| 0 | Afghanistan | 2.003822e+10 |
| 1 | Albania | 1.321151e+10 |
| 2 | Algeria | 2.135185e+11 |
| 3 | American Samoa | NaN |
| 4 | Andorra | NaN |
データ件数を確認すると、全部で214件あり、gdp_2014 が NaN で無いものは185件あることが分かります。このデータについて、要約統計量を計算していきます。
dd.count()
country 214
gdp_2014 185
dtype: int64
モーメントから求められる要約統計量
pandas の DataFrame には200近いメソッドが定義されていますが、この中にモーメントを計算する関数も含まれています。
- 平均:
mean() - 分散:
std() - 歪度:
skew() - 尖度:
kurtosis()
ds = pd.DataFrame(index=('gdp_2014',))
ds['count'] = dd.count()['gdp_2014']
ds['mean'] = dd.mean()['gdp_2014']
ds['std'] = dd.std()['gdp_2014']
ds['skew'] = dd.skew()['gdp_2014']
ds['kurtosis'] = dd.kurtosis()['gdp_2014']
ds
| count | mean | std | skew | kurtosis | |
|---|---|---|---|---|---|
| gdp_2014 | 185 | 4.138522e+11 | 1.595941e+12 | 8.197562 | 78.277232 |
年ごとに計算を繰り返すことで、分布の変化を追跡できるようになります。
wb.download() を呼び出すときの start パラメータの値を変えて、データ変換を再実行しながら試してみましょう。
DataFrame のその他のメソッドは API ドキュメントを参照しましょう。
順序から求められる要約統計量
pandas の DataFrame には describe() メソッドがあります。
先ほど計算したデータ件数、平均、分散に加えて、最小、最大、四分位数を計算します。四分位数の50%は中央値になります。
dd.describe()
| gdp_2014 | |
|---|---|
| count | 1.850000e+02 |
| mean | 4.138522e+11 |
| std | 1.595941e+12 |
| min | 3.785955e+07 |
| 25% | 7.962424e+09 |
| 50% | 3.299619e+10 |
| 75% | 2.178723e+11 |
| max | 1.741900e+13 |
最小と最大は6桁も異なりますので、ヒストグラムを描いたり、刈込平均を計算してみると、分布の様子をより理解出来るようになるでしょう。
ここでは、上位と下位のそれぞれ10件を棒グラフで描画してみます。
描画には Seaborn を使います。
import seaborn as sns
%matplotlib inline
sns.set_style('whitegrid')
上位10件を表示します。グラフにしてみると、上位2件（アメリカと中国）が突出していることが分かります。
sns.barplot(x='gdp_2014', y='country',
data=dd.sort_values('gdp_2014', ascending=False).head(10).reset_index(),
palette=sns.color_palette('Greens_r', n_colors=15))
同様にして下位10件を表示します。下位1件（ツバル）は少し上位よりも差が大きいことが分かります。
sns.barplot(x='gdp_2014', y='country',
data=dd.sort_values('gdp_2014').head(10).reset_index(),
palette=sns.color_palette('Reds_r', n_colors=15))
度数から求められる要約統計量
上位、下位ともに分布が大きく散らばっていることが分かりましたので、ヒストグラムを描きながら度数の多い範囲を探してみます。
とりあえずでヒストグラムを描画してみると、上位2件の値に引きづられて特徴の分からないグラフになってしまいます。
dd.hist(bins=30)
桁が大きく異なることが一因ですので、log10 をとってみます。
import math
lg10 = dd['gdp_2014'].map(lambda v: math.log10(v))
lg10.describe()
count 185.000000
mean 10.533801
std 1.048719
min 7.578175
25% 9.901045
50% 10.518464
75% 11.338202
max 13.241023
Name: gdp_2014, dtype: float64
lg10.hist(bins=30)
今度は先ほど確認した特徴が見えてきました。
下位1件が左側、上位2件が右側に外れ値として存在しています。また、log10 が10から11の部分に多くの値が含まれています。もう少し正確に四分位数も確認すると、9.9から11.3に約半数が含まれていると言えます。
散布図と組み合わせてみると、傾向をよりはっきりと確認できます。
d = dd['gdp_2014']
a = pd.DataFrame({'gdp_2014': d, 'log10': d.map(lambda v: math.log10(v))})
sns.jointplot(x='log10', y='gdp_2014', data=a, marginal_kws=dict(bins=30))
外れ値を除外してみると、もう少しグラフの値が読みやすくなり、log10 が 12 付近での非連続性が分かります。
d = dd[(lg10 > 8) & (lg10 < 13)]['gdp_2014']
a = pd.DataFrame({'gdp_2014': d, 'log10': d.map(lambda v: math.log10(v))})
sns.jointplot(x='log10', y='gdp_2014', data=a, marginal_kws=dict(bins=30))
また、中央値である10.5を中心として±1.5に範囲を限定すると、9.5付近と11.0付近の値が少ないことが分かります。
d = dd[(lg10 >= 9) & (lg10 <= 12)]['gdp_2014']
a = pd.DataFrame({'gdp_2014': d, 'log10': d.map(lambda v: math.log10(v))})
sns.jointplot(x='log10', y='gdp_2014', data=a, marginal_kws=dict(bins=30))
さて、最初に確認した Country Query のレスポンスには incomeLevel という属性が含まれていました。これを改めて集計してみると、以下の結果を得られます。Aggregates は集約結果なので別扱いとして、2種類の上位層、1種類の下位層、そして2種類の中間層になっています。ここまでで確認した度数密度と合わせて考えると、具体的な数値をイメージしやすくなったと思います。
df_countries.groupby('incomeLevel')['name',].count().reset_index()
| incomeLevel | name | |
|---|---|---|
| 0 | Aggregates | 49 |
| 1 | High income: OECD | 32 |
| 2 | High income: nonOECD | 48 |
| 3 | Low income | 31 |
| 4 | Lower middle income | 51 |
| 5 | Upper middle income | 53 |
追記 (2016/7/2)：
７月１日に世界銀行での分類が更新されました。いくつかのグループが変わっています。
分類の詳細としては、国民一人あたり総所得 (gross national income per capita) が $1,025 以下だと Low income 、 $1,026 から $4,035 は Lower middle-income 、 $4,036 から $12,475 は Upper middle-income 、 $12,476 以上だと High-income になるそうです。
/追記ここまで
まとめ
pandas を使って、モーメント、順序、度数に関する要約統計量を計算しました。
各国のGDPのように、正規分布などの特定の分布を想定することが難しい場合でも、個別の値を見たり対数を計算することで特徴を数値で表現出来るようになりました。その中では、ヒストグラム (histogram) と散布図 (scatter plot) を使い、グラフでデータ分布を確認することも有用です。
実際の利用シーンについては、PayPal の人が書いた記事（日本語訳がありがたい）が分かりやすいと思います。基本的な数値計算はライブラリがカバーしてくれますので、その数値を読んで意味を理解できるようになると、データ分析の内容にも深みが出てくると考えられます。