カテゴリー: プログラミング

ベクトル自己回帰のダミーデータを生成するために、多変量正規分布に従う乱数が必要なので、
Pythonで生成する方法を紹介します。

numpyとscipyにそれぞれ用意されています。
同じ名前の関数だったので、どちらかの実装をもう一方がラップしているのかと思っていたのですが、
引数の微妙な違いなどあり、どうやら個別に実装されているようです。

ドキュメントはそれぞれ次のページにあります。

numpy
numpy.random.multivariate_normal
(この記事の numpy は version 1.16を使っています。 numpy 1.17.0 のリリースノートを見ると、random moduleに変更が加えられており、どうやらこの関数にも影響が出てるようなのでご注意ください。)

scipy
scipy.stats.multivariate_normal

さて、実際に期待値と分散共分散行列を指定してそれぞれ乱数を生成してみましょう。

import numpy as np
from scipy.stats import multivariate_normal

# 期待値と分散共分散行列の準備
mean = np.array([3, 5])
cov = np.array([[4, -1.2], [-1.2, 1]])

# numpy を用いた生成
data_1 = np.random.multivariate_normal(mean, cov, size=200)

# データ型の確認
print(data_1.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_1, axis=0))
# [3.00496708 4.94669956]

# 分散共分散の確認
print(np.cov(data_1, rowvar=False))
"""
[[ 3.86542859 -1.31389501]
 [-1.31389501  0.93002097]]
"""

# scipyで生成する方法
data_2 = multivariate_normal(mean, cov).rvs(size=200)

# データ型の確認
print(data_2.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_2, axis=0))
# [2.81459692 5.10444347]

# 分散共分散の確認
print(np.cov(data_2, rowvar=False))
"""
[[ 4.46151626 -1.28084696]
 [-1.28084696  1.06831954]]
"""

それぞれきちんと生成できたようです。

分散共分散行列の正定値性のバリデーションなど細かなオプションを持っていますが、
あまり使う機会はなさそうです。
(きちんと行う場合も、事前に固有値を求めて確認しておけば大丈夫だと思います。)

前回の更新でPrestoでデータの縦横変換をする方法を紹介しましたが、
クエリで処理を完結させる必要がないときは、一旦pandasのデータフレームに格納してから処理をするのも便利です。

その場合、 pandas.pivot_table を使います。

使い方は簡単で、pd.pivot_tableに、変換したいデータフレーム、
列にするカラム、行にするカラム、集計する値、集計に使う関数を指定するだけです。
fill_value引数で欠損値を埋めるなどの細かい設定もできます。
ドキュメントの例を使ってやってみます。

import pandas as pd
# データ作成
df = pd.DataFrame(
    {
        "A": ["foo", "foo", "foo", "foo", "foo", "bar", "bar", "bar", "bar"],
        "B": ["one", "one", "one", "two", "two", "one", "one", "two", "two"],
        "C": ["small", "large", "large", "small",
              "small", "large", "small", "small", "large"],
        "D": [1, 2, 2, 3, 3, 4, 5, 6, 7],
        "E": [2, 4, 5, 5, 6, 6, 8, 9, 9]
    }
)
print(df)
"""
     A    B      C  D  E
0  foo  one  small  1  2
1  foo  one  large  2  4
2  foo  one  large  2  5
3  foo  two  small  3  5
4  foo  two  small  3  6
5  bar  one  large  4  6
6  bar  one  small  5  8
7  bar  two  small  6  9
8  bar  two  large  7  9
"""

table_0 = pd.pivot_table(
                df,
                values="D",
                index="A",
                columns="C",
                aggfunc="sum",
                fill_value=0,
        )
print(table_0)
"""
C    large  small
A
bar     11     11
foo      4      7
"""

# 行や列、集計関数は配列で複数指定することもできる
table_1 = pd.pivot_table(
                df,
                values="D",
                index=["A", "B"],
                columns="C",
                aggfunc=["sum", "count"],
                fill_value=0,
        )
print(table_1)
"""
          sum       count
C       large small large small
A   B
bar one     4     5     1     1
    two     7     6     1     1
foo one     4     1     2     1
    two     0     6     0     2
"""

昨日に続いて例外処理の話です。
ある関数内に、エラーが発生しうる処理がある時、その関数内でtry:〜except:〜処理を書いて
綺麗に例外を処理するけど、その例外を関数の呼び出し元にも伝えて例外を伝播させたいことがあります。

このような時も、 raise 文を使うことができます。

自分は最近まで raise 文は新規に例外を発生させいさせる機能しかないと思ってました。
こういう風に。

raise ValueError

"""
ValueError                                Traceback (most recent call last)
 in ()
----> 1 raise ValueError

ValueError: 
"""

しかし、raiseを単体で使用すると、そのスコープで有効になっている例外を再送出できます。

参考: 7.8. raise 文

試してみる前に、非常に単純な例なのですが次のようなケースを考えてみます。
関数 inv は 引数の逆数を返す関数で、0が渡されたら例外になるはずのものです。
そして、0が渡されたら内部で例外処理をしています。
そして、print_inv は inv を使って、与えられた数の逆数を表示します。

def inv(x):
    try:
        return 1/x
    except Exception as e:
        print(e)


def print_inv(x):
    try:
        print(inv(x))
    except Exception as e:
        print(e)
    else:
        print("print_invで例外は発生しませんでした")


print_inv(0)
"""
division by zero
None
print_invで例外は発生しませんでした
"""

ご覧の通り、 inv 内で例外処理しているので、呼び出し元では例外の発生を検知できていません。

ここで、 raiseを使って例外の再送出を入れてみます。

def inv(x):
    try:
        return 1/x
    except Exception as e:
        print(e)
        # 例外の再送出
        raise


def print_inv(x):
    try:
        print(inv(x))
    except Exception as e:
        print(e)
    else:
        print("print_invで例外は発生しませんでした")


print_inv(0)
"""
division by zero
division by zero
"""

division by zero が 2回表示されました。 inv と print_inv でそれぞれキャッチされた例外です。

jupyterでインタラクティブにPythonを使っているとあまり必要ないのですが、
本番コードを書くときなどは流石に例外処理を真面目に実装する必要があることがあります。
そこまで高頻度にあることではなく、すぐ忘れてしまうので、書き方をまとめておこうと思います。

参考になるドキュメントは次の2箇所です。
8. エラーと例外
組み込み例外

基本的に次のような書き方になります。
必須なのは、 try と except で、 exceptは複数書くこともできます。
except する例外には as e のように別名をつけることができ、
別名をつけておけば処理中で利用できます。
else と finally はオプションなので不要ならば省略可能です。

try:
    # ここに例外が発生しうるコードを書く

except [キャッチしたい例外クラス]:
    # 例外が発生した時に実行するコード

else:
    # 例外が発生なかった時に実行するコード

finally:
    # 必ず実行するコード

とりあえず定番の 0で割る演算で試してみましょう。

import numpy as np


def inv(data):
    try:
        inverse_data = 1/data
    except ZeroDivisionError as e:
        print(e)

    except TypeError as e:
        print(e)

    else:
        print("正常終了")
        return inverse_data
        print("このメッセージは表示されない")

    finally:
        print("finallyに書いた文は必ず実行されます")


print(inv(5))
"""
正常終了
finallyに書いた文は必ず実行されます
0.2
"""

print(inv(0))
"""
division by zero
finallyに書いた文は必ず実行されます
None
"""

print(inv("a"))
"""
unsupported operand type(s) for /: 'int' and 'str'
finallyに書いた文は必ず実行されます
None
"""

例外が発生した、0と”a”　については想定通りに動きました。

実は例外が発生しなかったinv(5)が僕にとっては少し驚きでした。
else: のブロックの中で、 return して関数を抜けているので、
それより後ろの finally: のブロックは流石に実行されないと思っていたのですが、
print関数がバッチリ実行されています。

改めてよく読んでみれば、ドキュメント中にもしっかりそう書いてありました。
この辺りはきちんと理解して使う必要がありそうです。

– もし try 文が break 文、 continue 文または return 文のいずれかに達すると、その:keyword:break 文、 continue 文または return 文の実行の直前に finally 節が実行されます。
– もし finally 節が return 文を含む場合、 try 節の return 文より先に、そしてその代わりに、 finally 節の return 文が実行されます。

今回はブログ記事用に書いたコードだったので、
ZeroDivisionError と TypeError を 分けて書きましたが、
Exception のようなキャッチできる範囲の広い例外を指定しておけばまとめて受け取ってくれます。
(本当はあまり良くないと思うのですが、便利なので大抵そうしています。)

def inv2(data):
    try:
        inverse_data = 1/data
        return inverse_data
    except Exception as e:
        print(e)
        return None


print(inv2(5))
"""
0.2
"""

print(inv2(0))
"""
division by zero
None
"""

print(inv2("a"))
"""
unsupported operand type(s) for /: 'int' and 'str'
None
"""

また、例外処理の中で例外の種類を区別する必要が全くない場合、
except Exception as e:
の代わりに、
except:
とだけ書いておけば、より簡単に全ての例外をキャッチしてくれます。

昨日の indicator の記事を書くためにドキュメントを読んでいて見つけた、 validate という引数の紹介です。

データフレーム通しを結合する時に、結合に使うキーのユニーク性が重要になることがあります。
事前に確認するようにコードを書いておけば済む話ではあるのですが、
pandasのmerge関数では、キーが一意でなかった時にエラーを上げてくれる引数があるようです。

ドキュメント: pandas.merge

引数 validate には、 one_to_one / one_to_many / many_to_one / many_to_many の4種類の文字列か、None(デフォルト)を
渡すことができます。
そして、 one_to_one なら 1:1, one_to_many なら 1:m　、という風にkeyが対応してなければエラーを上げてくれます。
many_to_many と None はノーチェックです。
なので、きちんとエラーをキャッチするように次のようにコードを書けます。
(データの準備等は省略)

try:
    df_merge = pd.merge(
            df_0,
            df_1,
            on="key",
            how="outer",
            validate="one_to_many",
        )
except pd.errors.MergeError as e:
    print(e)

# Merge keys are not unique in left dataset; not a one-to-many merge

発生するエラーは pd.errors.MergeError です。
データを何種類か用意して、validateの引数を変えながら動かすと色々動きがわかると思います。

どこで見かけたか忘れてしまった(TwitterかQiitaかその辺りのはず)のですが、
pandasのデータフレームのマージをする時に便利な引数を知ったので紹介します。

DataFrame同士を列の値で結合する時、pd.mergeを使います。

how=”inner”で利用する場合は何も問題ないのですが、
left/right/outerで使う場合、結果の中に、ちゃんと左右のデータフレームにレコードが存在してうまく結合できた行と、
一方にしか存在せず、結合はしなかった行が混在します。

left_on/right_on　を使って結合した場合はそこの欠損を見ればまだ見分けられるのですが、
同名列をonで結合すると見分けがつかず、少し不便です。

このような時、 indicator=True を指定しておくと、 結果に _merge という列が追加され、
各レコードが左右のデータフレームのどちらに起因しているか出力してくれます。

やってみたのがこちらです。

import pandas as pd
df_0 = pd.DataFrame(
            {
                "id": range(5),
                "key": [1, 5, 12, 7, 8],
                "value0": ["a", "b", "c", "d", "e"],
            }
        )
df_1 = pd.DataFrame(
            {
                "key": range(10),
                "value1": range(0, 100, 10),
            }
        )

df_merge = pd.merge(
        df_0,
        df_1,
        on='key',
        how="outer",
        indicator=True,
    )
print(df_merge)
"""
     id  key value0  value1      _merge
0   0.0    1      a    10.0        both
1   1.0    5      b    50.0        both
2   2.0   12      c     NaN   left_only
3   3.0    7      d    70.0        both
4   4.0    8      e    80.0        both
5   NaN    0    NaN     0.0  right_only
6   NaN    2    NaN    20.0  right_only
7   NaN    3    NaN    30.0  right_only
8   NaN    4    NaN    40.0  right_only
9   NaN    6    NaN    60.0  right_only
10  NaN    9    NaN    90.0  right_only
"""

both / left_only / right_only
で、 key の由来が確認できます。

単に何かの領域を塗りつぶしたり、時系列データの予測モデルの信頼区間の可視化などで使われたり、
関数のグラフとx軸の間を塗りつぶしたりするあいつです。

matplotlibでは、fill_between というメソッドが用意されており、これを使って実現できます。
ドキュメント: matplotlib.axes.Axes.fill_between

通常の plot は xとyの値をリストか何かで渡しますが、fill_betweenでは、y1とy2という風にyの値を２ペア渡します。
(なお、y2を省略すると、y1とx軸の間を塗りつぶしてくれます。)

また、 y1 と y2 の間を全て塗りつぶすのではなく、 where で、塗りつぶす領域を指定することもできます。
where に渡すのは x と同じ長さの True or False のリストです。
TrueとTrueの間が塗りつぶされます。
False, True, False のような孤立したTrueの分は塗りつぶされないので注意が必要です。

この他、 interpolate という引数が用意されています。
これは where が使われていて、かつ二つの曲線が閉じている場合に、はみ出さないように綺麗に塗ってくれるオプションです。
とりあえずTrue指定しておいて良いと思います。
この後サンプルコードを紹介しますが、最後の一つのグラフはあえて interpolate を指定せずに少しガタついてるグラフにしました。

import numpy as np
import matplotlib.pyplot as plt

# データ作成
x = np.linspace(0, 2*np.pi, 101)
y1 = np.sin(x)
y2 = np.sin(2*x)

fig = plt.figure(figsize=(12, 8), facecolor="w")
ax = fig.add_subplot(2, 2, 1, title="2線の間を全て塗りつぶす")
ax.plot(x, y1, label="sin(x)")
ax.plot(x, y2, label="sin(2x)")
ax.fill_between(
    x,
    y1,
    y2,
    alpha=0.3,
    interpolate=True,
)
ax.legend()

ax = fig.add_subplot(2, 2, 2, title="Whereで塗りつぶす領域を絞り込む")
ax.plot(x, y1, label="sin(x)")
ax.plot(x, y2, label="sin(2x)")
ax.fill_between(
    x,
    y1,
    y2,
    where=(y1 >= y2),
    alpha=0.3,
    interpolate=True,
    label="sin(x)>=sin(2x)"
)
ax.fill_between(
    x,
    y1,
    y2,
    where=(y1 < y2),
    alpha=0.3,
    interpolate=True,
    label="sin(x)<sin(2x)"
)
ax.legend()

ax = fig.add_subplot(2, 2, 3, title="y2を省略するとx軸との間を塗りつぶす")
ax.plot(x, y1, label="sin(x)")
ax.fill_between(
    x,
    y1,
    alpha=0.3,
    interpolate=True,
)
ax.legend()

ax = fig.add_subplot(2, 2, 4, title="interpolate=Trueを指定しないと隙間が発生しうる")
ax.plot(x, y1, label="sin(x)")
ax.plot(x, y2, label="sin(2x)")
ax.fill_between(
    x[::10],
    y1[::10],
    y2[::10],
    alpha=0.3,
)
ax.legend()

plt.show()

結果。

実用性は皆無なのですが、他にやっている人を見かけなかったのでやってみました。
前回の記事で紹介した matplotlibのxkcdスタイルの続きです。
ドキュメントを読めば明らかなのですが、 plt.xkcd()には３種類の引数を渡すことができます。

ドキュメント: matplotlib.pyplot.xkcd

３つの引数と意味はそのまま引用します。

scale : float, optional
The amplitude of the wiggle perpendicular to the source line.

length : float, optional
The length of the wiggle along the line.

randomness : float, optional
The scale factor by which the length is shrunken or expanded.

初期値は (scale=1, length=100, randomness=2) です。

色々試したところ、 scaleと randomness は 増やすと徐々にグラフが崩れていき、
length は減らすと崩れていくようです。

初期値と、それぞれ値を変更した3パターンをグラフ出力してみました。
(randomness はこれだけ変えても変化がわかりにくかったので、scaleも変更しています。)

import matplotlib.pyplot as plt
import numpy as np


# グラフを描く処理は共通化
def graph_plot(ax):
    X0 = np.linspace(0, 2*np.pi, 200)
    Y_sin = np.sin(X0)+2
    Y_cos = np.cos(X0)+2
    X1 = np.arange(7)
    Y1 = (X1 ** 2)/36
    ax.plot(X0, Y_sin, label="$y=\\sin(x)$")
    ax.plot(X0, Y_cos, label="$y=\\cos(x)$")
    ax.bar(X1, Y1, alpha=0.3, color="g")
    ax.legend()


fig = plt.figure(figsize=(12, 10), facecolor="w")
# 間隔調整
fig.subplots_adjust(hspace=0.3, wspace=0.3)
# xkcd オプションの影響を局所化するため with で使う。
with plt.xkcd(scale=1, length=100, randomness=2):
    ax = fig.add_subplot(2, 2, 1, title="default")
    graph_plot(ax)

with plt.xkcd(scale=2, length=100, randomness=2):
    ax = fig.add_subplot(2, 2, 2, title="scale=2")
    graph_plot(ax)

with plt.xkcd(scale=1, length=50, randomness=2):
    ax = fig.add_subplot(2, 2, 3, title="length = 50")
    graph_plot(ax)

with plt.xkcd(scale=2, length=100, randomness=6):
    ax = fig.add_subplot(2, 2, 4, title="scale=2, randomness=6")
    graph_plot(ax)
plt.show()

出力されるのがこちらです。

結構雰囲気変わりますね。
とはいえ、あまりやりすぎるとくどくなるので、初期設定だけで困ることもなさそうです。
(そもそもこのスタイルが必要になる場面も基本的に無いのですが。)

全くの余談ですが、matplotlibのドキュメントページのULRにxkcdをつけるとドキュメントのスタイルが変わります。
(よくみると内容も変わってっています。)

お暇な時に見比べてみてください。
https://matplotlib.org/
https://matplotlib.org/xkcd/

xkcdってなんだ？って方はこちらをどうぞ。
https://xkcd.com/
wikipedia: xkcd

誰が何の目的で実装されたのか不明ですが、matplotlibにはグラフをxkcdのコミック風に出力する機能があります。
面白いので僕はこういう機能は結構好きです。

ドキュメント: matplotlib.pyplot.xkcd

使い方は簡単で、グラグを書く前、要はplotやbarなどの関数を使う前に、plt.xkcd()を差し込むだけ。
ただ、この手軽さに落とし穴がありました。一回呼び出すと戻せなくなるのです。
(調査にかなり手こずったので、この記事もどちらかというとxkcdの使い方より元への戻し方を伝えたい。)

ドキュメントにも、pcParamsを上書きしてしまうと書いてあります。

Notes
This function works by a number of rcParams, so it will probably override others you have set before.
f you want the effects of this function to be temporary, it can be used as a context manager, for example:

context manager ってのは要は with句の事のようです。(あとでちゃんと調べたい。)

要するに、 with plt.xkcd(): で有効な範囲をあらかじめ絞りましょう
参考に、以下のコードでは二つのグラフをxkcd風に書いた後に、普通のグラフを2つ作成しました。

import matplotlib.pyplot as plt
import numpy as np

data = [5, 4, 3, 2, 1]
label = [f"item {i}" for i in range(5)]
X = np.linspace(0, 2*np.pi, 200)
Y_sin = np.sin(X)
Y_cos = np.cos(X)

fig = plt.figure(figsize=(12, 9), facecolor="w")

# xkcd オプションの影響を局所化するため with で使う。
with plt.xkcd():
    ax = fig.add_subplot(2, 2, 1, title="xkcd pie chart")
    ax.pie(
        data,
        labels=label,
        autopct='%3.1f%%',  # 割合をグラフ中に明記
        counterclock=False,  # 時計回りに変更
        startangle=90,  # 開始点の位置を変更
    )
    ax = fig.add_subplot(2, 2, 2, title="xkcd plot")
    ax.plot(X, Y_sin, label="$y=\\sin(x)$")
    ax.plot(X, Y_cos, label="$y=\\cos(x)$")
    ax.legend()

ax = fig.add_subplot(2, 2, 3, title="normal pie chart")
ax.pie(
    data,
    labels=label,
    autopct='%3.1f%%',  # 割合をグラフ中に明記
    counterclock=False,  # 時計回りに変更
    startangle=90,  # 開始点の位置を変更
)

ax = fig.add_subplot(2, 2, 4, title="normal plot")
ax.plot(X, Y_sin, label="$y=\\sin(x)$")
ax.plot(X, Y_cos, label="$y=\\cos(x)$")
ax.legend()

plt.show()

出力されたグラフがこちら。

場面を選べば、プレゼンや資料などで使いやすそうなグラフですね。
あと、注意点としては、対応したフォントがないので日本語文字が使えません。

もし、 with を使わずに plt.xkcd() してしまったら、それ以降のグラフは全部、
xkcdモードで出力されてしまいます。
jupyter notebookであれば戻すためにカーネルの再起動が必要になるので気をつけましょう。

自分の分析のために円グラフを描きたい場面というのがあまりなく、誰かのためのダッシュボードや、稀に自分でも必要になるときはTableauで作成するので、
Pythonで円グラフを作成することはあまりないのですが、最近機会があったので方法メモしておきます。

円グラフは英語で pie chart というので、
matplotlibでも pie という名前の関数で作成できます。
ドキュメント: matplotlib.pyplot.pie

最小構成で作成するなら、データとラベルを渡してあげれば、それだけで描いてくれます。
データも割合ではなく数量で渡しても、勝手に合計100%になるように描いてくれるので楽です。
ただ、初期設定だと(個人的に)少し見慣れないデザインになるのでいくつかオプション設定します。
比較用にほぼ何も設定しないバージョンとそれぞれ出力したのが次のコードです。

import matplotlib.pyplot as plt

data = [5, 4, 3, 2, 1]
label = ["項目1", "項目2", "項目3", "項目4", "項目5"]

fig = plt.figure(figsize=(12, 7), facecolor="w")
ax = fig.add_subplot(1, 2, 1, title="円グラフ (初期設定)")
ax.pie(
    data,
    labels=label,
)

ax = fig.add_subplot(1, 2, 2, title="円グラフ (見た目修正)")
ax.pie(
    data,
    labels=label,
    autopct='%3.1f%%',  # 割合をグラフ中に明記
    counterclock=False,  # 時計回りに変更
    startangle=90,  # 開始点の位置を変更
)
plt.show()

結果がこちら。

どこまで工数を使うかにもよりますが、色とかももう少し工夫した方が使いやすそうですね。
(ただ、あまり凝ったことをするなら別のBIツールに任す方がオススメです。)