pandasのデータフレームのexplodeメソッドの紹介

以前の記事で、pandasのDataFrameの列中に含まれている配列を行に展開する方法について書きました。
参考: pandasのDataFrameのappendは遅い

お恥ずかしながらこの記事を書いた時は知らなかったのですが、pandasの version 0.25.0 から、専用のメソッドが準備されていて、
先述の記事のような面倒なことはしなくて良くなっています。

そのメソッドが pandas.DataFrame.explode です。

使い方は簡単で、行に展開したい、つまり配列を含んだ列の名前を渡すだけです。

実際にやってみます。


import pandas as pd

df = pd.DataFrame(
    {'A': [[1, 2, 3], 'foo', [], [3, 4]], 'B': 1, "C": [1, 2, 3, 4]})
print(df)
"""
           A  B  C
0  [1, 2, 3]  1  1
1        foo  1  2
2         []  1  3
3     [3, 4]  1  4
"""

df_explode = df.explode('A')
print(df_explode)
"""
     A  B  C
0    1  1  1
0    2  1  1
0    3  1  1
1  foo  1  2
2  NaN  1  3
3    3  1  4
3    4  1  4
"""

めちゃくちゃお手軽ですね。
前の記事で書いていたようにappendの遅さに文句を言ったり、自分で配列を作るコードを書いたりと言ったことは全くしなくて良くなりました。

boto3でEC2インスタンスを起動してIPアドレスを取得する

EC2のインスタンスを操作するとき、(会社ではCLIのバッチファイルを手元に持ってるのですが)私物の環境ではいつもコンソールにログインして起動/停止していたので、
これもバッチ化することにしました。
AWS CLIではなく、 Pythonライブラリの boto3で作ってみます。

職場の環境と異なり、自分の学習環境のEC2はIPアドレスを固定していません。(お金かかるから。)
なので、毎回コンソールから起動して、IPアドレスを確認し、そのIPアドレスに接続する、という作業を行っていました。
これを短縮するため、起動したらIPアドレスを取得して表示できるようにします。

ドキュメントはこの辺りが参考になります。
ec2.html#instance

以下のように、ec2インスタンスのオブジェクトを用意します。


import boto3

ec2 = boto3.resource('ec2')
instance = ec2.Instance('{インスタンスID}')

そして、startで起動してwait_until_running()で起動完了を待ち、
public_ip_addressを取得すれば良いです。

完成したバッチファイルの中身は次のようになります。
これを.pyファイルとして保存して、実行できるように権限を744にします。


#!/usr/bin/env python
import boto3
instance_id = "{インスタンスID}"

ec2 = boto3.resource('ec2')
instance = ec2.Instance(instance_id)
print("インスタンス起動開始")
instance.start()
instance.wait_until_running()
print("インスタンス起動完了")
ip = instance.public_ip_address
print(f"IPアドレス: {ip}")

実行するときちんと、
インスタンス起動開始
インスタンス起動完了
IPアドレス: xx.xxx.xxx.xxx
が表示されました。

ついでですが、停止バッチは次のようになります。


#!/usr/bin/env python
import boto3
instance_id = "{インスタンスID}"

ec2 = boto3.resource('ec2')
instance = ec2.Instance(instance_id)
print("インスタンス停止開始")
instance.stop()
instance.wait_until_stopped()
print("インスタンス停止完了")

Amazon Aurora Serverlessを使ってみる

ずっと前から気になっていたまま、使ったことがなかったのですが重い腰を上げて Aurora Serverless を触ってみることにしました。

用途は仕事ではなく、自分の趣味と勉強でとあるデータを蓄積するDBです。
AWSのアカウント持っていてDBが必要なら、RDSを使うのが順当なのですが、
通常のRDSは個人で使うにはちょっとお高いので、これまではEC2に導入したMySQLを使っていました。
そこにAurora の Serverless が登場ということで、自分の分析作業時のみの課金で使えるとなれば非常にお得な気がします。
完全に移行できるかとコストはトータルどうなるか、といった懸念は残りますが、まずは一度触ってみることにしました。

Aurora Serverless にはいくつか制限事項があるようです。
それを満たせるように事前に準備を進めていきます。

詳細はこちらを参照: Aurora Serverless の制約事項

1. VPC (Virtual Private Cloud)

Aurora Serverless にパブリックなIPアドレスを割り当てることはできず、おなじVPC内からしかアクセスできないそうです。
なので、ローカルの端末から踏み台等使わずに直接繋ぐことはできなさそうですね。
また、EC2インスタンスををクライアントにする場合は、おなじVPCの中に置いておく必要があります。
(僕はVPCはデフォルトの1個しか使っていないので、この制限は特に意識することはありませんでした。
VPCを複数使っている人は注意が必要です。)

2. AWS PrivateLink エンドポイント

僕がこの辺りの用語に疎いので、ドキュメントをそのまま引用します。
各 Aurora Serverless DB クラスターには、2 つの AWS PrivateLink エンドポイントが必要です。VPC 内の AWS PrivateLink エンドポイントが制限に達した場合、その VPC にそれ以上 Aurora Serverless クラスターを作成することはできません。
要は、VPCにサブネットを2つ以上用意しておく必要があります。
これも自分の場合は元々東京リージョンに3個持っていたので特に作業の必要はありませんでした。

3. セキュリティグループ

ドキュメントの通り、3306番のポートを利用するので、「インバウンド」で、3306番ポートが開通したセキュリティグループが必要になります。
セキュリティグループの設定画面で「タイプ」に「MYSQL/Aurora」を選択すると自動的に必要な設定が入るので作っておきます。

4. クライアント

実際の利用にはDBに接続するクライアントが必要です。
MySQLか、その互換のMariaDBクライアントが入ったインスタンスを同じVPC内に立てておきましょう。

さて、準備が終わった後、実際に DBを作成するのですが、これは簡単でした。
手順にすると多いですが画面に沿って進むだけなのでほぼ迷いません。

1. AWSのコンソールにログインし、RDSの管理画面に移動する。
2. [データベースの作成]をクリック
3. 標準作成を選択 (簡単作成でも良い)
4. エンジンのタイプ は [Amazon Aurora]を選択
5. エディション は [MySQLとの互換性を持つ Amazon Aurora]
6. MySQLのバージョンを選ぶ。
  これは、ドキュメントで指定されているバージョンでないと、次の選択肢からサーバレスが消えるので注意が必要です。

7. データベースの機能 は [サーバーレス]
8. DBクラスター識別子 は何か名前をつける。
9. マスターユーザー名 を指定 (デフォルトは admin)
10. パスワードの設定。 (自分はパスワードの自動作成にしました)
11. キャパシティーの設定を変更。(予算を抑えるため気持ち低めにしました。)
12. [データベースの作成]をクリック

ここまで進むとDBの作成が始まり、数分でDBが出来上がります。
出来上がったら、[認証情報の詳細を表示] から、 adminのパスワードを入手します。
また、同時に エンドポイント もわかるのでこれも記録しておきます。

DBが作成できたら、クライアントの環境から以下のコマンドで接続できます。


mysql -h {エンドポイント} -u {ユーザー名} -p

Enter password: 
と聞かれるのでパスワードを入力する。

ポートはデフォルトの3306を使っているので、 -P 3306 はつける必要ありません。(つけても大丈夫ですが。)

これで、通常のDBとして、使用できるようになりました。

もしDBが作成できているのに接続できなかった場合、セキュリティグループを見直してみてください。
(僕は用意していたセキュリティグループと違うグループが付与されていてしばらく詰まりました。)

あとは、残る課題はお値段ですね。
しばらく使ってみて、どのくらいの金額になるのか様子をみようと思います。

※ 以下、追記
2~3日ほど放置した後、様子を見たらメキメキと課金されていました。触っていない間も稼働しっぱなしだったようです。

設定項目の中に、「数分間アイドル状態のままの場合コンピューティング性能を一時停止する」
というがあり、これにチェックを入れないと、期待してた使っていない間は停止する効果が得られないようです。
てっきりデフォルトだと思っていたので危なかったです。
(請求金額アラートに救われました)

matplotlibでレーダーチャート(メモリも多角形)を描写する

Pythonでレーダーチャートを書きたくなり、matplotlibでやってみたのでそのメモです。
公式にサンプルがあるのですが、2020年09月 現在うまく動きません。
参考: api example code: radar_chart.py
実装自体も、PolarAxesクラスを継承してメソッドを書き換えるかなり仰々しいものですし、
メモリが円形のままで、僕が望む形ではなかったのでゼロベースでやってみました。

まずライブラリをインポートして、適当にデータを作っておきます。
データはレーダーチャートで可視化する値とそれぞれのラベルがあれば良いです。


import matplotlib.pyplot as plt
import numpy as np

values = np.array([31, 18, 96, 53, 68])
labels = [f"データ{i}" for i in range(1, len(values)+1)]

さて、レーダーチャートですが、見栄えにこだわりがなければ簡単に書くことができます。
matplotlibで極座標のグラフを作り、一周ぐるっとplotするだけです。
簡易版ですがそのコードを先に紹介します。


# 多角形を閉じるためにデータの最後に最初の値を追加する。
radar_values = np.concatenate([values, [values[0]]])
# プロットする角度を生成する。
angles = np.linspace(0, 2 * np.pi, len(labels) + 1, endpoint=True)

fig = plt.figure(facecolor="w")
# 極座標でaxを作成。
ax = fig.add_subplot(1, 1, 1, polar=True)
# レーダーチャートの線を引く
ax.plot(angles, radar_values)
# レーダーチャートの内側を塗りつぶす
ax.fill(angles, radar_values, alpha=0.2)
# 項目ラベルの表示
ax.set_thetagrids(angles[:-1] * 180 / np.pi, labels)

ax.set_title("レーダーチャート", pad=20)
plt.show()

出力がこちらです。

これでも最低限の要件は満たせますね。ただ、データを表してる青の線が真っ直ぐなのに、
その目安となるメモリ線が円形なのが気になります。
また、普通の曲座標と違って、ラベルを上を始点にして時計回りにしたいです。

このラベルの開始位置と回転方向を変えるのは簡単なのですが、メモリを多角形にするのはまともに取り組むと非常に大変です。
(公式サンプルの様なClassを継承しての大がかりな改修が必要になります。)

なので、アプローチを変えてみました。
matplotlibの機能で引いてくれるメモリ線は全部消してしまいます。
そして、定数値のレーダーチャートとして、灰色の線を自分で引きました。

出来上がったコードがこちらです。


# 多角形を閉じるためにデータの最後に最初の値を追加する。
radar_values = np.concatenate([values, [values[0]]])
# プロットする角度を生成する。
angles = np.linspace(0, 2 * np.pi, len(labels) + 1, endpoint=True)
# メモリ軸の生成
rgrids = [0, 20, 40, 60, 80, 100]


fig = plt.figure(facecolor="w")
# 極座標でaxを作成
ax = fig.add_subplot(1, 1, 1, polar=True)
# レーダーチャートの線を引く
ax.plot(angles, radar_values)
# レーダーチャートの内側を塗りつぶす
ax.fill(angles, radar_values, alpha=0.2)
# 項目ラベルの表示
ax.set_thetagrids(angles[:-1] * 180 / np.pi, labels)
# 円形の目盛線を消す
ax.set_rgrids([])
# 一番外側の円を消す
ax.spines['polar'].set_visible(False)
# 始点を上(北)に変更
ax.set_theta_zero_location("N")
# 時計回りに変更(デフォルトの逆回り)
ax.set_theta_direction(-1)

# 多角形の目盛線を引く
for grid_value in rgrids:
    grid_values = [grid_value] * (len(labels)+1)
    ax.plot(angles, grid_values, color="gray",  linewidth=0.5)

# メモリの値を表示する
for t in rgrids:
    # xが偏角、yが絶対値でテキストの表示場所が指定される
    ax.text(x=0, y=t, s=t)
    
# rの範囲を指定
ax.set_rlim([min(rgrids), max(rgrids)])

ax.set_title("レーダーチャート", pad=20)
plt.show()

出力がこちら。

自分がイメージしていたのに近いものが作れました。
コード中にコメントを全部入れたので、ここからさらに見た目を変える場合はすぐ改良できると思います。

バートレット検定

前回の記事の中で、SciPyにはF検定の関数が実装されていないという話を書きました。
参考: PythonでF検定を実装する

しかし、その代わりにSciPyには分散が等しいことを検定する別の方法が実装されています。
その一つが、バートレット検定です。
関数はこれです。
scipy.stats.bartlett

複数(2個だけでなく3個以上も可能)のサンプルを渡してあげると、統計量とp値を計算してくれます。

今回の記事では、使い方の前にバートレット検定の統計量を紹介します。
参考にしたのは、東京大学出版会の自然科学の統計学 (青い本) の 94ページです。

サンプルの数を$a$とし、それぞれのサンプルのサイズを$n_i \ \ (j=1,\dots, a)$とします。
また、サンプルサイズの合計を$n=\sum_{j}n_j$と置いておきます。
各サンプルの不偏分散を
$$
V_i = \sum_{j} \frac{(y_{ij}-\bar{y_i})^2}{n_i-1}
$$
として、さらにこれらを併合したものを
$$
V_e = \sum_i \frac{(n_i-1)V_i}{n-a} = \sum_i \sum_j \frac{(y_{ij}-\bar{y_i})^2}{n-a}
$$
とします。
さらに、
$$
B = (n-a)\log{V_e} -\sum_i (n_i-1)\log{V_i}
$$
と置いて、
$$
B’ = \frac{B}{
1+\frac{1}{3(a-1)}\{\sum_i \frac{1}{n_i-1} – \frac{1}{n-a}\}
}
$$
と補正すると、この$B’$は自由度$a-1$の$\chi^2$分布に近似的にしたがいます。
それを利用して、検定を行うのがバートレット検定です。

数式に沿ってNumPyで実装してみたのが次の関数です。ランダム生成した3サンプルで試してみました。


from scipy.stats import norm
from scipy.stats import chi2
import numpy as np

# スクラッチで実装したバーレット検定量の算出関数
def bartlett_value(*args):
    # サンプル数
    a = len(args)
    # 各サンプルの サンプルサイズと不偏分散
    ni = np.zeros(a)
    Vi = np.zeros(a)
    for j in range(a):
        ni[j] = len(args[j])
        Vi[j] = np.var(args[j], ddof=1)

    # サンプルサイズの合計
    N = np.sum(ni)
    # Veの計算
    Ve = np.sum((ni-1)*Vi) / (N-a)

    # 統計量の計算
    B = (N-a) * np.log(Ve) - np.sum((ni-1)*np.log(Vi))
    C = 1+1/(3*(a-1))*(np.sum(1/(ni-1))-1/(N-a))

    return B/C


# サンプルを3種類生成する
data_a = norm(loc=3, scale=5).rvs(20)
data_b = norm(loc=7, scale=8).rvs(30)
data_c = norm(loc=-5, scale=5).rvs(15)

# 統計量を算出
b_value = bartlett_value(data_a, data_b, data_c)

# p値を算出
p_value = chi2(3-1).sf(b_value)

print("B値:", b_value)
# B値: 8.26319721031961
print("p値:", p_value)
# p値: 0.016057189188779606

$\alpha=0.05$とすると帰無仮説が棄却され、分散が等しくないことが検出されていますね。

もちろん、SciPyに関数は用意されているので、普段はこの様にスクラッチで実装する必要はなく、
呼び出すだけで良いです。結果が一致することを見ておきます。


from scipy.stats import bartlett


b_value, p_value = bartlett(data_a, data_b, data_c)
print("B値:", b_value)
# B値: 8.26319721031961
print("p値:", p_value)
# p値: 0.016057189188779606

バッチリ一致しました。

PythonでF検定を実装する

最近、とある二つのサンプルの分散が異なることを確認する機会があり、F検定を行う必要がありました。
いつもみたいにSciPyですぐにできるだろうと考えていたのですが、
Statistical functions をみる限りでは、SciPyにF検定は実装されていなさそうでした。

その代わり、バートレット検定とルビーン検定が実装されているのですが。
この二つの検定については別途勉強するとして、とりあえずF検定を行いたかったので自分で実装しました。
幸い、F分布自体はSciPyにあるので簡単です。
t検定のメソッドを参考にし、サンプルを二つ渡したらF統計量とp値を返してくれる関数ように実装しました。

※F検定自体の説明は今回は省略します。
興味のあるかたには、東大出版会の統計学入門(いわゆる赤本)の244ページ、12.2.4 母分散の比の検定 あたりが参考になります。


from scipy.stats import f
import numpy as np


def ftest(a, b):
    # 統計量Fの計算
    v1 = np.var(a, ddof=1)
    v2 = np.var(b, ddof=1)
    n1 = len(a)
    n2 = len(b)
    f_value = v1/v2

    # 帰無仮説が正しい場合にFが従う確率分を生成
    f_frozen = f.freeze(dfn=n1-1, dfd=n2-1)

    # 右側
    p1 = f_frozen.sf(f_value)
    # 左側
    p2 = f_frozen.cdf(f_value)
    # 小さい方の2倍がp値
    p_value = min(p1, p2) * 2

    # 統計量Fとp値を返す
    return f_value, p_value

きちんと実装できているかどうか確認するため、統計学入門の練習問題を一つ解いておきます。
252ページの練習問題、 12.2 の iii) が分散の検定なのでやってみます。

次のコードのdata_1とdata_2の分散が等しいというのが帰無仮説で、等しくないというのが対立仮説です。


# 問題文のデータを入力
data_1 = np.array([15.4, 18.3, 16.5, 17.4, 18.9, 17.2, 15.0, 15.7, 17.9, 16.5])
data_2 = np.array([14.2, 15.9, 16.0, 14.0, 17.0, 13.8, 15.2, 14.5, 15.0, 14.4])

# F統計量とp値を計算
f_value, p_value = ftest(data_1, data_2)
print(f"F統計量: {f_value:1.3f}")
# F統計量: 1.564
print(f"p値: {p_value:1.3f}")
# p値: 0.516

F統計量は正解と一致しましたし、帰無仮説が棄却できないのも確認できました。

GCPのテキストの感情分析をPython SDKでやってみる

前回の記事ではAPIを直接叩きましたが、現実的にはSDKを使った方が手軽なのでその方法を紹介します。

準備として、 クイックスタート: クライアント ライブラリの使用 をみながら以下の内容を行っておきます。
1. サービスアカウントの作成。
2. jsonファイルのローカル保存。
3. jsonファイルのパスを環境変数 GOOGLE_APPLICATION_CREDENTIALS に設定。
4. SDK google-cloud-language のインストール

サービスアカウントを使わずに既存のAPIキーを使ってどうにかできないか方法を探してたのですがそれは無理そうですね。
諦めてサービスアカウントを作りましょう。

準備ができたら実行してみます。 ハローワールド的なコード例があるのでほとんどそのまま動かしました。


# ライブラリのインポート
from google.cloud import language
from google.cloud.language import enums
from google.cloud.language import types

# クライアントのインスタン作成
client = language.LanguageServiceClient()

# 分析対象のテキスト
text = u'Hello, world!'

document = types.Document(
    content=text,
    type=enums.Document.Type.PLAIN_TEXT
)

# テキストの感情分析
sentiment = client.analyze_sentiment(document=document).document_sentiment

print(f"Text: {text}")
print(f"Sentiment: {sentiment.score}, {sentiment.magnitude}")

"""
Text: Hello, world!
Sentiment: 0.6000000238418579, 0.6000000238418579
"""

お作法に慣れるまで少しだけかかりそうなのですが、それでも簡単に動かせましたね。
Hello, world! はポジティブな言葉のようです。

GCPのCloud Natural Language API でセンチメント分析

せっかくアカウントを作ったのでGCPでテキストのセンチメント分析(ポジネガ)を試してみます。
Python ライブラリも用意されている様なのですが、最初なのでAPIを直接叩いてみました。

こちらに、 curl コマンドを使う例があり、これを参考に、pythonで書きなおしてみます。


curl -X POST \
     -H "Authorization: Bearer "$(gcloud auth application-default print-access-token) \
     -H "Content-Type: application/json; charset=utf-8" \
     --data "{
  'encodingType': 'UTF8',
  'document': {
    'type': 'PLAIN_TEXT',
    'content': 'Enjoy your vacation!'
  }
}" "https://language.googleapis.com/v1/documents:analyzeSentiment"

今回は認証はAPIキーで行うので、URLにクエリパラメーター(key={APIキー})で流けることに注意しながら、pythonでやってみました。


# サンプルテキスト
text = (
    "メロスは激怒した。必ず、かの邪智暴虐の王を除かなければならぬと決意した。"
    "メロスには政治がわからぬ。メロスは、村の牧人である。笛を吹き、羊と遊んで暮して来た。"
    "けれども邪悪に対しては、人一倍に敏感であった。"
    "きょう未明メロスは村を出発し、野を越え山越え、十里はなれた此のシラクスの市にやって来た。"
)

# 自分のAPIキーに置き換える
apl_key = "{自分のAPIキー}"

# APIのurl情報
url = 'https://language.googleapis.com/v1/documents:analyzeSentiment?key=' + apl_key

# ヘッダーとデータの設定
header = {'Content-Type': 'application/json'}
data = {
    "document": {
        "type": "PLAIN_TEXT",
        "content": text
    },
    "encodingType": "UTF8"
}

response = requests.post(url, headers=header, json=data).json()

# 結果表示
response
"""
{'documentSentiment': {'magnitude': 1.2, 'score': 0},
 'language': 'ja',
 'sentences': [{'text': {'content': 'メロスは激怒した。', 'beginOffset': 0},
   'sentiment': {'magnitude': 0.1, 'score': -0.1}},
  {'text': {'content': '必ず、かの邪智暴虐の王を除かなければならぬと決意した。', 'beginOffset': 27},
   'sentiment': {'magnitude': 0.3, 'score': 0.3}},
  {'text': {'content': 'メロスには政治がわからぬ。', 'beginOffset': 108},
   'sentiment': {'magnitude': 0.5, 'score': -0.5}},
  {'text': {'content': 'メロスは、村の牧人である。', 'beginOffset': 147},
   'sentiment': {'magnitude': 0, 'score': 0}},
  {'text': {'content': '笛を吹き、羊と遊んで暮して来た。', 'beginOffset': 186},
   'sentiment': {'magnitude': 0, 'score': 0}},
  {'text': {'content': 'けれども邪悪に対しては、人一倍に敏感であった。', 'beginOffset': 234},
   'sentiment': {'magnitude': 0, 'score': 0}},
  {'text': {'content': 'きょう未明メロスは村を出発し、野を越え山越え、十里はなれた此のシラクスの市にやって来た。',
    'beginOffset': 303},
   'sentiment': {'magnitude': 0, 'score': 0}}]}
"""

ドキュメント全体としての、magnitudeとscoreに加えて、文ごとにもそれぞれ表示してくれているのがありがたいですね。
(わかりにくいですが、textは全部連結されて一つの文字列になっているので、GCPが自動的に分けてくれたものになります。)

‘メロスは激怒した。’ の score -0.1 より、
‘メロスには政治がわからぬ。’ の score -0.5 の方がネガティブと判定されているのですね。

結果の解釈については、
感情分析の値の解釈
という文書があるのでこちらを読んでおきましょう。
スコアが0.0付近のものの解釈などかなり重要な情報も含まれます。