投資のためのデータサイエンス

個人の投資活動に役立つデータ分析にまつわる話題を綴ります。

Pythonによる分散分析のデモ

ある変量の変動の要因を分析する場合には、母集団のサブグループ別に分布や平均値を比較して、グループ間で差があるかどうかをみる、ということはよく行われます。これをより系統的に行う手法が分散分析です。分散分析は、原理的にはダミー変数を用いた回帰分析と同じですが、例えば「性別」や「年齢階級」といった要因単位で指標値に有意に影響を及ぼしているかを見るような場合に便利です。ここでは都道府県別の各種指標値を例題として、「(地理的に隣接した)地方」や「大都市があるか」といった要因で異なるかどうかを見ることにします。
まず、必要なライブラリをインポートします。

# ライブラリのインポート
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font=['IPAexGothic'])
import japanize_matplotlib # matplotlib日本語
from google.colab import files # 画像ファイルダウンロード用

ここでは、日本の47都道府県をクロスセクションデータとして用いることにします。まず、エクセルファイルからPandasデータフレームに都道府県の属性データを読込みます。

# ファイルからデータを読み込む
df01 = pd.read_excel('/content/drive/My Drive/Colab Notebooks/データ解析講座一日目演習.xlsx',sheet_name='都道府県')
df01.head()

f:id:nicjps230:20210902182435p:plain
上記のデータ列のうち、「都市化」は政令指定都市がある道府県から北海道を除いて東京都と香川県(四国経済の拠点が多く存在)を加えたものを「都市県」、その他を「郊外県」としました。また、「メディア」は以下の記事にあるメディア接触の特徴に関する分類です。 https://news.livedoor.com/article/detail/9064551/
同様に、公的ウェブサイトに掲載されている、47都道府県の家計指標データを読込みます。

# ファイルからデータを読み込む
df02 = pd.read_excel('/content/drive/My Drive/Colab Notebooks/データ解析講座一日目演習.xlsx',sheet_name='家計指標')
df02.head()

f:id:nicjps230:20210902182854p:plain
そして、読み込んだ2つのデータフレームを、県名をキーとして横に連結します。

# 県名をキーにして横連結する
df1 = pd.merge(df01, df02, left_on='県名', right_on='県名')
df1.drop(columns=['コード_x','コード_y'],inplace=True)

一元配置分散分析
ここでは都道府県の食料費割合の変動に及ぼす「地域」等の属性要因の影響をみるために、一元配置分散分析を実行します。statsmodelsライブラリを用います。

# 一元配置分散分析
import statsmodels.formula.api as smf
import statsmodels.api as sm
oneway_anova = smf.ols(formula='食料費割合~地域+都市化+メディア',data=df1).fit()
sm.stats.anova_lm(oneway_anova,typ=2).round(2)

これにより、以下のように各要因についてのF検定統計量が出力されます。
f:id:nicjps230:20210902183641p:plain
尚、交互作用項を加える場合はモデル式に"地域*都市化"のように*でつないだ項を追加します。
有意な主効果の様子を見る
ここではseabornライブラリのcatplot()を用いて、各水準毎の指標値の平均を棒グラフで比較します。

# 有意差のある効果についてその中身をみる
import matplotlib.pyplot as plt
import seaborn as sns
p = sns.catplot(x='食料費割合',y='都市化',data=df1,kind='bar',height=4,aspect=2)
p.set(title='食料費割合の都市化による差')

以下のような形で平均の違いを比較することができます。
f:id:nicjps230:20210902184340p:plain
交互作用の様子を見る
交互作用(ある要因による指標への影響が、他の要因の水準によって異なるかどうか)の様子を見るためには、statsmodelsライブラリのinteraction_plotを用います。

from statsmodels.graphics.factorplots import interaction_plot
fig=interaction_plot(df1.地域,df1.都市化,df1.食料費割合,colors=['red','blue'],ms=10)
plt.title('地域別食料費割合と都市化')
plt.legend(bbox_to_anchor=(1.05,1.0),loc='upper left')

これにより交互作用の様子がひと目でわかるようになります。近畿地方において、都市部と郊外部で食料費割合の差が大きいことがわかります。
f:id:nicjps230:20210902184835p:plain

データサイエンスのおすすめオンライン記事(8月9日付)

前回のこのシリーズの投稿からバタバタしている間にオリンピックが終わってしまった。これは最近始まったことではないが、バレーボールではコーチがタブレット端末片手にデータ分析結果に基づいてリアルタイムに戦略を指南していたし、サッカーでも日本代表チームには優れた分析スタッフがいて、自チームや対戦する相手チームの直近の試合のデータを徹底的に分析しているという話をしていた。スポーツのデータ活用は予算と人材さえあればどこでも始められるようになってきている。
今回は「機械学習」と「スクレイピング」の記事を中心に紹介する。

Pythonコードを高速化する方法
www.kdnuggets.com
インタープリター言語であるPythonを少しでも早く走らせるプログラミング上のコツをいくつか挙げている。

ウェブページのスクレイピングに回帰ベースの機械学習モデルを用いる
www.datasciencecentral.com
ユーザエージェントとブロックされる可能性の関連をモデル化してスクレイピングの効率化を図る。

225個のPythonによる機械学習プロジェクト
medium.com
Python言語を用いた機械学習プロジェクト(問題・データ・コードが掲載されており自分で試して学べる)が225個も紹介さている良記事。初級者向けプロジェクトは最初の方でまとめて紹介されている。

Amazonの商品カタログデータをスクレイピングする方法
www.datasciencecentral.com
Pythonのscrapyライブラリを用いてAmazonの商品カタログデータをスクレイピングする具体的な方法の解説。

顔認識:最新の方法と最良のツール
www.datasciencecentral.com
「知識ベース」「テンプレート照合」「特徴量ベース」「外見ベース」といった顔認識の最新技術の解説と有料/無料で利用できる顔認識ツールの紹介。

PythonによるECカタログデータの分析小技集(下)

前回に引き続き、「楽天ブックス書籍検索API」を利用してダウンロードした書籍情報データを用いて、ECサイトのカタログデータを処理・分析するためのpythonの小技をまとめました。
まず前回と同様に必要なライブラリを読み込み、データを保存先からロードします。

# ライブラリのインポート
import numpy as np
import pandas as pd
import re
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font=['IPAexGothic'])
# parquetファイルからロード
df01 = pd.read_parquet('rakuten_books.parquet')
# 列を絞り込む
focus_cols=['title','subTitle','seriesName','contents','author','publisherName','size',\
	'itemCaption','salesDate','itemPrice','reviewCount','reviewAverage']
df1 = df01.loc[:,focus_cols]

定型テキストの出現頻度を求める
テキスト情報の列の中でも、出現する文字列のパターンが限定されている場合があります。いわゆる「カテゴリカルデータ」に相当するものです。その出現頻度を知りたいことはよくあります。これにはcollectionsライブラリのCounter()メソッドを用います。ここでは列'publisherName'(出版社)の頻度分布を求め、上位10社を出力します。

# 頻度分布を求める列名を指定する
col_name = 'publisherName'
import collections
c_count = collections.Counter(df1[col_name])
# カウント上位のリストを得る
p_names, p_counts = zip(*c_count.most_common(10))
print(p_names)

尚、出力結果は以下のようになりました。
('講談社', 'KADOKAWA', '集英社', '小学館', '宝島社', 'スクウェア・エニックス', '白泉社', '秋田書店', '幻冬舎', '朝日新聞出版')

正規表現による検索
決まった文字列ではなく、ある規則性を持った文字列を検索したい場合には「正規表現」による検索を用います。ここでは、コミックセットのように複数の書籍が一括して販売されているケースを検索し、「○冊セット」の○にあたる数字を取り出すことを考えます。ここではpythonライブラリ're'(regular expressionを略したもの)のsearchメソッド(と内容を取り出すためのgroupメソッド)を一件につき2回用いています。1回目の適用で'○冊セット'の文字列をとりだし、2回目の適用でその中の○にあたる数字の部分を取り出して、int()で数値化しています。

def func_bulk(s):
    """タイトル文字列sの中の冊数を表すパターンを検索して冊数の数値を返す"""
    nitem = 0
    blkstr1 = re.search('[0-9]+巻セット',s)
    if not blkstr1 is None:
        nitem = int(re.search('[0-9]+',blkstr1.group()).group())
        return nitem
    blkstr2 = re.search('[0-9]+冊セット',s)
    if not blkstr2 is None:
        nitem = int(re.search('[0-9]+',blkstr2.group()).group())
        return nitem
    blkstr3 = re.search('全[0-9]+巻',s)
    if not blkstr3 is None:
        nitem = int(re.search('[0-9]+',blkstr3.group()).group())
        return nitem
    blkstr4 = re.search('全[0-9]+冊',s)
    if not blkstr4 is None:
        nitem = int(re.search('[0-9]+',blkstr4.group()).group())
        return nitem
    return nitem
# セット書籍の冊数を取り出す
df1 = df1.assign(n_bulk = df1['title'].apply(func_bulk))
df81 = df1[df1['n_bulk'] > 0]
dfi.export(df81,'kensaku_bulk.png')

結果は以下のようになりました。
f:id:nicjps230:20210723155132p:plain

データフレーム内の複数の列の演算により新たな列を作成する
これまでにもデータフレーム内のある一つの列から、ある演算によって新しい列を作り出す例が出てきました。演算が簡単なものならばその列を作り出す文の中に式を埋め込めばいいのですが、式が複雑な場合は別途関数を定義して、元の列に関してapplyメソッドを使う、という形でした。それでは、データフレーム内の複数の列を使って新しい列を作り出す場合はどうするのでしょうか。その場合は、データフレームそのものについてappyメソッドを適用します。ここでは列'reviewCount'と列'reviewAverage'を用いて、新しい列'review_pi'を作成しています。

# 新しい列を作成するための関数
def func_review_pi(df):
  return np.log(df.loc['reviewCount']+1) + df.loc['reviewAverage']
# データフレーム内の複数の列の演算により新たな列を作成する
df51 = df1.assign(review_pi = df1.apply(func_review_pi, axis=1))

新しい列を作成するための関数に他の引数がある場合
これまでの例では、データフレーム内の新たな列を作るための関数には「元の列」や「元のデータフレーム」以外の引数はありませんでした。今度は、その関数が別途他の引数を必要とする場合の記述を見てみます。 ここでの例は、書籍のタイトルにある指定したキーワードが含まれていたら、あるジャンルの書籍であると判定するための関数の作成です。一つのジャンルについてキーワードは複数あるものとし、キーワードは2次元のリストで与えられるとします。

# ジャンル別のキーワードのリストの定義
# 0:経済, 1:歴史, 2:科学, 3:医学, 4:旅行, 5:辞典, 6:人生, 7:食事, 8:情報, 9:入試, 10:漫画
kws = [['経済','財務','お金','家計','会計','金融'],
       ['史','古代'],
       ['科学','サイエンス','宇宙'],
       ['医','癌','内科','外科','歯科','手術','疾患'],
       ['旅','紀行'],
       ['辞典','広辞苑'],
       ['人生','生き方','生きかた','ライフ'],
       ['料理','食事','グルメ','レシピ'],
       ['Excel','Word','Python','iPhone','iPad','プログラミング','アプリ','HTML','ソフトウェア','情報技術','Adobe'],
       ['入試','中学受験','高校受験','大学受験'],
       ['鬼滅の刃','リベンジャーズ','るろうに剣心','スライムだった件','青の祓魔師','呪術廻戦','進撃の巨人']]
# キーワード判定関数の定義
def func_genre(s, ww, k):
  """文字列sに2次元リストwwの第kジャンルのキーワードが
  含まれていれば1, 含まれていなければ0を返す"""
  q = 0
  for j in range(len(ww[k])):
    if ww[k][j] in s:
      q = 1
  return q

演算の元になる列以外の引数がある場合、その引数はapply()メソッドの中で引数名=値の形で指定します。

# タイトルに経済分野のキーワードを含むかどうかのダミー変数を生成する
df1= df1.assign(G経済=df1['title'].apply(func_genre, ww=kws, k=0))

PythonによるECカタログデータの分析小技集(上)

一部のECサイトでは、販売商品のカタログデータをAPIにより取得することができます。このようなデータには、商品名、商品の補足説明、価格などの情報が含まれています。ここでは、このようなECサイトのカタログデータを処理・分析するためのpythonの小技をまとめました。ここで用いるデータは「楽天ブックス書籍検索API」を利用してダウンロードした書籍情報データです。
まず必要なライブラリを読み込み、データを保存先からロードします。

# ライブラリのインポート
import numpy as np
import pandas as pd
import re
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font=['IPAexGothic'])
# parquetファイルからロード
df01 = pd.read_parquet('rakuten_books.parquet')
# 列を絞り込む
focus_cols=['title','subTitle','seriesName','contents','author','publisherName','size',\
	'itemCaption','salesDate','itemPrice','reviewCount','reviewAverage']
df1 = df01.loc[:,focus_cols]

テキスト項目の文字列検索
まず、テキスト情報からあるキーワードを含む行を抽出するケースを考えます。pandasのSeries.str.contains()を用います。これにより、指定した文字列を含む行のみを取り出すことができます。またこのコードでは、Jupyter Notebookの出力イメージをそのまま画像ファイルに変換するdataframe_imageライブラリを用いています。

import dataframe_image as dfi
# テキスト項目で特定の文字列を検索する
col_name = 'title'
ch_kw = '資産'
df201 = df1[df1[col_name].str.contains(ch_kw)]
dfi.export(df201,'kensaku_sisan.png')

以下のような検索結果が得られます。
f:id:nicjps230:20210709161416p:plain
次に、2つのキーワードについて、その両方を含む行を取り出したい(and条件)場合は、一つのキーワードで絞りこんだデータフレームに対してもう一つのキーワードでさらに絞りこみをかけることでできます。【注】一つの命令文にブール演算を埋め込んで検索することは、下に記載したquery()メソッドではできますが、それ以外ではエラーが出るためここでは避けています。

# テキスト項目で特定の複数の文字列をand検索する
col_name = 'itemCaption'
ch_kw1 = '成功'; ch_kw2 = '金'
df202 = df1[df1[col_name].str.contains(ch_kw1)]
df203 = df202[df202[col_name].str.contains(ch_kw2)]
print('該当するデータ件数 = ',len(df203))

一方、2つのキーワードのどちらかが含まれる行を抽出したい(or条件)場合は、元のデータフレームを個別のキーワードで絞りこんだ結果を縦結合することでできます。これは検索対象の列が異なる場合にも応用できます。どちらかの検索結果がゼロ行でも特に問題はありません。

# テキスト項目で特定の複数の文字列をor検索する
col_name = 'title'
ch_kw1 = '資産'; ch_kw2 = '財産'
df204 = df1[df1[col_name].str.contains(ch_kw1)]
df205 = df1[df1[col_name].str.contains(ch_kw2)]
df206 = pd.concat([df204,df205])
print('該当するデータ件数 = ',len(df206))

特定の数値項目で条件検索する
今度は数値データが入った列について、ある条件を満たす行を取り出すことを考えます。例えば値がゼロより大きい行を抽出するには、df[df['列名']>0]とする書き方と、df.query('列名>0')という書き方があります。

# 数値変数で条件検索する
ch_col='reviewCount'; s_val=500
df31 = df1[df1[ch_col] >= s_val]
dfi.export(df31,'top_review_no.png')
#数値検索の別の記述法
ch_col='reviewCount'; s_val=500
df32 = df1.query('{0} >= @s_val'.format(ch_col))
dfi.export(df32,'top_review_no2.png')

結果は以下のようになります。
f:id:nicjps230:20210719000425p:plain
特定の列項目が空でないケースを抽出する
次に、テキストからなる特定の列項目が空でないケースを抽出する方法を考えます。テキスト情報で何も書かれていないのは空文字列('')なので、まず空文字列を欠損値(NaN)に置き換えます。その後で欠損値を含む行を削除します。

# 列名を指定する
col_name = 'contents'
# データフレームをコピーする
df211 = df1.copy()
#空白をまずNaNに置き換え
df211[col_name].replace('', np.nan, inplace=True)
#Nanを削除 inplace=Trueでdfが上書きされる。
df211.dropna(subset=[col_name], inplace=True)
print('該当するデータ件数 = ',len(df211))

(下巻へ続く)

機械学習モデルの数学的記述と計算(その1):バイナリ学習

機械学習は、データからそれらが生成されている法則性を学んで、予測や分類に役立てるものである。
ここでは、スタンフォード大学のAndrew Ng教授の講義ノートをベースに、機械学習の数学的バックグラウンドを整理する。

教師ありバイナリ学習のデータの表現

教師ありバイナリ学習に用いられるデータを、特徴量とターゲットのペアとして考える。これを、(x,y) により表現する。ここで、 x \in \mathbb{R}^{n_{x}} は特徴量のベクトルを表し、 y \in \{0,1\} はターゲット変数を表す。
m個の学習用データの組は以下で表される


(x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\cdots,(x^{(m)},y^{(m)})

特徴量とターゲット変数をそれぞれまとめて以下のように記述する


x=(x^{(1)},x^{(2)},\cdots,x^{(m)}),\ \ x \in \mathbb{R}^{n_{x} \times m}

y=(y^{(1)},y^{(2)},\cdots,y^{(m)}),\ \ y \in \mathbb{R}^{1 \times m}

ロジスティック回帰モデル
x \in \mathbb{R}^{n_{x}} が与えられたとき、 以下の式により \hat{y}=P(y=1|X),0 \le \hat{y} \le 1を得たいとする。
(パラメータとアウトプット)


w \in \mathbb{R}^{n_{x}},\ \ \ b \in \mathbb{R}

\hat{y}=\sigma(w ^{T} x+b)
ここで、

\sigma (z)=\frac{1}{1+\exp (-z)}

ロジスティック回帰モデルの損失関数
上記の定式化で、\{(x^{(1)},y^{(1)}),\cdots,(x^{(m)},y^{(m)})\} が与えられたとき、 \hat{y}^{(i)} \approx y^{(i)} を得たいとする。
損失(誤差)関数を以下のような「交差エントロピー誤差関数」として定義する。


J(w, b)=\frac{1}{m}\sum_{i=1}^{m}\mathcal{L}(\hat{y},y)=-\frac{1}{m}\sum_{i=1}^{m}(y^{(i)}\log \hat{y}^{(i)}+(1-y^{(i)})\log (1-\hat{y}^{(i)}))

最急降下法
上の式の J(w,b)を最小化する w,bを見つけたいとする。例えば J wの関数としてみた場合、 Jの勾配を計算して、勾配が急な方向に沿って進む。つまり以下を反復計算する


w:=w-\alpha \frac{dJ(w)}{dw}
ここで \alphaは学習率(learning rate)である。

計算グラフ
機械学習では以下のような計算グラフを考える。そして、その微分を逆伝搬で考える。これにより計算を局所化し、計算効率が高まる。
ここでは、 J(a,b,c)=3(a+bc) を例として考える。
f:id:nicjps230:20210711130639p:plain

上のロジスティック回帰モデルの計算グラフとその微分(概要のみ)は以下のようになる。
f:id:nicjps230:20210716140912p:plain

m個のデータ例でのロジスティック回帰モデル
損失関数の式を再掲する


J(w, b)=\frac{1}{m}\sum_{i=1}^{m}\mathcal{L}(a^{(i)},y^{(i)})
ここで、

a^{(i)}=\hat{y}^{(i)}=\sigma (w^{T}x^{(i)}+b)
いま、パラメータが w_{1},w_{2},bのみであるとする。上記の損失関数の、パラメータ w_{1}に関する傾きを求めるためには、以下を計算する。

\frac{\partial }{\partial w_{1}}J(w,b)=\frac{1}{m}\sum_{i=1}^{m}\frac{\partial }{\partial w_{1}}\mathcal{L}(a^{(i)},y^{(i)})

上記のように微分を逆伝搬で考える場合、各ノードの変数が1単位変化した時に最終出力変数がどのくらい変化するか、すなわち

\frac{\partial (final output var)}{\partial (var)}=dvar
を考えることになる。例えば、 w_{1}についての勾配は、上の図の微分の逆伝搬を伝っていき、

dw_{1}=\frac{\partial \mathcal{L}}{\partial w_{1}}=\frac{\partial \mathcal{L}}{\partial a}\cdot \frac{\partial a}{\partial z}\cdot x_{1}
と求められる。いま、

\mathcal{L}(a,y)=-(y\log a+(1-y)\log (1-a))
 aについての「勾配」を求めると、

da=\frac{\partial \mathcal{L}(a,y)}{\partial a}=-\frac{y}{a}+\frac{1-y}{1-a}
となる。
また、

a=\sigma (z)=\frac{1}{1+\exp (-z)}
 zについて解くと、

z=-\log (\frac{1}{a}-1)
となり、逆関数微分法より、

\frac{\partial a}{\partial z}=\frac{1}{\frac{\partial z}{\partial a}}=a(1-a)
となる。それゆえ、

dz=\frac{\partial \mathcal{L}}{\partial z}=\frac{\partial \mathcal{L}}{\partial a} \cdot \frac{\partial a}{\partial z}=(-\frac{y}{a}+\frac{1-y}{1-a}) \cdot a(1-a)=a-y
となる。
さらに、

z=w_{1}x_{1}+w_{2}x_{2}+b
より、

dw_{1}=\frac{\partial \mathcal{L}}{\partial w_{1}}=\frac{\partial \mathcal{L}}{\partial z} \cdot \frac{\partial z}{\partial w_{1}}=dz \cdot x_{1}
となる。同様に、 dw_{2}=x_{2} \cdot dz,   db=dzとなる。

アルゴリズムpythonによる実装
上記のアルゴリズムでロジスティック回帰の係数を求める計算をpythonで実装した。データはKaggle上にあるロジスティック回帰用のデータを用いた。ここでの留意点は、特徴量の線型結合の絶対値が大きくなる(数十程度)と、シグモイド関数極値に振れてしまい、対数が計算できなくなることである。そのため、ここではあらかじめダミー変数以外の変数を平均0、標準偏差1で基準化した。

# ライブラリのインポート
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib # matplotlib日本語
# データの読み込み
df01 = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Social_Network_Ads.csv')
# 性別に関するダミー変数の作成
df02 = pd.get_dummies(df01, drop_first=True, columns=['Gender'])
# ダミー変数以外の特徴量を標準化する
df02 = df02.assign(salary_std = (df02['EstimatedSalary']-df02['EstimatedSalary'].mean())/df02['EstimatedSalary'].std())
df02 = df02.assign(age_std = (df02['Age']-df02['Age'].mean())/df02['Age'].std())
# 分析に用いるデータをNumpy配列に書き出す
features_list = ['salary_std','Gender_Male','age_std']
X_train = df02[features_list].to_numpy()
y_train = df02['Purchased'].to_numpy()
# バイナリ学習のアルゴリズム
w1 = 1; w2 = 1; w3 = 1; b = 1 # 回帰係数の初期化
alpha = 1 # 学習率
m = len(y_train)
z = np.zeros(m); a=np.zeros(m); dz=np.zeros(m)
for j in range(250): # 指定した回数の反復計算
  J = 0; dw1 = 0; dw2 = 0; dw3 = 0; db = 0 # 損失関数と勾配の初期化
  for i in range(m): # データ行についての処理
    z[i] = w1 * X_train[i,0] + w2 * X_train[i,1] + w3 * X_train[i,2] + b
    a[i] = 1 / (1+np.exp(-z[i]))
    J += -(y_train[i] * np.log(a[i]) + (1-y_train[i]) * np.log(1-a[i]))
    dz[i] = a[i] - y_train[i]
    dw1 += X_train[i,0] * dz[i]
    dw2 += X_train[i,1] * dz[i]
    dw3 += X_train[i,2] * dz[i]
    db += dz[i]
  J /= m; dw1 /= m; dw2 /= m; dw3 /= m; db /= m
  w1 -= alpha * dw1
  w2 -= alpha * dw2
  w3 -= alpha * dw3
  b -= alpha * db
  if j % 20 == 0:
    print('j={}, J={:.3g}, w1={:.3g}, w2={:.3g}, w3={:.3g},b={:.3g}'.format(j,J,w1,w2,w3,b))

損失関数Jの収束の様子から学習率αを調節していく。αを調節した結果は以下のようになった。収束には100回以上の反復計算を要している。
f:id:nicjps230:20210715140252p:plain

データサイエンスのおすすめオンライン記事(6月25日付)

緊急事態宣言はほぼ解除され、限定的ながら飲食店での酒類提供も復活した。当方にもようやくワクチン接種券が届いたが、現状で唯一可能性のある自衛隊の大規模接種も満杯で予約できない。さらにワクチン接種ペースが加速し、安心して活動できる日が一日も早く到来することを願うのみである。
今回もここ1~2ヶ月の間に掲載された新しい英文記事からピックアップして紹介する。

20倍早い格子探索交差検証法
towardsdatascience.com
機械学習におけるハイパーパラメータのチューニングで用いられる「格子探索交差検証法」の高速版である「二等分格子探索交差検証法」(halving grid search cross validation)についての解説。

初学者のためのトップデータサイエンスプロジェクト
www.kdnuggets.com
「データ収集」「探索的データ解析」「データ可視化」などのスキル別に、実データとチュートリアルからなる「プロジェクト」を紹介している。OJT以外でデータサイエンススキルを鍛えたい人には必見の良記事。

ズームイン:回路入門
distill.pub
顕微鏡で物体の細部を見るように、表面的にはブラックボックス化しているニューラルネットワークの結合と重みを理解する。

7ステップの機械学習プロセス
www.datasciencecentral.com
Data Science CentralサイトのオーナーであるVincent Granville氏の最新の投稿。機械学習プロジェクトの一般的なプロセスを要約したもの。同氏の卓越した見識で書かれた記事は読むに値する。

新技術はいかにウェブ開発に影響を及ぼしているか(2021年版)
www.datasciencecentral.com
人工知能bot」「加速モバイルページ(AMP)」「単一ページアプリケーション(SPA)」など最近のウェブ開発にまつわる新技術の動向をまとめたもの。

Pythonによるログデータの分析小技集

最近のビッグデータのうちの多くのものが、アクセスログ、イベントログなどのログデータです。その特徴として以下の2つがあります。
(1)多数の個体についての時系列データである(タイムスタンプを持つ)
(2)収集しているデータの大部分が(定型/非定型の)テキスト情報である
そこで今回はこのようなログデータをPythonで加工・分析するためのTip(小技・コツなど)をまとめました。
データの事例として、Kaggleにアップロードされたイベントログのデータ(https://www.kaggle.com/mehulkatara/windows-event-log) を用いました。
ここでは、「基本的なデータクリーニングや構文解析は終わっており、個体識別IDとタイムスタンプのマルチインデックスを持つPandasデータフレームとなっている」状態から始めます。
まずデータをデータフレームに読み込み、特定の個体のケースに絞り込みます。ここではdf.index.isin(キーワードのリスト)でインデックスに基づく検索を行っています。

import pandas as pd
df12 = pd.read_parquet('eventlog.parquet')
# 特定の個体のケースに限定する
ch_machine = 'LAPTOP-1MKMTVPM'
df2 = df12[df12.index.isin([ch_machine], level=0)]

定型テキストの出現頻度の推移をグラフ化する
次に、タイムスタンプから1時間毎といった「丸め時刻」を抽出し、ある定型テキストの出現頻度を丸め時刻毎に集計してグラフ化します。丸め関数は四捨五入のround()でなく切り捨てのfloor()を用います。演算によりデータフレームに新しい列を作る操作は、df['列名']への代入(警告が出ることがある)でなく、df.assign()メソッドを用いています。定型テキストの種別毎の合計度数が縦方向に並びますが、df.pivot()でそれらを横方向に展開しています。グラフ化は、pandasに用意されているdf['列名'].plot()メソッドを用いています。

# 丸め時刻毎に定型テキストの出現頻度を集計してグラフ化する
# 日単位の丸め時間変数を作成する
df21=df2.assign(time_round=df2.index.to_series().str[1].dt.floor('D'))
# 日とEntryType毎に行数をカウントする
df3 = df21.groupby(['time_round','EntryType']).count()
df31 = df3.reset_index(); df32 = df31.set_index('time_round')
# 各EntryTypeを列とする変数を作成する
df4 = df32.pivot(columns='EntryType',values='ch_machine_name')
df41 = df4.fillna(0) # 欠損値にゼロを代入する
import matplotlib.pyplot as plt # エラーの発生の推移をグラフ化する
from matplotlib import rcParams; rcParams['font.family']='MS Gothic' 
df41['Warning'].plot(figsize=(9,6)); df41['Error'].plot(figsize=(9,6))
plt.title('日別の警告とエラーの数   マシン: '+ch_machine)
plt.legend(); plt.show()

以下のようなグラフが表示されます。
f:id:nicjps230:20210620232110p:plain
時間範囲を指定して検索する
次は、タイムスタンプ情報をもとに、ある期間に入っているデータの検索を行う事例です。ここではタイムスタンプをインデックスから一般の列に書き出した上で、別途設定した日付値との比較演算を行っています。尚、タイムスタンプにタイムゾーン情報がある場合は、ここでのquery()の比較演算はエラーになりますが、その場合はdf['列名'].replace(tzinfo=None)によりタイムゾーンの情報を消しておきます。

# 時刻(範囲)を指定したデータの検索
import datetime
from datetime import timedelta
# 日付、時、時間範囲を指定する
query_date=20201023; query_time_h=8; query_hours=21
# 日付関連値を求める
c_year=query_date//10000
c_month=(query_date//100)%100
c_day=query_date%100
# データ開始時刻の日付値
dt1=datetime.datetime(c_year,c_month,c_day,query_time_h,0,0)
# データ終了時刻の日付値
dt2 = dt1 + timedelta(hours=query_hours)
# インデックスから時間変数を取り出す
df51 = df2.assign(c_time = df2.index.to_series().str[1])
 # 指定した時間範囲内のデータに絞り込む
df52 = df51.query('c_time >= @dt1 & c_time <= @dt2')
df53 = df52.loc[:,['EntryType','Message','Source']]

この結果、以下のように行が絞り込まれたデータフレームが得られます。
f:id:nicjps230:20210620234301p:plain
定型テキストの出現傾向を表に要約する
数値データの場合はある時間間隔毎の集計値を見て傾向をつかみますが、テキストデータの場合は別の要約方法が必要となります。ここでは、ある時間間隔毎に、最も頻繁に出現した文字列をその時間間隔における「代表文字列」として、それを表の形式に集約し、テキストの出現傾向を把握しやすいようにします。ここでは(系列名).value_counts()にて文字列の出現頻度表を取得し、さらに.index[0]で最頻度の文字列を抽出しています。またここでもdf.pivot()を用いてデータを横方向に展開し、さらにmatplotlibのtable関数を用いて表を作成しています。メッセージなどのテキスト情報は一般的に長いので、表に表示するにあたって、文字数の切り詰めと改行の強制挿入を実施しています。

# 丸め時間区間毎に最も多く登場した定型テキストを表形式で表示する
import textwrap
# strftime()で日本語を使用するのに必要
import locale; locale.setlocale(locale.LC_ALL, '')
def func_part_of_month(tround):
    '''月の上旬・中旬・下旬の判定関数'''
    if tround.day <=10 :
        return '上旬'
    elif tround.day <= 20 :
        return '中旬'
    else :
        return '下旬'
df61 = df2.reset_index()
# 年月を表す変数を導入
df62=df61.assign(yr_month=df61['timestamp'].dt.strftime('%Y年%m月'))
# 上/中/下旬を表す変数を導入
df62=df62.assign(part_of_month=df62['timestamp'].apply(func_part_of_month))
# グループ毎に最頻出の文字列を求める
df63=df62.groupby(['yr_month','part_of_month'])['Message'].apply(lambda x:x.value_counts().index[0])
df64 = pd.DataFrame(df63)
def func_wrap(s):
    '''長い文字を切り詰め、途中に改行を入れる関数'''
    s_trunc = textwrap.shorten(s,100)
    s_wrap_list = textwrap.wrap(s_trunc, 25)
    return '\n'.join(s_wrap_list)
# 文字列切り詰めと改行挿入
df64['f_msg'] = df64['Message'].apply(func_wrap)
df65=df64.reset_index(); df66 = df65.set_index('yr_month')
# 各part_of_monthを列とする変数を作成する
df67 = df66.pivot(columns='part_of_month',values='f_msg')
# 欠損値の表示方法の指定
df68 = df67.fillna('記録なし')
# データの行数に応じて表の高さを決める
tbl_height = len(df68) * 1.0
# 図形の描画
fig=plt.figure(figsize=(8,tbl_height))
ax=fig.add_subplot(111); ax.axis('off')
# matplotlibのtable関数で表を描く
tbl=ax.table(cellText=df68.values,bbox=[0,0,1,1],colLabels=df68.columns,rowLabels=df68.index)
plt.plot()

以下のような表が出力されます。
f:id:nicjps230:20210620235235p:plain