20201208のPythonに関する記事は30件です。

Git/GitHubの初期操作に触れてみた(@Python/Django)

はじめに

今回は、Python/DjangoにてWebアプリを作成する過程で扱う
Git/GitHubの初歩的な初期操作について簡潔に確認していきます。

【前工程】
・djangoプロジェクト作成はこちら
・Git初期設定はこちら

GitHubの初期操作

基本的には以下のように、公式のsetupに推奨される手順通り、コマンドを実行していきます。これらを実行することで、ローカルのプロジェクトを保存することから、Web上のGitHubへUPするところまで行います。

※事前にGitHubにサインインし、「+」マークから「New repository」を選択し、任意のリポジトリ名を付けて「Create」しておきます。

①リポジトリの作成

$ git init

②コミット対象を指定
 今回の場合、"."とすることで、カレントディレクトリ以下の全てのファイルを対象としています。

$ git add .

③ローカルに保存する
 "first commit"という名前でローカルに保存します。

$ git commit -m "first commit"

④ローカルのブランチ名を変更する

$ git branch -M main

⑤事前に作成しておいたGitHubのリポジトリを使用してリモートリポジトリを追加する

$ git remote add origin https://github.com/指定のリポジトリ名/

誤って追加した場合の削除は以下の通りです。

$ git remote rm origin

⑥GitHub上にプッシュする(Web上にUPする)

$ git push -u origin main

ここでうまく実行されない場合は、
初期設定がうまくいってない可能性があるため、
再確認し再度実行してみます。

GitHub上で"Success"と表示され完了しました。

まとめ

今回は、Webアプリを作成する際、事前のプロジェクト作成後のGit/GitHub初期操作を確認しました。

Python/djangoによるWebアプリ作成の過程で、さらにGit/GitHubについて掘り下げていきたいと思います。

ありがとうございました。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Git/GitHubに触れてみた(@Python/Django)

はじめに

今回は、Python/DjangoにてWebアプリを作成する過程で扱う
Git/GitHubの初歩的な初期操作について簡潔に確認していきます。

【前工程】
・djangoプロジェクト作成はこちら
・Git初期設定はこちら

GitHubの初期操作

基本的には以下のように、公式のsetupに推奨される手順通り、コマンドを実行していきます。これらを実行することで、ローカルのプロジェクトを保存することから、Web上のGitHubへUPするところまで行います。

※事前にGitHubにサインインし、「+」マークから「New repository」を選択し、任意のリポジトリ名を付けて「Create」しておきます。

①リポジトリの作成

$ git init

②コミット対象を指定
 今回の場合、"."とすることで、カレントディレクトリ以下の全てのファイルを対象としています。

$ git add .

③ローカルに保存する
 "first commit"という名前でローカルに保存します。

$ git commit -m "first commit"

④ローカルのブランチ名を変更する

$ git branch -M main

⑤事前に作成しておいたGitHubのリポジトリを使用してリモートリポジトリを追加する

$ git remote add origin https://github.com/指定のリポジトリ名/

誤って追加した場合の削除は以下の通りです。

$ git remote rm origin

⑥GitHub上にプッシュする(Web上にUPする)

$ git push -u origin main

ここでうまく実行されない場合は、
初期設定がうまくいってない可能性があるため、
再確認し再度実行してみます。

GitHub上で"Success"と表示され完了しました。

まとめ

今回は、Webアプリを作成する際、事前のプロジェクト作成後のGit/GitHub初期操作を確認しました。

Python/djangoによるWebアプリ作成の過程で、さらにGit/GitHubについて掘り下げていきたいと思います。

ありがとうございました。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

cronから仮想環境venvのスクリプトを実行する

https://github.com/kazuho/kaztools/blob/master/cronlog
cron で > /dev/null して椅子を投げられないための3つの方法 - 酒日記 はてな支店

$ cd tmp
$ python3 -m venv venv
$ pip install requests
$ cat p.py
import requests
$ cat /etc/cron.d/test
*/1 * * * * k8uwall cd /home/k8uwall/tmp; /usr/local/bin/cronlog --timestamp -- venv/bin/python p.py

# tail -f /var/log/cron
Dec  8 21:53:01 k8uwallsv CROND[6402]: (k8uwall) CMD (cd /home/k8uwall/tmp/c/; /
usr/local/bin/cronlog --timestamp -- venv/bin/python p.py)
  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Mac版Arduino IDE で M5Stackがビルドできない(import Serialエラー)

まとめ

  • M5Paperをゲットし、MacでArduinoIDEのボードマネージャ(M5Stackライブラリ: v1.0.6)でM5Paperを選んでビルドしたところ、以下の問題に遭遇した。
    • Serial接続関連のパッケージがインストールされてないよ、と言う事らしい。
    • これはM5Stackライブラリに含まれるesptoolとMac環境の問題なので、このバージョンでM5Paper以外のボードを選択しても同じくビルドが通らなくなる。
  • Python3を何度かインストールし直したり、色々やった挙句、きれいに解決したのでここにメモしておく。
Traceback (most recent call last):
  File "/Users/nabeshin/Library/Arduino15/packages/m5stack/tools/esptool_py/3.0.0/esptool.py", line 39, in <module>
    import serial
ImportError: No module named serial
exit status 1
ボードM5Stack-Paperに対するコンパイル時にエラーが発生しました。

スクリーンショット 2020-12-06 23.07.50.png

実行環境

  • OS: MacOS Catalina
  • ArduinoIDE: v1.8.13
  • M5Stackライブラリ: v1.0.6
    • M5Paperが含まれているライブラリは以下のボードマネージャのURLから入手した。

スクリーンショット 2020-12-06 23.11.55.png

Python環境

nabeshin@iMacNabeshin ~ % python2

WARNING: Python 2.7 is not recommended. 
This version is included in macOS for compatibility with legacy software. 
Future versions of macOS will not include Python 2.7. 
Instead, it is recommended that you transition to using 'python3' from within Terminal.
  • MacOSでは、初めからPython2がインストールされているが、Python3への移行が推奨されているようなので別途python3をインストールしてある。
    • ~/.zprofile でpython3にパスを通し、デフォルトにしている。
nabeshin@iMacNabeshin ~ % env | grep PATH 
PATH=/Library/Frameworks/Python.framework/Versions/3.9/bin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin
nabeshin@iMacNabeshin ~ % python --version
Python 3.9.0
nabeshin@iMacNabeshin ~ % 
  • python2で>>import serialを打って見たところ、ImportError: No module named serialとなる。
  • >pip install serialしてあるputhon3だとエラーにならない。

Arduino IDEは強制的にPython2を使ってしまう

  • シェルではバッチリpython3にパスが通っており、>>import serialも成功するが、Arduino IDEでビルドすると、どうしてもPython2が使われてしまう。

検証その1

  • 以下のようにpython3へのシンボリックリンクを張ってみたが、変わらず。。。
nabeshin@iMacNabeshin 3.0.0 % pwd
/Users/nabeshin/Library/Arduino15/packages/m5stack/tools/esptool_py/3.0.0
nabeshin@iMacNabeshin 3.0.0 % ls -la
total 7992
drwxr-xr-x  6 nabeshin  staff      192 12  7 21:37 .
drwxr-xr-x  4 nabeshin  staff      128 11 28 23:58 ..
-rw-r--r--@ 1 nabeshin  staff     6148 11 29 17:23 .DS_Store
-rwxr-xr-x  1 nabeshin  staff  3910331 11 28 21:47 esptool
-rwxr-xr-x@ 1 nabeshin  staff   169633 12  6 23:02 esptool.py
lrwxr-xr-x  1 nabeshin  staff       22 12  7 21:37 python -> /usr/local/bin/python3
nabeshin@iMacNabeshin 3.0.0 % 

検証その2

  • 以下のようにesptool.pyにpythonのバージョンをファイルに書き出すようにしてみた。

スクリーンショット 2020-12-06 23.07.02.png

  • 結果
    • だめだ、やはり2.7が強制的に使わされてしまう。。。
    • この呪縛から逃れる方法はないものか。。。
nabeshin@iMacNabeshin tmp % pwd
/Users/nabeshin/tmp
nabeshin@iMacNabeshin tmp % cat DEBUG_Version.txt
2.7.16 (default, Jun  5 2020, 22:59:21) 
[GCC 4.2.1 Compatible Apple LLVM 11.0.3 (clang-1103.0.29.20) (-macos10.15-objc-
nabeshin@iMacNabeshin tmp % 

と言うわけで、
なんとかしてpython2側にserialパッケージを追加しなければならない。と言うことがわかったが、
> pip2 install serial を素直にやりたくてもpip2がインストールできない。。。

解決方法

  • 全く同じ問題について、公式のM5Stack Communityで話されていた。
  • 解決した、とされている手順に沿って以下を実施。
  • これにより、Arduino IDEでもビルドできるようになった。
nabeshin@iMacNabeshin ~ % python2

WARNING: Python 2.7 is not recommended. 
This version is included in macOS for compatibility with legacy software. 
Future versions of macOS will not include Python 2.7. 
Instead, it is recommended that you transition to using 'python3' from within Terminal.

Python 2.7.16 (default, Jun  5 2020, 22:59:21) 
[GCC 4.2.1 Compatible Apple LLVM 11.0.3 (clang-1103.0.29.20) (-macos10.15-objc- on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> print sys.path
['', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/Library/Python/2.7/site-packages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC']
>>> ^D
nabeshin@iMacNabeshin ~ % sudo pip3 install --target /Library/Python/2.7/site-packages pyserial
Password:
WARNING: The directory '/Users/nabeshin/Library/Caches/pip' or its parent directory is not owned or is not writable by the current user. The cache has been disabled. Check the permissions and owner of that directory. If executing pip with sudo, you may want sudo's -H flag.
Collecting pyserial
  Downloading pyserial-3.5-py2.py3-none-any.whl (90 kB)
     |████████████████████████████████| 90 kB 8.8 MB/s 
Installing collected packages: pyserial
Successfully installed pyserial-3.5
nabeshin@iMacNabeshin ~ % 






nabeshin@iMacNabeshin ~ % python2

WARNING: Python 2.7 is not recommended. 
This version is included in macOS for compatibility with legacy software. 
Future versions of macOS will not include Python 2.7. 
Instead, it is recommended that you transition to using 'python3' from within Terminal.

Python 2.7.16 (default, Jun  5 2020, 22:59:21) 
[GCC 4.2.1 Compatible Apple LLVM 11.0.3 (clang-1103.0.29.20) (-macos10.15-objc- on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import serial
>>> serial
<module 'serial' from '/Library/Python/2.7/site-packages/serial/__init__.py'>
>>> 

pip3からpython2側にパッケージを追加するコマンド

  • > sudo pip3 install --target /Library/Python/2.7/site-packages pyserial
    • なるほど、こうやって実行すればいいのか。
    • よく覚えておきます。
  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

DNNにおけるClassificationの種類まとめ

概要

"AI"とかいうやつの課員用教育ネタの続き
深層学習の出口レイヤを変えると簡単に問題を切り替えられることを自習させようかと…
備忘録も兼ねて残すこととした

  • 実施期間: 2020年12月
  • 環境:Keras

ここで説明する各モデルの出口の全結合レイヤについて記述する
それより前のレイヤ構成はケースバイケース

Number of Nodes : ノード数 [1]
Activation : 活性関数 [2]
Loss : ロスの算出方法 [3]
Accuracy metrics: Fit中の精度評価指標 [4]
Output: 出力の型

Kerasだと、それぞれ下記の引数がそれらとなる
Optimizerはadamである必要はない
Dense, compileはそれぞれオフィシャルサイト参照のこと

model.add(Dense([1], activation='[2]'))
model.compile(loss='[3]', optimizer='adam‘, metrics=['[4]'])

Regression

連続値である目的変数yを予測するモデルで回帰問題と呼ばれる

Regression

基本形がこれ
Training済みモデルにXを入力すると予測値のy_hat(ここでは0.539)が出力される
KerasではDenseの引数にActivationを指定しなければデフォルトでLinearとなる

image.png

Multi-output regression

入力Xに対して複数のyがあるモデルがこれ
出口の全結合レイヤのNode数はそのyの数(ここでは3つ)でTrainingさせておく
image.png

Classification

離散値である目的変数yを予測するモデルで分類問題と呼ばれる
モデル出口とyの型以外はRegressionと変わらず、ビビるに及ばない
以下に説明する3タイプは日本語でなんと呼ぶか知らない

Binary classification

Training済みモデルにXを入力すると予測値のy_hat(ここでは1(True))が出力される
yが連続値から離散値になり、出口の全結合レイヤの引数がやや変わっただけ
あとは上述のRegressionと同じ
不良品か否かを当てるようなモデル
image.png

Multilabel binary classification

Multi-output regressionとほぼ同じ
実際は各Nodeから出力されるのは連続値で、これを0(False)か1(True)に丸める
ただ小生はHeuristicに、例えば0.3以上ならTrueとか閾値は柔軟にしてよいと思っている
Training時のy(例えば[1, 1, 0])を各Nodeで予測するので、それぞれの出力は 0 <y_hat< 1 となる
image.png

Multiclass classification

これだけ少し毛色が違う
Activationがsoftmaxとなっており、全Nodeからの出力中Trueは1つである
実際は各Nodeから出力されるのは確率値となり、従い全Nodeの出力合計は1.0になる
最大のものを1(True)、それ以外を0(False)とする
犬、猫、人の画像からそのどれかを当てるようなモデルに使う
image.png

まとめ

モデル # of nodes Activation Loss Accuracy metrics Output
Regression 1 Linear MSE, etc. accuracy 連続値
Multi-output regression 複数 Linear MSE, etc. accuracy 連続値
Binary classification 1 sigmoid binary_crossentropy binary_accuracy 離散値(True/False)
Multilabel binary classification 複数 sigmoid binary_crossentropy binary_accuracy 離散値(複数True/False)
Multiclass classification 複数 softmax categorical_crossentropy categorical_accuracy 離散値(単数True/False)

“Class”と”Label”の違い

下図の場合、Classは3つありLabelは[犬, 猫, 人]
image.png

以上

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Pytorchでパディングして順伝搬計算まで(LSTM)

はじめに

こんにちは、とあるIT企業の人間です。
自然言語処理の勉強のためにPytorchを使用して、LSTMモデルを実装していたところ、「お、系列長バラバラやけど、どうすんの?」と思い、自分なりに上手く書けたのでこの記事を書きました。
あくまでも我流の一手法のため、参考にしていただけるとありがたいです。

なお、関数などの詳細は公式ページを参考にしてください。

何かご指摘・質問等ございましたら、ご遠慮なくコメント欄にお願い致します。

目的

自然言語をIDで表現した状態のデータがある場合に以下の処理をスマートに実装したい
1. 系列長がバラバラのため、パディングする。
2. Embedding層、LSTM層に入力する(順伝搬計算)
3. それぞれの系列長の最終出力を抽出する。

つまり以下のような系列長が違うデータをスマートな実装で一気に処理したい、ということです。
データ例:
inp = [[ 7, 3, 5], [ 9, 5], [ 7, 1, 8, 4, 4]]

環境

  • Python 3.7.6
  • pytorch 1.5.0

実装

データ

系列長がバラバラのリストと、それぞれの系列長が入ったリストを準備します。
リストで各要素をtensor型に変換してください。今回では以下のデータを準備します。

# 使用データの定義
word_id = [T.tensor([ 7, 3, 5]),        # 3
           T.tensor([ 9, 5]),           # 2
           T.tensor([ 7, 1, 8, 4, 4])]  # 5

こちらの系列長は、$[3,2,5]$です。こちらもリストで用意します。

# 使用データの定義
seq_list = [3, 2, 5]

1. 系列長がバラバラのため、パディングする。

pytorchには、パディングしてくれる関数pad_packed_sequenceがあるので、そちらを使いますが、その入力のためにまずpack_sequenceを使用します。

# パディング
## ①PackedSequenceという一つのオブジェクト型に変換する。
word_id_packed = nn.utils.rnn.pack_sequence(word_id, enforce_sorted=False)
## ②0埋めする
word_id_padded, _ = nn.utils.rnn.pad_packed_sequence(word_id_packed,batch_first =True,padding_value=0)

word_id_paddedの出力結果です。
お、0埋めできてますね。

tensor([[7, 3, 5, 0, 0],
        [9, 5, 0, 0, 0],
        [7, 1, 8, 4, 4]])

「0以外でパディングしたい」「最大系列長を指定したい」と思う方は公式ページに引数の説明があるのでそちらを見ていただくと良いかと思います。

2. Embedding層、LSTM層に入力する(順伝搬計算)

ここでは最終的に各系列長の出力が欲しいので、パディングした状態を保ったままLSTM層に入力する必要があります。そのため、LSTM層への入力前にpack_padded_sequenceという関数を使用します。

また「"0"はパディングしたもの」として扱いたいので、nn.Embeddingpadding_idx=0を指定します。"-1"など他の値でパディングした場合は、padding_idxの値を変更してください。

# 各層のパラメータ定義
## Embedding層(辞書の単語数:10、次元数:5)
emb_layer = nn.Embedding(10, 5, padding_idx=0)
## LSTM層(入力次元数:5、出力次元数:2)
lstm_layer = nn.LSTM(5, 2) 


# 順伝搬計算
## Embedding層へ入力
word_emb = emb_layer(word_id_padded)
## パディングした状態を保持したまま、PackedSequenceオブジェクト型に変換する。
word_emb_packed = nn.utils.rnn.pack_padded_sequence(word_emb, seq_list, batch_first =True, enforce_sorted=False)
## LSTM層へ入力する。この時、第二返り値も受け取るようにセットする。(ht, ct)
lstm, (ht, ct)  = lstm_layer(word_emb_packed)

lstmの出力結果です。
各系列長分の計算だけができていますね。Good。

tensor([[[ 0.2595,  0.0771],
         [ 0.2072,  0.2954],
         [ 0.3157,  0.0541],
         [ 0.0000,  0.0000],
         [ 0.0000,  0.0000]],

        [[ 0.3509,  0.2226],
         [ 0.2906,  0.0636],
         [ 0.0000,  0.0000],
         [ 0.0000,  0.0000],
         [ 0.0000,  0.0000]],

        [[ 0.2595,  0.0771],
         [ 0.2382,  0.0205],
         [ 0.2692,  0.1638],
         [ 0.2793,  0.0356],
         [ 0.3055, -0.0184]]]
          , grad_fn=<IndexSelectBackward>)

3. それぞれの系列長の最終出力を抽出する

最後に先ほどのhtを見てみましょう。

print(ht[-1])

出力結果です。これだけ見ても「お?」という感じなので先ほどの出力を見てみましょう。

tensor([[ 0.3157,  0.0541],
        [ 0.2906,  0.0636],
        [ 0.3055, -0.0184]], grad_fn=<SelectBackward>)

おお!!lstmの出力からうまく抽出できてますね!!

tensor([[[ 0.2595,  0.0771],
         [ 0.2072,  0.2954],
         [ 0.3157,  0.0541],  <------
         [ 0.0000,  0.0000],
         [ 0.0000,  0.0000]],

        [[ 0.3509,  0.2226],
         [ 0.2906,  0.0636],  <------
         [ 0.0000,  0.0000],
         [ 0.0000,  0.0000],
         [ 0.0000,  0.0000]],

        [[ 0.2595,  0.0771],
         [ 0.2382,  0.0205],
         [ 0.2692,  0.1638],
         [ 0.2793,  0.0356],
         [ 0.3055, -0.0184]]] <------
          , grad_fn=<IndexSelectBackward>)

全体のコード

個人的にはすっきりしました。

padding_lstm.py
# 使用データの定義
word_id = [T.tensor([ 7, 3, 5]),        # 3
           T.tensor([ 9, 5]),           # 2
           T.tensor([ 7, 1, 8, 4, 4])]  # 5

seq_list = [3, 2, 5]


# パディング
## ①PackedSequenceという一つのオブジェクト型に変換する。
word_id_packed = nn.utils.rnn.pack_sequence(word_id, enforce_sorted=False)
## ②0埋めする
word_id_padded, _ = nn.utils.rnn.pad_packed_sequence(word_id_packed,batch_first =True,padding_value=0)


# 各層のパラメータ定義
## Embedding層(辞書の単語数:10、次元数:5)
emb_layer = nn.Embedding(10, 5, padding_idx=0)
## LSTM層(入力次元数:5、出力次元数:2)
lstm_layer = nn.LSTM(5, 2) 


# 順伝搬計算
## Embedding層へ入力
word_emb = emb_layer(word_id_padded)
## パディングした状態を保持したまま、PackedSequenceオブジェクト型に変換する。
word_emb_packed = nn.utils.rnn.pack_padded_sequence(word_emb, seq_list, batch_first =True, enforce_sorted=False)
## LSTM層へ入力する。この時、第二返り値も受け取るようにセットする。(ht, ct)
lstm, (ht, ct)  = lstm_layer(word_emb_packed)

# 
print(ht[-1])

参考

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

ネットワーク分析3 発展編

Aidemy 2020/12/8

はじめに

 こんにちは、んがょぺです!バリバリの文系ですが、AIの可能性に興味を持ったのがきっかけで、AI特化型スクール「Aidemy」に通い、勉強しています。ここで得られた知識を皆さんと共有したいと思い、Qiitaでまとめています。以前のまとめ記事も多くの方に読んでいただけてとても嬉しいです。ありがとうございます!
 今回は、ネットワーク分析の3つ目の投稿になります。どうぞよろしくお願いします。

*本記事は「Aidemy」での学習内容を「自分の言葉で」まとめたものになります。表現の間違いや勘違いを含む可能性があります。ご了承ください。

今回学ぶこと
・JSONファイルの扱い方
・実際のデータを使ってネットワーク分析

JSONファイル

JSONファイルとは

・Chapter3では実際のデータを使ってネットワーク分析を行う。pythonで実際のデータを取得する方法の一つとして、「JSONファイルを読み込む」ということがある。
・JSONファイルを読み込むには、「open()」関数を使ってファイルを開き、「load()」関数を使う必要がある。こうすることで、python上では辞書型のオブジェクトとして扱うことができるようになる。
・以下では、人間データが入ったjsonファイルについて、年齢をkey名前をvalueとした辞書に再格納するものである。

・コードスクリーンショット 2020-12-05 12.22.51.png

実際のデータを使ってネットワーク分析

データの読み込みと整理

・ここでは、A_社の社内SNSの記録をまとめたデータ「sample_data.json」を使って、このファイルの読み込み、整理を行なっていく。具体的に、このデータには「A氏の番号:{B氏の番号:値,C氏の番号...}」と格納されている。例えば「"1":{"3":2}」と格納されていた場合、これは「1氏が3氏に送ったメンション数が2回である」ことを示す。
・このデータを、次項以降でグラフにしたいのであるが、この時、
人を頂点、人同士の繋がりを辺とする。このようにする場合、前処理として、「誰とも繋がりのない人」を除外しておく。読み込んだjsonファイルを「data」とすると、その各keyについて、data[key]の長さが0であれば、誰にもメンションしていないということになるので、これを抽出し、「del」_で削除する。

・コードスクリーンショット 2020-12-05 12.46.42.png

データの保存

・前項で修正したデータを保存する。この場合は、「open()」関数で「'w'」と指定して書き込み専用で呼び出す必要がある。ちなみに、この時存在しないファイル名を呼び出したときは、そのファイルが作成される。次に「json.dump()」で、辞書型からjson型の文字列に変換することで、書き込むことができる。
・json.dump()の各引数は次のようになる。「ensure_ascii」にはTrueorFalseを指定し、入力された文字列をそのまま出力するかを指定する。「indent」には数値を指定し、インデントの字数を指定する。「sort_keys」もTorFで、辞書の出力にkeyでソートするかを指定する。「separators」(',',':')と指定することで、「keyとvalueのセットごとに「,」で区切る」「keyとvalueは「:」で区切る」ことを指定する。

・コードスクリーンショット 2020-12-05 22.43.30.png

グラフの描画

・ここでは、重み付きグラフを作成する。辺の重みがメンション数となればよい。
・まずは、頂点がdataのkeyがstr型になっているので、dataをint型に変える。その後、この二重辞書の、keyの中のkey(sub_key)を指定することで、メンション数を取得できる。これを「G.add_edge()」の引数「weight」に指定することで、重みを指定できる。辺を結ぶ頂点については、そのkey(発信者)とsub_key(受信者)とすれば良い。
・あとは、この重みを辺の太さにして、グラフを描画する。辺の太さは「nx.draw_networkx_edges()」の引数「width」に指定すれば良い。G.edges()で取得できる「頂点1,頂点2,重み」の3つ目の「weight」をwidthとし、それを基準に辺の太さを決定する。

・コードスクリーンショット 2020-12-05 14.29.52.png

・結果スクリーンショット 2020-12-05 14.32.38.png

グラフの見た目を整える

・次に、頂点の大きさを、固有ベクトル中心性の大きさに応じて変化させる。固有ベクトル中心性は「nx.eigenvector_centrality_numpy()」で算出できる。また、今回のようにedgeに「weight」という重みを与えている場合は、引数に「weight="weight"」と指定する。これを与えない場合、edgeは全て等しいものとみなされる。
・ノードの大きさは、今回は、固有ベクトル中心性「eigv_cent」の各値をリストとし、それに10000をかけたものとする。これを「nx.draw_networkx_nodes()」の引数「node_size」に指定する。

・コードスクリーンショット 2020-12-05 19.13.00.png

・結果スクリーンショット 2020-12-05 19.14.07.png

クラスタリングを行う

・ここではクラスタリングを行う。前Chapterの復習であるが、クラスタリングは連結されたグラフをいくつかのサブグラフに分けることであり、コードとしては「community.best_partition(G)」で行える。これにより、各頂点が何番目のコミュニティに属しているのかが辞書型で返される。
・上記コードについては、今までのコードに付け加えるだけで良い。また、ここでもコミュニティごとに色分けを行って描画する。コミュニティはクラスタリングを行った「partition」の、ノードごとの値を参照すれば良い。

・コード(その他の部分は前項までと同じ)スクリーンショット 2020-12-05 19.39.13.png

・結果スクリーンショット 2020-12-05 19.38.51.png

中心性が高い人物の抽出

・ここでは、各コミュニティごとに最も固有ベクトル中心性が高い人物を抽出する。前項の方法でpartitionを作成し、このvalue(コミュニティ)について、その番号に位置に空のリストがある、(二重)リスト「part_com」を作成する。わかりにくいがつまり「[[コミュニティ0のnode],[コミュニティ1のnode],...]」のようになるリストである。そして、この空のリストの中に、コミュニティごとのnodeを入れていく。nodeはpartitionのkeyに格納されているので、これを入れる。
・そして、それぞれのコミュニティのリストごとに固有ベクトル中心性の最大値の頂点を求める。まずはグラフを初期化し、新たなグラフ「G_part」を作成する。固有ベクトル中心性の算出にはnode同士の繋がりが必要であり、辺を構成する2つのnodeが同じコミュニティにあることが前提条件となるので、edgeの両端(edge[0]とedge[1])が同じコミュニティ(part)内にある場合に、重みを含んだ辺edgeをグラフに追加する。
・そして、このグラフのnodeについて、固有ベクトル中心性をkeyとして、最大値を取得する。

・コードスクリーンショット 2020-12-05 22.07.29.png

コミュニティの指定とモジュラリティ指標

・前項までで扱った「community.best_partition()」は、引数「partition」に、「キーが頂点番号で、値がコミュニティの番号」である辞書を指定すると、この辞書を参考にしてコミュニティが形成される。
・また、「community.modularity(partition,G)」で、良い分割かを表す指標「モジュラリティ指標」を算出できる。
・以下では、コミュニティの分け方を指定しないでクラスタリングを行う「partition1」と、頂点番号が1桁のものをコミュニティ「0」、2桁のものを「1」としてコミュニティを分ける「partition2」について、それぞれのモジュラリティ指標を算出している。

・コードスクリーンショット 2020-12-05 22.31.09.png

まとめ

・実際のデータでネットワーク分析を行う場合も、扱うファイルがJSONファイルであるなどの違いを除けば、基本的にはkarateclubの時と変わらない
・今回のJSONデータは、pythonで扱えるように読み込むと、{A:{B:value}}のような二重辞書になっている。そのため、中の辞書のkeyは「sub_key」として取り出すと、データを楽に扱える。また、グラフの描画についても、中の辞書のvalue(メンション数)を重みとして、辺を作成すれば良い。
クラスタリングを行うことで、各頂点がどのコミュニティに属するかを知ることができる。これを使って色分けしたり、各コミュニティごとの中心性が最大の頂点を抽出したりできる。

今回は以上です。最後まで読んでいただき、ありがとうございました。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

漫画データセットで遊ぶ2

ご挨拶

こんにちは、マンボウです。
漫画データセット1で遊ぶに引き続き、kaggleのデータを分析していきます。

使うデータ(漫画レビューサイトのデータを前処理した結果)

TOP RANKED MANGAS MyAnimeList (MAL)
スクリーンショット 2020-11-26 11.20.42.png

どのジャンルが付いていると高いスコアが付きやすい?

海外の方の好みをみてみると、例えば以下のような感想です。

  • 「音楽」漫画は安定して非常に評価が高い
  • 「侍」漫画は概ね評価が高いが玉石混交感

スクリーンショット 2020-12-08 20.57.01.png
参考までに、作品数です。
スクリーンショット 2020-12-08 21.00.16.png
音楽ジャンルの作品一覧も出してみました。
スクリーンショット 2020-12-08 21.07.14.png

どういうジャンルの組み合わせだとスコアが高くなる?

実は各漫画、複数のジャンルを持ってます。
なので、上記のように単体ジャンルでの評価というよりかは、組み合わせの視点も必要かと思いました。

  • スコア上位50作品
  • スコア下位50作品

のジャンルの組み合わせ傾向を、アソシエーション分析の結果で見てみます。

スコア上位のジャンル組み合わせ傾向

「青年」「ドラマ」とが掛け合わさると強力のようです
スクリーンショット 2020-12-08 21.13.17.png
ヒットする作品TOP10はこんな感じ。
スクリーンショット 2020-12-08 21.32.49.png

スコア下位のジャンル組み合わせ傾向

「少年」「アクション」「ファンタジー」弱いですね。。。
スクリーンショット 2020-12-08 21.13.47.png
「少年」「アクション」にヒットするワースト10はこんな感じ。
スクリーンショット 2020-12-08 21.39.36.png

感想

ジャンル単一で見ても、組み合わせで見ても、日本人とちょっと違う感性が垣間見えます。
特に、リアリティがある漫画のほうが、海外では好まれるのかもしれない、と感じる部分がありました。
このあたりに考察にお詳しい方は、ぜひ感想を教えて下さい。

全部ノンプログラミングで実施

これらは全部、プログラミング不要の分析ツールnehanで実施しました。
下記実際にデータ処理を行ったノードプロセスです。
無料トライアルをしたい!興味がある!場合はぜひお問い合わせください!

nehan分析ワークフロー

スクリーンショット 2020-12-08 21.44.39.png

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

顔認証AIのFaceNetは顔のどこに注目しているのか

はじめに

この記事は顔学2020アドベントカレンダーの6日目の記事です.
前回投稿したFaceNetを使って,顔認証AIが顔のどこに注目して人を判別しているのかを可視化してみます.

FaceNet

顔認証のための顔認識モデルです.512次元の空間上に顔を埋め込むいわゆるEmbedding Modelで,空間上のユークリッド距離を使って顔の類似度を計算することで顔認証を行います.

前回の記事でFaceNetの使い方を紹介しています.

GradCAM

画像認識モデルがどこに注目をしているのかをヒートマップで可視化してくれるモデルです.画像の顕著性(Saliency)の可視化と似たような分野です.画像の顕著性については人間の視覚特性(Treismanの特徴統合理論)を根拠に画像処理を行うのですが,GradCAMは画像認識モデルの出力を参考にします.

最終層の出力は基本的にクラスラベルになっているので,クラスラベルの予測値の中で最も大きな値を認識結果とします.GradCAMでは,最も大きい値の出力(予測値)を計算する上で,寄与の大きかった画像箇所を逆算することで注目箇所を推定します.

スクリーンショット 2020-12-08 204620.png
スクリーンショット 2020-12-08 204733.png

FaceNet(Embedding Model)にどう導入するか

FaceNetの出力は空間上の座標なので,何か特定のクラスを表す認識結果を出力するわけではありません.なので,どの出力の値が重要なのか判断するのが非常に難しいです.

すごく厳密に寄与した出力チャネルを決める研究もあるみたいなんですが,ここでは簡単のために全出力の絶対値の中で最も大きかった値を採用しました.

Adapting Grad-CAM for Embedding Networks

DCGANInterFaceGANなどの論文が示唆するように,深層学習のモデルが構築する中間表現(潜在変数)には意味的軸があり,また現実世界の法則を反映したような軸同士の相関関係を持っています.(年齢を変化させると自然と眼鏡が出てきたり,白髪になったりする.つまり,意味的軸同士のcos類似度が高い.)

無題.png

FaceNetも多くの顔画像をきれいに分布させる顔空間を構築しているはずです.なので,なんらかの意味的定量化がなされた空間上であれば,各要素にも固有の意味を含んでいると仮定し,各座標値の絶対値を評価すればそれなりの結果になるだろうと見当をつけました.

実装(コード)

動作環境はGoogle Colaboratoryです.前回のノートブックに追記で書いています.

データ

前回の記事で使った首相データを再利用します.

ライブラリのインストール

!pip install tf-explain tensorflow==2.0.0 mira keras-facenet

FaceNetで顔ベクトルを取得

from mira.detectors import MTCNN
from keras_facenet import FaceNet
from PIL import Image
import numpy as np
import cv2
import matplotlib.pyplot as plt
import os

detector = MTCNN() # 顔領域検出器
embedder = FaceNet() # FaceNetのモデルを持つクラス

img = img = cv2.imread(FILE_PATH) # 画像読み込み
img_rgb = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) # RGB形式に変換
face = detector.detect(img_rgb)[0] # 画像中に複数の顔が検出されることもある.先頭をとってくる.
embedding = embedder.embeddings([face.selection.extract(img_rgb)])

GradCAMで注目箇所を可視化

from tf_explain.core.grad_cam import GradCAM

# 
model = embedder.model # FaceNetのモデル(KerasのModelクラス)
model.summary()

# 前処理
img = face.selection.extract(img_rgb)
img = cv2.resize(img, (160, 160))
X = np.float32([embedder._normalize(img)])
data = (X, None)

# 出力の重要要素を決定
abs = np.abs(embedding) # 絶対値
top_channel = np.argmax(abs) # 絶対値が最も大きい要素番号

# GradCAMで可視化
explainer = GradCAM()
grid = explainer.explain(data, model, class_index=top_channel, layer_name="Block8_6_Conv2d_1x1")
explainer.save(grid, ".", "grad_cam.png")

結果の例をいくつか貼っておきます.
abe.pngsuga.png

重要要素の決定の妥当性は置いておき,結果自体は顕著性っぽいものがとれてそうですね.主に鼻から目にかけての領域を一番重要視していると今回は判断できそうです.

自分が人の顔を見る際に真っ先に目が行くのが鼻や目のあたりなので,直感的にも正しい気がしています.

ちなみに,西洋人と東洋人で顔を観察するときの視点移動方法は違うらしいですね.今回の結果は鼻を中心に顕著性が分布したので,どちらかというと東洋人っぽい観察の特性だと思いました.
Culture Shapes How We Look at Faces

さいごに

今回はFaceNetが顔認識を行う際にどこを重視してみているのか,独自の指標で可視化してみました.この指標が本当に正しいかはちょっと担保できませんが,少し試してみる程度ならいい結果が得られたと思っています.

コードはこちらから利用できます.前回の記事のノートブックに追記する形で書いてあるので,後半まで順番に実行してください.

結局,修論のためにお休みをもらうといっておきながら記事のネタを思いついたので帰ってきてしまいました.また気が向いたら更新するくらいの気持ちでカレンダーを埋めていこうと思います.

参考

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

LeetCodeに毎日挑戦してみた 83. Remove Duplicates from Sorted List(Python、Go)

はじめに

無料英単語サイトE-tanを運営中の@ishishowです。

プログラマとしての能力を上げるために毎日leetcodeに取り組み、自分なりの解き方を挙げていきたいと思います。

Leetcodeとは

leetcode.com
ソフトウェア開発職のコーディング面接の練習といえばこれらしいです。
合計1500問以上のコーデイング問題が投稿されていて、実際の面接でも同じ問題が出されることは多いらしいとのことです。

golang入門+アルゴリズム脳の強化のためにgoとPythonで解いていこうと思います。(Pythonは弱弱だが経験あり)

19問目(問題83)

83. Remove Duplicates from Sorted List

問題内容

Given a sorted linked list, delete all duplicates such that each element appear only once.

(日本語訳)

ソートされたリンクリストを指定して、各要素が1*回*だけ表示されるように、すべての重複を削除します。

Example 1:

Input: 1->1->2
Output: 1->2

Example 2:

Input: 1->1->2->3->3
Output: 1->2->3

考え方

  1. headの参照をcurにコピーします。

  2. curをwhileループで回していき、cur.valとcur.next.valを比較して次のvalともし同じだったら次のノードを上書きします。

  3. 違う場合はそのままノードを進めていきます。

  4. 戻り値にはheadを返します。

  • 解答コード
class Solution(object):
    def deleteDuplicates(self, head):
        cur = head
        while cur and cur.next:
            if cur.next.val == cur.val:
                cur.next = cur.next.next
            else:
                cur = cur.next
        return head
  • Goでも書いてみます!
func deleteDuplicates(head *ListNode) *ListNode {
    cur := head
    for cur != nil && cur.Next != nil {
        if cur.Next.Val == cur.Val {
            cur.Next = cur.Next.Next
        } else {
            cur = cur.Next
        }
    }
    return head
}
  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

SCAPYでTCPを頑張ってみる(その1:準備編)

この文書の目的

この文書は SCAPY をインストールして TCP 通信を試してみよう、と思った人向けのものです。というのも、Ping (ICMP Echo) や DNS つまりコネクションレスなパケット操作については SCAPY はシンプルに、予想したとおりに書いて素直に処理できるのですが、TCP についてはいろいろと引っ掛かることが多いからです。
特にこの「準備編」では、TCP のために ACK などフラグを設定したパケットの操作で OS (カーネル)と衝突してしまい、コネクションをリセットされてしまう問題に焦点を当て、ARP Spoofing によってそれに対処する方法についてまとめました。
なお、主に仕組みを理解する過程に沿わせて書いていますので、TCP でデータを取得したいだけの人は、「その3:StreamSocketを使う」に直行するのが良いです。

前提

以下のことを前提に書いています。

  • SCAPY がインストールされている
  • TCP について一通りの理解がある・調べることができる

実験:3 Way Handshake 処理を書いてみる

以下に簡単に TCP の 3 Way Handshake 処理を書いてみました。適当な Web サーバに対して、TCP connection を開けてみます。以下のコードでは example.com としていますが、あなたがアクセスできる適当なホストに修正して実行すると良いです。

sport = random.randint(30000,60000)
seq = random.randint(1000,2000)
ip = IP(dst='example.com')
tcp = TCP(sport=sport, dport=80, flags='S', seq=seq)

syn = ip/tcp
syn_ack = sr1(syn)

tcp.seq = syn_ack.ack
tcp.ack = syn_ack.seq + 1
tcp.flags = 'A'
ack = ip/tcp
send(ack)

正常なケース

SCAPY に与えて実行すると、本来ならば以下のような表示になるはずです。赤文字箇所に注目してください。

Begin emission:
Finished sending 1 packets.
..............................................................................................*
Received 95 packets, got 1 answers, remaining 0 packets
.
Sent 1 packets.

つまり SYN 送信、SYN-ACK 受信、ACK 送信、の 3 パケットの動きが見えます。

このときの tcpdump によるモニタリングの結果は以下のようになるでしょう。Flags に続く [ ] 内の記号に注目してください。なお見やすさのために適度に不要な表示を除き、行番号がつけてあります。

 1: IP 192.168.1.2.44506 > example.com.http: Flags [S], seq 1080, length 0
 2: IP example.com.http > 192.168.1.2.44506: Flags [S.], seq 3401191913, ack 1081, length 0
 3: IP 192.168.1.2.44506 > example.com.http: Flags [.], ack 1, length 0

Flags の表示を見れば、正しく SYN, SYN-ACK, ACK の 3 パケットでハンドシェイクが行われていることがわかります。あなたの環境がこのように正しく 3 Way Handshake できるのであれば、この文書をこれ以降読む必要は無いかもしれません。次の「その2:はじめてのGET編」に進んでください。

うまくいかないケース

ところが全く同じコードを与えた私のもう一台の環境(MacOS 10.15.7)では異なる挙動がtcpdump で観察できます。次は RST が手元のマシンのOSから送られてしまう例です。

 1: IP 192.168.1.2.57947 > example.com.http: Flags [S], seq 1008, length 0
 2: IP example.com.http > 192.168.1.2.57947: Flags [S.], seq 3656582299, ack 1009, length 0
 3: IP 192.168.1.2.57947 > example.com.http: Flags [R], seq 1009, length 0
 4: IP 192.168.1.2.57947 > example.com.http: Flags [.], ack 1, length 0
 5: IP example.com.http > 192.168.1.2.57947: Flags [R], seq 3656582300, length 0

1, 2 パケットめまでは正常なケースでの SYN, SYN-ACK と同様ですが、3 パケットめは OS カーネルから RST フラグ [R] 付きパケットが送られています。4パケットめが一番目の事例の 3 パケットめと同じもの、つまりSCAPYのコードによって送られたACKです。
5パケットめは幾らか複雑で、つまり4パケットめによって相手ホスト example.com 側では通信は既に切断されたものとされ、そこに届いたACKに対して「いまさら」その相手はできない、とRSTで返されたことを示しています。

あなたの環境でもしこの [R] フラグが出るようであれば、これ以降を読んで対策すると良いです。

TCP reset 問題

この、何故か自分の手元の機械が TCP Connection をリセットしてしまう問題については、 ScapyによるTCP通信 に丁寧に解説されています。

なお、筆者は複数台の MacOS 10.15.7(19H2)マシンを使って実験していたが、このときRSTパケットを送るものと、そうでないものが発生し、何に原因があるか見いだせず混乱した状態に陥りました。(一時期、MacOS の Firewall の ON/OFF によって現象が分かれる状態になったため、これか、と思われたのですが、何度か試すうちに再現性が失われました。)

Workaround - 回避策

上に示した解説ではファイアウォールを用いてOSによる介入を迂回する対策が示されていますが、他にも解決方法はあります。たとえば(上の解説でも参照されている)こちらの What happens if you write a TCP stack in Python? 文書(日本語版)では ARP Spoofing を用いた迂回策が提示されています。

ここでは SCAPY による解決策として、ARP Spoofing 的な方法を試します。Firewallに手を入れるのは恒久的に作用して良いのですが、もともと SCAPY を使ってパケット処理を学ぼうとしているのですから、そのための教材としてこの作業は悪くありませんからね。

下図に元々の機器構成を示します。
structure.png

  • router の IP アドレスは 192.168.1.1
  • その先につながれている Mac のアドレスは 192.168.1.2
  • この Mac がインターネットのどこかにある example.com と通信する
  • router から送られてきた 192.168.1.2 宛てのパケットは、Mac の OS が受信する

これを以下のような構成だと router に信じさせるのです。
structure2.png

  • router の先には 192.168.1.180 が居る
  • それは Mac の中にある(っぽい)
  • この 192.168.1.180 がインターネットのどこかにある example.com と通信する
  • router から送られてきた 192.168.1.2 宛てのパケットは、Mac の OS が受信する
  • router から送られてきた 192.168.1.180 宛てのパケットは、SCAPY が受信する

つまり router は 192.168.1.180 宛てのパケットを Mac に接続されているネットワーク・インタフェイスに送りつけますが、Mac の OS は、これが「自分宛だ」とは思わず無視します(RSTパケットを送ることもない)。しかし Mac のインタフェイスを見張っている SCAPY はそのパケットを観測することができます。そして SCAPY の sr1( ) のような関数は「自分が送ったパケットへの返答だ」と考えて、これを「受信」処理に入れるのです。

手法の確認

先に乱暴なコードでの解決方法を示します。以下を SCAPY で実行することで、ターゲットのデバイス(router: 192.168.1.1) に、IP アドレス 192.168.1.180 が SCAPY プログラムが動作している手元の機器(Mac)のものだと教えることができます。

send(IP(dst='192.168.1.1', src='192.168.1.180')/ICMP())
time.sleep(0.3)
send(ARP(op='is-at', psrc='192.168.1.180', pdst='192.168.1.1', hwdst="ff:ff:ff:ff:ff:ff"))

動作としては、以下のようになることを期待しています。

  1. ターゲットに向けて ICMP echo request を出す
  2. ターゲットに ARP request を出させる
  3. ターゲットに向けて ARP response を送る
  4. すると、そこに向けて送られる ICMP echo reply が送り出される

ただし上記の乱暴なコードでは、4 番目のパケットについて受信処理を用意していません。投げっぱなしで、返信があっても捨てています。

以下にモニタリングの結果を示します。例によって見やすさのために適度に不要な表示を除き、行番号がつけてあります。

1: 3c:22:fb:03:55:c6 > c0:25:aa:2a:5e:42, IPv4, length 42: 192.168.1.180 > 192.168.1.1: ICMP echo request, id 0, seq 0, length 8
2: c0:25:aa:2a:5e:42 > ff:ff:ff:ff:ff:ff, ARP, length 60: Request who-has 192.168.1.180 tell 192.168.1.1, length 46
3: 3c:22:fb:03:55:c6 > c0:25:aa:2a:5e:42, ARP, length 42: Reply 192.168.1.180 is-at 3c:22:fb:03:55:c6, length 28
4: c0:25:aa:2a:5e:42 > 3c:22:fb:03:55:c6, IPv4, length 60: 192.168.1.1 > 192.168.1.180: ICMP echo reply, id 0, seq 0, length 8

行番号ごとに動きを説明します。先述の期待したとおりの結果になっています。

  1. SCAPY によって、source IP アドレスに 192.168.1.180 をつけた ICMP echo request が出る(ひょっとするとこれより前に router 向けの ARP 応答があるかもしれません)
  2. 1. に返信するために、router から 192.168.1.180 向けの ARP リクエストが出る
  3. (0.3sec 待ってから)SCAPY によって 192.168.1.180 はここ(自身の MAC アドレス)に居るぞ、と ARP response を合成して送り出す
  4. router はこの返答を信じて、ICMP echo reply を送出する

これが成功してからしばらく(ターゲットのARPキャッシュに情報が残っている間)は、単に ping を送っても返ってきます。これでターゲットの ARP テーブルに載っていることが確実ですね。(下の記述は send() ではなく sr1() を使って正しく受信処理をしています。)

>>> sr1(IP(dst='192.168.1.1', src='192.168.1.180')/ICMP())
Begin emission:
Finished sending 1 packets.
...*
Received 4 packets, got 1 answers, remaining 0 packets
<IP  version=4 ihl=5 tos=0x0 len=28 id=52956 flags= frag=0 ttl=64 proto=icmp chksum=0xbec src=192.168.1.1 dst=192.168.1.180 |<ICMP  type=echo-reply code=0 chksum=0xffff id=0x0 seq=0x0 |<Padding  load='\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00' |>>>
>>> 

コードの整理

以下にこれまでの内容を整理して、少しましにしたコードを示します。冒頭の、ターゲットに IP アドレスを教え込む処理に注目してください。実行すると、ICMP echo reply を受信したことが以下のメッセージで確認できます。
IP / ICMP 192.168.1.1 > 192.168.1.180 echo-reply 0 / Padding

これで OS によるリセット処理の介入なしに TCP connection を開けて、閉じることができると思います。

def arp_spoof(targetIP, injectIP):
    send(ARP(op='is-at', psrc=injectIP, pdst=targetIP, hwdst="ff:ff:ff:ff:ff:ff"))

fakeIP = '192.168.1.180'                # fake src IP address
gwaddr = conf.route.route('0.0.0.0')[2] # target gateway address

### inject IP address to target
t = threading.Timer(0.3, arp_spoof, args=(gwaddr, fakeIP, ))
t.start()
res = sr1(IP(dst=gwaddr, src=fakeIP)/ICMP())
print(res.summary())

### open TCP connection
sport = random.randint(30000,60000)
seq = random.randint(1000,2000)
ip = IP(dst='example.com', src=fakeIP)   # need to set injected IP
tcp = TCP(sport=sport,dport=80,flags='S',seq=seq)

syn = ip/tcp
syn_ack = sr1(syn)

tcp.seq = syn_ack.ack
tcp.ack = syn_ack.seq + 1
tcp.flags = 'A'
ack = ip/tcp
send(ack)

### INSERT YOUR ACTION HERE

### close TCP connection
tcp.flags = 'FA'
fin_ack = sr1(ip/tcp)

tcp.seq += 1
tcp.ack = fin_ack.seq + 1
tcp.flags = 'A'
send(ip/tcp)

この、TCP コネクションを開いて閉じるまでの間に、あなたのパケット処理を書けば良いのです。

「SCAPYでTCPを頑張ってみる(その2:はじめてのGET編)」に進んでください。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Markdownのテスト仕様書を管理したい!

忙しい方へ

おそらく、大抵のソフトウェアの仕様書はExcelで管理することが多いと思います。しかし、バージョン管理や差分がGitHubと比べてイマイチだと思っていました。

そのような問題を解決するために、まずは(個人的に)取り組みやすいテスト仕様書をMarkdown+GitHubで管理するという試みをしてみました。バージョン管理や差分管理は確かにやりやすくなりました。ついでに、一覧で眺めたり集計したいときに変換するpythonスクリプトも作成しました。

変換ソースコードは GitHub で公開しています。

image-20201130175253467.png

不要なExcelを撲滅したい

未だにソフトウェアの仕様書にExcelを管理する場合が多いと思います。私はこのExcel管理に嫌気がさしていました。

  • バイナリデータのため、GitHubなどのバージョン管理ツールとの相性が悪い。"〇〇_最終版.xlsx", "〇〇_最終版_new.xlsx", "〇〇_最終版(これを使う).xlsx"のようなファイルが同一ディレクトリ内にあると禿げそうになります。
  • 誰がいつどのように編集したのか、という差分管理もイマイチ。
  • ファイルが大きくなると、オープンや保存にも異様に時間がかかる。
  • そもそも仕様書がどこにあるかが、誰かに聞かないとわからない

等、挙げればキリがありません。ソースコードはGitHubで管理しているのだから、同じ場所・同じツールで仕様書を統一管理したいという思いが強くなっていきました。

とはいってもいきなり全部を置き換えるのは難しいため、まずは個人的に取り組みやすかったテスト仕様書Excelの撲滅に取り組みました。

Markdown+GitHubによるテスト仕様書管理

Markdown+GitHubならば上記の悩みをほぼ解決してくれましたが、Markdownはテキストベースのため、一覧表示ができなかったり、集計作業はできませんでした。そこで、Markdown型テスト仕様書をExcelに変換するpythonスクリプトを作成しました。ソースコードは GitHub で公開しています。

実際に使ってみて

実際のプロジェクトで運用してみて、良い点もあれば課題も見えてきました。

良かった点

  • ソースコードと同じディレクトリ階層で仕様書を管理できるため、仕様書を探す手間が省けた。
  • 仕様書内まで一括検索・置換ができるようになった(Excelは開くまで中身がわからなかった)。
  • バージョン管理や差分管理がものすごく快適になった。
  • オープンにもほとんど時間がかからず、生産性が向上した。
  • IDEだけで作業が完結するようになった。IDE/エディタとExcelのにらめっこ作業がなくなった。

課題

  • 図表現が弱い。テスト結果のスクリーンショットの貼り付けができない。
  • 集計作業ができない。ここはさすがにExcelに軍配。

さいごに

Markdownによる仕様書管理はまだまだ発展途上ではありますが、正直従来のExcel管理よりかは恩恵が大きいと感じました。旧式の方法を一新するのは中々骨が折れますが、できる範囲で続けていこうと思いました。次はもっと別の仕様書で。

参考文献

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Python+Redisでcredentialsとaccess_tokenの管理をする

この記事はJSL(日本システム技研) Advent Calendar 2020年 7日の記事です。
毎年のことながら遅刻してる… :innocent:

今回はPython+Redisの記事です。
自分はRDBMSばかり使っていて、NoSQL設計時にもどうも考え方の違いに混乱した身です。

この記事の目的

PythonとRedisで

  • credentialsの作成と管理
  • access_tokenの作成と管理
  • credentiasが作成したaccess_tokenの取得

をする。

DynamoDBを使えばセカンダリインデックスを使用してすぐできる問題なのですが、Redisを使うとなると設計周りで少し頭を使います。
意外にこの辺にターゲットを絞ってコミットした記事を見かけなかったため、同じように悩んでる人のために少しでも参考になれと思います:hand_splayed:

まずはRedisライブラリを使用してみる

前提

  • Redis Serverインストール済み
  • Pythonはvirtualenvにactivateしredis-pyをインストール済み

  1. Redis起動

    $ redis-server
    
  2. Pythonをインタラクティブモードで開く

    $ python
    
  3. clientを生成して操作する

    >>> from redis import Redis
    >>> client = Redis.from_url('redis://localhost:6379/1')
    >>> client.set('hoge', 'fuga')
    True
    >>> client.get('hoge')
    b'fuga'
    

こんな感じでRedis操作が可能です。
clientのメソッドはだいたいredis-cliと同じなので、もともとRedisを触っていた人ならすんなり入れるかと思います。

credentialsとaccess_tokenの関係

ココから本題です。
今回扱うcredentialsとaccess_tokenの構成は以下

credentials
key [primary]
secret
max_token_count

※credentialsのmax_token_countは、最大発行できるtokenの数が入ります。

access_token
token [primary]
expires_at

単純に作成するとこのようになります。(コピペしやすいよう>>>を排除しています)

import json
from redis import Redis
client = Redis.from_url('redis://localhost:6379/1')
credentials_params = {'secret': 'aaaaaa_secret', 'max_token_count': 3}
client.set('credentials:::aaaaaa_key', json.dumps(credentials_params))
> True
access_token_params = {'expires_at': '2021-05-08T09:46:10.363778+00:00'}
client.set('access_token:::aaaaaa_token', json.dumps(access_token_params))
> True
client.get('credentials:::aaaaaa_key')
> b'{"secret": "bbbbbbbb", "max_token_count": 3}'
client.get('access_token:::aaaaaa_token')
> b'{"expires_at": "2021-05-08T09:46:10.363778+00:00"}' 

同一DBを使用していることによるkeyの重複を防ぐため、set時にaccess_token:::credentials:::とprefixにつけています。
また、msetを使ってHash形式(辞書)で保存することもできるのですが、valueをjsonにしたほうがデータの加工がしやすいため、ここではsetを使っています。

問題点

簡単に図式化すると、今Redisは以下のような感じになっています。
redis01
視覚的にCredentialsとAccessTokenを分けていますが、実際にパーティションが区切られているということはありません。
Redis DBに対してcredentials:::aaaaaa_keyで取り出すかaccess_token:::aaaaaa_tokenで取り出すかの違いです。

ここで問題になるのはaccess_tokenをcredentials_keyで引っ張ってこれないため、max_token_countが役に立たないということです。
RDBMS脳な私は「access_tokenにcredentials_key列を追加して、そこから引っ張ればいいじゃん」と思っていたのですが、そこはNoSQL。そんなことはできません。
Redisへ格納できる形式を色々と調べてみましたが、不可能でした…

解決策

途方にくれる私に天啓が舞い降りた。
https://tylerstroud.com/2014/11/18/storing-and-querying-objects-in-redis/
< 「セット型使おうぜ」

...access_token生成/破棄でcredentialsに紐づくセットを増減させて手動でindexingさせればよいのでは!?

ということで対応

...
# set型でcredentialsに紐づく値=access_tokenを登録
client.sadd('credentials_has_access_tokens:::access_token:::aaaaaa_key', 'access_token:::aaaaaa_token')
> 1
# set型の検索
client.smembers('credentials_has_access_tokens:::access_token:::aaaaaa_key')
> {b'access_token:::aaaaaa_token'}
# もう一つトークン発行
client.set('access_token:::bbbbbb_token', json.dumps(access_token_params))
> True
client.get('access_token:::bbbbbb_token')
> b'{"expires_at": "2021-05-08T09:46:10.363778+00:00"}'
client.sadd('credentials_has_access_tokens:::access_token:::aaaaaa_key', 'access_token:::bbbbbb_token')
> 1
# 2つのaccess_tokenがcredentialsをキーに取得できる
client.smembers('credentials_has_access_tokens:::access_token:::aaaaaa_key')
> {b'access_token:::aaaaaa_token', b'access_token:::bbbbbb_token'}

これで当初の目的は達成できそうですね!
この対応で、下図CredentialsHasAccessTokensというcredentialsとaccess_tokenの集合ができました。
redis02.png

なお、今回の設計ではsmembersで引っ張ってくるのがあくまでaccess_tokenと定義しているため、expires_atなどの内部情報はcredentials_keyから取得することはできません。

Pythonコード化

上述はあくまでインタラクティブモードでの動作です。
ここから実際にPythonのコードに落とし込んでいきます。

import redis
import json
import random
import secrets
import string

from typing import Optional, Tuple
from datetime import datetime
from dateutil import relativedelta

# keyを引数にしてcredentialsを取得する(access_token発行時にcredentialsを確認する時等)
def get_credentials(client: redis.Redis, key: str) -> Optional[dict]:
    if not (credentials := client.get(f'credentials:::{key}')):
        return None
    return json.loads(credentials)

# access_tokenを引数にしてtokenに紐づく情報(=ユーザー情報)を取得する
def get_user(client: redis.Redis, access_token: str) -> Optional[dict]:
    if not (user := client.get(f'access_token:::{access_token}')):
        return None
    return json.loads(user)

# access_tokenを生成すると同時にkey - access_tokenの集合に追加し、access_tokenとexpires_atを返却
def create_access_token(client: redis.Redis, key: str, expires_days: int = 7) -> Tuple[str, str]:
    access_token = secrets.token_urlsafe(64)
    expires_at = (datetime.now() + relativedelta.relativedelta(days=+expires_days)).isoformat(timespec='seconds')
    client.set(
        f'access_token:::{access_token}',
        json.dumps({'key': key, 'expires_at': expires_at}),
    )
    client.sadd(f'credentials_has_access_tokens:::{key}', access_token)
    return access_token, expires_at

# 引数のcredentials_keyで生成したaccess_tokenの集合を返却
def issued_access_tokens(client: redis.Redis, key) -> set:
    return client.smembers(f'credentials_has_access_tokens:::{key}')

# credentialsを生成し、key/secretを返却
def create_credentials(client: redis.Redis) -> Tuple[str, str]:
    key = ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(32))
    params = {
        'secret':  secrets.token_urlsafe(64),
        'max_token_count': 3,
    }
    client.set(f'credentials:::{key}', json.dumps(params))
    return key, params['secret']

こんな感じになりました。
exampleなのでところどころ実用化する際には修正が必要になるかと思います。

まとめ

  • credentialsの作成と管理
  • access_tokenの作成と管理
  • credentiasが作成したaccess_tokenの取得

問題となったのは3つ目に関しては、集合を使うことで実現できましたが、今回の問題はPythonに限定されるわけではありません。
NoSQLは基本的に正規化の概念がないため、RDBMSで普通に行っていたことができず苦労するかもしれませんが、
考え方を変えることによっていくらでも実現は可能です。

redis-pyを使った良きPython+Redisがありますよう :pray:

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

SeleniumでハマったMaxRetryErrorの解決法

環境

  • MacOS : Catalina - 10.15.7
  • Python : Python 3.6.8
  • Selenium : 3.141.0

概要

業務で、競合他社のGoogle検索順位を調べるバッチ作成をしていたときのお話。
数個の「キーワード」をループでそれぞれ検索して、指定する他社のドメインが出てくるかをチェックするようなバッチだった。 以下、ソースコードに似たコード例を載せる。

sample.py
from selenium import webdriver
import chromedriver_binary
from selenium.webdriver.chrome.options import Options

# Chromeに接続して引数で検索する
def search(driver, word):
    driver.get("https://www.google.com")
    search = driver.find_element_by_name('q')
    search.send_keys(word)
    search.submit()
    return driver.page_source

def parse():
    ...省略...
    options = Options()
    options.add_argument('--headless')
    options.add_argument('--no-sandbox')
    options.add_argument('--disable-gpu')
    driver = webdriver.Chrome(options=options)

    for kw in keyword:
        # 検索欄にキーワードを入れて検索した後に現れる結果画面のソース
        source = search(driver, kw)
        result = analyze(source)
        ...省略...

キーワードが文字列のリストになっているから、一つずつ取り出して検索をかけて結果を得るって感じのコード。
しかし、実行すると、

エラー発生

urllib3.exceptions.MaxRetryError: HTTPConnectionPool(host='127.0.0.1', port=51691):
 Max retries exceeded with url: /session/fd963d7b05fb4b6344e0e8034c73ea67/url
(Caused by NewConnectionError('<urllib3.connection.HTTPConnection object at 0x7fc4e2773f60>:
Failed to establish a new connection: [Errno 111] Connection refused',))

上の例外が吐き出された。

試したこと

  • urllib3を最新にアップグレード
  • time.sleep(10)で検索の間隔を空けた

最初はurllib3のバージョンが悪いのかなと思ったけど違うっぽい

解決方法

結局、「設定したドライバーを検索ごとに再設定」したらいけた。
どういうことかというと、

sample.py
def parse():
    ...省略...
    options = Options()
    options.add_argument('--headless')
    options.add_argument('--no-sandbox')
    options.add_argument('--disable-gpu')
    driver = webdriver.Chrome(options=options)

    for kw in keyword:
        # 検索欄にキーワードを入れて検索した後に現れる結果画面のソース
        source = search(driver, kw)
        result = analyze(source)
        ...省略...

上のoptions = Options()からdriver = webdriver.Chrome(...)のところを、

sample.py
def parse():
    ...省略...
    for kw in keyword:
        # 検索欄にキーワードを入れて検索した後に現れる結果画面のソース
        options = Options()
        options.add_argument('--headless')
        options.add_argument('--no-sandbox')
        options.add_argument('--disable-gpu')
        driver = webdriver.Chrome(options=options)
        source = search(driver, kw)
        result = analyze(source)
        ...省略...

for文の中に入れて、検索するときはその都度ヘッドレスに設定することにした。

備考

  • Selenium のバージョン確認方法

Terminalから python3 を叩いて、下記を実行する。
この書き方じゃダメなモジュールもある。普通にpip3 list | grep "module name" の方が簡単かもしれない。

>>> import selenium
>>> selenium.__version__
'3.141.0'
  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

seabornで複数のグラフを重ねてプロットする

seabornのcatplotを使った時に複数のグラフを重ねるやり方がすぐ分からなかったのでメモしておく
ついでにcatplotを使わない場合もメモ

実行環境

jupyter notebook
python 3.8
seaborn version 0.11.0

catplotを使わずにグラフを重ねる方法

まずcatplotを使わずにシンプルに重ねる方法

sns.boxplot(data=df, x=metrics, y='Group', orient='h', color='white')
sns.swarmplot(data=df, x=metrics, y='Group', orient='h', hue='Classification', size=pointsize, palette='Set2')

これはそのままグラフを描くコードを別々に書けば良い.

jupyterではこれで複数のグラフを重ねられる.(グラフの形状が似たようなものの場合)

今回はboxplotswarmplotを重ねている.

jupyterでグラフを書き,plt.savefig('hoge.png')で保存すると以下のようになる.

casp14_gdtts_loss_group_boxplot.png

グラフを重ねることには成功.

しかしjupyterでは正しくグラフが表示されているが,そのまま保存するとlabelが見切れてしまった状態の画像が保存される.

これを直すためには

sns.boxplot(data=df, x=metrics, y='Group', orient='h', color='white')
sns.swarmplot(data=df, x=metrics, y='Group', orient='h', hue='Classification', size=pointsize, palette='Set2')
plt.tight_layout()
plt.savefig('hoge.png')

のように保存する前にplt.tight_layoutを使うか,

sns.boxplot(data=df, x=metrics, y='Group', orient='h', color='white')
sns.swarmplot(data=df, x=metrics, y='Group', orient='h', hue='Classification', size=pointsize, palette='Set2')
plt.savefig('hoge.png', bbox_inches='tight')

のように保存する時にbbox_inches='tight'を渡すと以下のようにグラフ全体を保存することができる.

casp14_gdtts_loss_group_boxplot_tight.png

このように見切れていないグラフの画像を保存することができた.

catplotでグラフを重ねる方法

catplotを使って画像を保存したい場合,以下のようにcatplotを2つ書くと2つのグラフが生成されてしまう.

sns.catplot(data=df, x=metrics, y='Group', orient='h', hue='Classification', size=8, kind='swarm')
sns.catplot(data=df, x=metrics, y='Group', orient='h', color='white', kind='box')

catplotでグラフを重ねたいときは以下のようにする.

g = sns.catplot(data=df, x=metrics, y='Group', orient='h', hue='Classification', size=8, kind='swarm')
g.map(sns.boxplot, data=df, x=metrics, y='Group', orient='h', color='white')

まず1行目は普通にcatplotで良い.

この時に返り値を受け取っておく.(ここではgとしている.seaborn.FaceGridが返ってくる)

続いて返り値(seaborn.FaceGrid)に重ねたいグラフをmapする.

例えばboxplotを重ねたいとき

sns.boxplot(data=df, x=metrics, y='Group', orient='h', color='white')

普通に書く場合は上記の書き方となるが,catplotに重ねる場合は

g.map(sns.boxplot, data=df, x=metrics, y='Group', orient='h', color='white')

のように第一引数に重ねたいグラフの名前(ここではsns.boxplot),それ以降の引数にはグラフを書く時に引数として与えたい引数(ここではdata=df, x=metrics, y='Group', orient='h', color='white')を渡す.

このようにするとcatplotを重ねることができた.

catplotで作成した画像は以下のようになった.

casp14_gdtts_loss_group_boxplot_catplot.png

参考

seaborn公式ドキュメント
https://seaborn.pydata.org/generated/seaborn.FacetGrid.html#seaborn.FacetGrid

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

簡単なFlaskの始め方【Python】

この記事では、
簡単なFlaskの始め方
について書いていきます。

インストールからWEBに反映させるところまでやります。

Flaskとは

Flask(フラスク)は、プログラミング言語Python用の、軽量なウェブアプリケーションフレームワークである。標準で提供する機能を最小限に保っているため、自身を「マイクロフレームワーク」と呼んでいる。
Flask - Wikipedia

やり方

まずはインストールをしていきます。

インストール

ターミナルで下のコードを実行してください。

$ pip install flask

ディレクトリ作成

次はディレクトリを作成します。
下のコードもターミナルで実行してください。

$ cd Desktop
$ mkdir flask_study

flask_studyの部分はなんでも良いです。
これを実行すると、デスクトップに新しいフォルダが作成されます。
それだけの事なので、普通に作っても変わらないと思います。

次は、今作ったディレクトリに移動します。

$ cd flask_study

移動したら、下のコードを実行してください。
名前の横に(venv)が出てきたらできています。

$ python3 -m venv venv
$ . venv/bin/activate            

Pythonを書く

作ったディレクトリ内でPythonを書いてください。

app.py
from flask import Flask

app = Flask(__name__)

@app.route("/")
def index():
    return "Hello World!"

if __name__ == "__main__":
    app.run(debug=True)

ターミナルでapp.pyを実行します。

$ python app.py

実行結果

...(省略)
 * Debug mode: on
 * Running on http://127.0.0.1:5000/ (Press CTRL+C to quit)
...(省略)

ここに記載されている http://127.0.0.1:5000/ にアクセスしてみましょう。
スクリーンショット 2020-12-08 18.39.20.png
このように、Hello World! が表示されました。

まとめ

今回は、
簡単なFlaskの始め方
について書いていきました。

次回は、HTMLを反映させる方法について書いていきたいと思います。

ありがとうございました。

参考

Flask - Wikipedia

注意

この記事は、プログラミング初学者が書いており、内容に誤りがある場合がございます。
あしからずご了承願います。
また、誤りにお気づきの場合にはご指摘いただけると幸いです。
よろしくお願いいたします。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

局所線形回帰を完全に理解したい

この記事は九工大古川研 Advent Calendar 2020 8日目の記事です.
本記事は古川研究室の学生が学習の一環として書いたものです.
内容が曖昧であったり表現が多少異なったりする場合があります.

まえがき

局所線形回帰を完全に理解した (「い」が大事)というモチベーションで
局所線形回帰について勉強してまとめました.

研究室に入りたての私でも理解できるような、
ちょっとクドいくらい丁寧に説明できたらと思っております。
(サークルの後輩も読んでくれるそうなので)

今回は線形回帰の拡張として解説します.

問題設定

${(x_i, y_i)}_{i=1}^{n}$ が与えられているとします.($ x, y$ のペアの集合,$i = 1, \cdots, n$)

ここで,$x$ は入力, $y$ は出力に相当するものです.

与えられた ${(x_i, y_i)}_{i=1}^{n}$ をもとに,
新規の入力 $x^{\ast}$ に対する出力値を予測したいとします.

qiita_pictures (1)-02.png

このときに,$y = f(x)$ となるモデル $f$ をデータから学習しましょうという話です.

線形回帰とは

いわゆる最小二乗法です.
「そんなもんわかっとるわい」という方は
読み飛ばしてもらって大丈夫です.

線形回帰では,$ f $ を一次関数として問題を解きます.
具体的には,
$ f(x) = a_0 + a_1 x $ などとおいて,データから $a_0$ と $a_1$ を推定します..
$a_0$ が切片(0次のパラメータ),$a_1$ が傾き(1次のパラメータ)です.

qiita_pictures (1)-04.png

このモデル $f$ を用いて推定した出力を $\hat{y}$ として,
直線の式を以下のように書きます.

\begin{align}
\hat{y} &= a_1 x + a_0 \\
&= \left( \begin{array}{cc} 1 & x \\ \end{array} \right)
\left( \begin{array}{c} a_0 \\ a_1 \end{array} \right) \\
\end{align}

ここで,

\begin{align} 
\boldsymbol{x} := \left( \begin{array}{c} 1 \\ x \end{array} \right), \,\,\, 
\boldsymbol{a} := \left( \begin{array}{c} a_0 \\ a_1 \end{array} \right)
\end{align}

として,以下のように表します.

y = \boldsymbol{x}^T \boldsymbol{a}

さて,パラメータ $\boldsymbol{a}$ をどうやって決めるのかという話ですが,
ある入力 $x_i$ に対する「出力の推定値 $\hat{y_i}$」 と「実際の出力 $y_i$」 の二乗誤差を
最小にするように $\boldsymbol{a}$ を決めます.
まあこの辺りは調べればいくらでもヒットするような内容ですので多くは語りません.

qiita_pictures (1)-06.png

『ある入力 $x_i$ に対する「出力の推定値 $\hat{y_i}$」 と「実際の出力 $y_i$」 の二乗誤差』の総和を $E$ とします.
いわゆる目的関数です.

\begin{align}
E &= \sum_i^n ( y_i - \hat{y}_i )^2 \\
&= \sum_i^n ( y_i - \boldsymbol{x}_i^T\boldsymbol{a} )^2 \\
\end{align}

$E$ を最小にするようなパラメータ $\boldsymbol{a}$ を求めていきます.

ここで,敢えて $\sum$ を外した形で書いてみます.

E =  ( y_1 - \boldsymbol{x}_1^T\boldsymbol{a} )^2 + \cdots + ( y_i - \boldsymbol{x}_i^T\boldsymbol{a} )^2

そしてこれを行列積の形で表現してみると以下のようになります.

E = \left( \begin{array}{cccc} y_1 - \boldsymbol{x}_1^T\boldsymbol{a} & \cdots & y_n - \boldsymbol{x}_n^T\boldsymbol{a}  \end{array} \right)
\left( \begin{array}{c} y_1 - \boldsymbol{x}_1^T\boldsymbol{a} \\ \vdots \\ y_n - \boldsymbol{x}_n^T\boldsymbol{a} \\ \end{array} \right) \\

ここで,

\boldsymbol{y} := \left( \begin{array}{c} y_1 \\ \vdots \\ y_n \\ \end{array} \right), \,\,\, \boldsymbol{X} := \left( \begin{array}{c} \boldsymbol{x}_1^T \\ \vdots \\ \boldsymbol{x}_n^T \\ \end{array} \right)

とすることで,以下のような変形ができます.

\left( \begin{array}{c} y_1 - \boldsymbol{x}_1^T\boldsymbol{a} \\ \vdots \\ y_n - \boldsymbol{x}_n^T\boldsymbol{a} \\ \end{array} \right) = \left( \begin{array}{c} y_1 \\ \vdots \\ y_n \\ \end{array} \right) - \left( \begin{array}{c} \boldsymbol{x}_1^T \\ \vdots \\ \boldsymbol{x}_n^T \\ \end{array} \right) \boldsymbol{a} = \boldsymbol{y} - \boldsymbol{X} \boldsymbol{a}

すなわち先程の $E$ 式は以下のように書けます.

\begin{align} 
E &= \left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right)^T
\left(\boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right) \\
&= \left\|\boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right\|^2
\end{align}

$\boldsymbol{y} - \boldsymbol{X} \boldsymbol{a}$ のユークリッドノルムの二乗の形ですね.

上で変形した $E$ の式について,
右辺を展開してから $\boldsymbol{a}$ で微分します.

\begin{align} 
\frac{\partial}{\partial \boldsymbol{a}} E
&= \frac{\partial}{\partial \boldsymbol{a}}\left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right)^T
\left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right) \\
&=  \frac{\partial}{\partial \boldsymbol{a}}
\left( \boldsymbol{y}^T - \boldsymbol{a}^T\boldsymbol{X}^T \right)
\left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right) \\
&= \frac{\partial}{\partial \boldsymbol{a}}
\left( \boldsymbol{y}^T\boldsymbol{y}
- \boldsymbol{y}^T\boldsymbol{X}\boldsymbol{a}
- \boldsymbol{a}^T\boldsymbol{X}^T\boldsymbol{y}
+ \left(\boldsymbol{X} \boldsymbol{a} \right)^T \left(\boldsymbol{X} \boldsymbol{a} \right)
\right) \\
&= 0 - \boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{X}^T \boldsymbol{y} + 2 \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{a} \\
&= - 2 \boldsymbol{X}^T \boldsymbol{y} - 2 \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{a}
\end{align}

となります.
そして,これを最小とする $\boldsymbol{a}$ ,すなわち

\frac{\partial}{\partial \boldsymbol{a}} E = 0

となる $\boldsymbol{a}$ を求めますと,

\begin{align}
- 2 \boldsymbol{X}^T \boldsymbol{y} - 2 \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{a} &= 0 \\
2\boldsymbol{X}^T \boldsymbol{X} \boldsymbol{a} &= 2 \boldsymbol{X}^T \boldsymbol{y} \\
\boldsymbol{a} &= \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1}\boldsymbol{X}^T \boldsymbol{y}
\end{align}

となります.
めでたしめでたし.

また,ここで求まった $\boldsymbol{a}$ は
ムーア・ペンローズの一般化逆行列を用いて以下のようにも表現できます.

\boldsymbol{a} = \boldsymbol{X}^{\dagger} \boldsymbol{y}

ただし,$\boldsymbol{X}^{\dagger}$ は以下のように定義されます.

\boldsymbol{X}^{\dagger} := \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1}\boldsymbol{X}^T

局所線形回帰に拡張

ここからが本題です.

上で説明した線形回帰では
推定した出力 $\hat{\boldsymbol{y}_i}$ と実際の出力 $\boldsymbol{y}_i$ の誤差を
すべてのデータに対して平等に重み付けしてパラメータを計算していました.

局所線形回帰では,
新規の入力に近いデータに対しては誤差を重く
逆に遠いデータに対しては誤差を軽く見て
パラメータを計算します.

イメージは以下の図の通りです.

qiita_pictures (1)-08.png

qiita_pictures (1)-09.png

qiita_pictures (1)-10.png

$x^\ast$ の値によってパラメータが異なるのも
なんとなくわかるかと思います.

このときの誤差の重みを決める関数(カーネル関数)$ k(x, x^\prime) $ は
以下の式で表されるとします.

k(x, x^\prime) = \exp \left(- \frac{1}{2  \sigma^2} ( x - x^\prime )^2 \right)

$\sigma$ はカーネル幅と呼ばれる(ハイパー)パラメータで,
ざっくり言いますと
「新規データに対する近い or 遠いを決める境目」を決めるものです.
これを大きくするとより広い範囲のデータを近傍と見なします.
これを究極に大きくすると先程の線形回帰と同じ結果を得られます.

\begin{align}
E &= \sum_i^n k(x^\ast, x_i) ( y_i - \boldsymbol{x}_i^T\boldsymbol{a} )^2 \\
&= \left( \begin{array}{ccc} y_1 - \boldsymbol{x}_1^T\boldsymbol{a} & \cdots & y_n - \boldsymbol{x}_n^T\boldsymbol{a} \end{array} \right)
\left( \begin{array}{ccc}
k(x^\ast, x_1) & & \huge{0} \\
 & \ddots  & \\
\huge{0} & & k(x^\ast, x_n) \\
\end{array} \right)
\left( \begin{array}{c} y_1 - \boldsymbol{x}_1^T\boldsymbol{a} \\ \vdots \\ y_n - \boldsymbol{x}_n^T\boldsymbol{a} \\ \end{array} \right) \\
\end{align}

ここで

\begin{align} 
\boldsymbol{y} &:= \left( \begin{array}{c} y_1 \\ \vdots \\ y_n \\ \end{array} \right)\\[10pt]
\boldsymbol{X} &:= \left( \begin{array}{c} \boldsymbol{x}_1^T \\ \vdots \\ \boldsymbol{x}_n^T \\ \end{array} \right) \\[10pt]

\boldsymbol{W}_{x^\ast} &:= \left( \begin{array}{ccc}
k(x^\ast, x_1) & & \huge{0} \\
 & \ddots  & \\
\huge{0} & & k(x^\ast, x_n) \\
\end{array} \right)
\end{align}

とおくと以下のように書けます.

E = \left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right)^T
\boldsymbol{W}_{x^\ast} \left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right)

線形回帰の時と似たような式が出てきましたね.
実際,この式は「$\boldsymbol{W}_{x^\ast}$ をメトリックとするノルム」として

E = \left\| \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a}  \right\|_{\boldsymbol{W}_{x^\ast}}^2

とも表現できる(らしい)です.(要出典)

そして,線形回帰の時と同様に
$\boldsymbol{a}$ で偏微分すると以下のようになります.

\begin{align}
\frac{\partial}{\partial \boldsymbol{a}} E &= \frac{\partial}{\partial \boldsymbol{a}} \left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right)^T
\boldsymbol{W}_{x^\ast} \left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right)\\
&= \frac{\partial}{\partial \boldsymbol{a}} \left( \boldsymbol{y}^T - \boldsymbol{a}^T\boldsymbol{X}^T \right) 
\boldsymbol{W}_{x^\ast} \left( \boldsymbol{y} - \boldsymbol{X}\boldsymbol{a} \right)\\
&=  \frac{\partial}{\partial \boldsymbol{a}}\left(
\boldsymbol{y}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y}
- \boldsymbol{y}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X}\boldsymbol{a}
- \boldsymbol{a}^T\boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y} 
+ (\boldsymbol{X}\boldsymbol{a})^T \boldsymbol{W}_{x^\ast} \boldsymbol{X}\boldsymbol{a} \right)\\
&= 0 - \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y} - \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y}  + 2 \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X} \boldsymbol{a} \\
&= - 2 \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y} + 2 \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X} \boldsymbol{a} \\
2 \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X} \boldsymbol{a} &= 2 \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y} \\
\boldsymbol{a} &= \left(  \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y}
\end{align}

結構大変ですね.

目的関数 $E$ をパラメータ $\boldsymbol{a}$ で微分した式が求まりましたので,

\frac{\partial}{\partial \boldsymbol{a}} E = 0

となるような $\boldsymbol{a}$ を求めていきます.

\begin{align}
2 \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X} \boldsymbol{a} &= 2 \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y} \\
\boldsymbol{a} &= \left(  \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{y}
\end{align}

ということで,目的関数 $E$ を最小とする
パラメータ $\boldsymbol{a}$ が求められました.
この式も線形回帰の時に求めたものと似ていますね.

「$\boldsymbol{W}_{x^\ast}$ をメトリックとする $\boldsymbol{X}$ の一般化逆行列 」を

\boldsymbol{X}_{\boldsymbol{W}_{x^\ast}}^{\dagger} :=  \left(  \boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \boldsymbol{W}_{x^\ast}

のように定義することで

\boldsymbol{a} =\boldsymbol{X}_{\boldsymbol{W}_{x^\ast}}^{\dagger}  \boldsymbol{y}

とも書けます.

実際に回帰してみる

Python で実際に回帰させてみます.
コードは最後に載せます.

学習データは $y = sin(x) + \epsilon$ ($\epsilon$ はガウシアンノイズ)としました.

data.png

線形回帰

まず線形回帰で推定した場合です.

lr (1).png

局所線形回帰

局所線形回帰で推定した場合です.
カーネル幅 $\sigma$ を $\sigma = 0.5, 1.0, 2.0, 10.0$ と変えて比較もしてみました.
($\sigma$ の値は図中に記しております.)

llr_05.png
llr_1.png
llr_2.png
llr_10.png

カーネル幅 $\sigma$ が小さいほど,より狭い領域で直線を推定します.
一方で,カーネル幅 $\sigma$ を極端に大きくすると
線形回帰と(ほぼ)同じ推定結果を得られます.

おまけ

局所線形回帰では,ある $\boldsymbol{x}$ に対する出力 $\boldsymbol{\hat{y}}$について,
局所的に一次関数で近似します.
ただ,一次関数にこだわる必要はなくて,
二次関数でも三次関数でもゼロ次(定数)関数でも良いわけです.
このように,局所線形回帰に対し
「局所的に n 次多項式で回帰する」というように
一般化した回帰モデルを「局所多項式回帰」と呼びます.

ちなみに,Nadaraya-Watson カーネル回帰は
局所多項式回帰のゼロ次関数版です.
Nadaraya-Watson is 何?という話は
後日 @nakashima1125 氏がわかりやすく解説してくれると思います.

付録:実装まとめ

Jupyter Notebook のままの方が見やすいかなと思いましたので
Gist のリンクを貼ることにしました.

local_linear_regression.ipynb

良かったら Google Colab で開いて動かしてみてください.
カーネル幅(最後のセルの sigma)を色々いじって遊んでみると楽しいかも知れません.
カーネル幅を小さくしすぎると
$\boldsymbol{X}^T \boldsymbol{W}_{x^\ast} \boldsymbol{X}$ の逆行列が計算できなくなってコケるので注意です.


ここまで読んでくださってありがとうございました.
編集リクエストもお待ちしてます().

参考文献

  • The Matrix Cookbook
  • Trevor Hastie et al.(2014)「統計的学習の基礎 ―データマイニング・推論・予測―」共立出版
  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Anacondaが有償化されて困っている人に贈る、Pythonのパッケージ管理

Anacondaの商用有償化

こちらの記事でご存知の方も多いかと思いますが、
2020/4より、従業員数200名以上の営利企業では、Anacondaの利用が有償化されたようです
(無料のIndividual Editionではなく、Commercial Editionを利用するよう規約改定)

この改正により、
業務でAnacondaを利用中で代替手段に困っている方
が巷にあふれているかと思います! (私もその一人です‥泣)

何とかならないものかと調査し、代わりとなる手段の比較と導入法をまとめてみました。

私が独自に調べてまとめた記事ですので、間違い等ございましたら指摘頂けるとありがたいです!

代替手段まとめ

いきなり結論ですが、WindowsにおけるAnaconda無償版の代替ツールを表にまとめてみました。
大きく下記3通りの代替手段がありそうです。

・各ツールのメリットまとめ

venv+pip pipenv Anaconda有償版
導入・学習の手軽さ 〇 (インストール不要) × (pipfileの学習が必要) × (一見手軽だが実は複雑)
ライブラリ依存関係の解決 ×
行列計算速度 × ×
コスト 〇 (無料) 〇 (無料) × (14.95$/月)



全てを読むのが面倒くさい!
という方は、下記フローを参考に実際の使用法まで飛んでください

1. 手軽に環境切替したい or 仮想環境自体が不要な方
Python組み込みのvenv+pipでOK (リンク先に使用法)

2. チームで環境を共有したい、ライブラリの依存関係も管理したい、でもタダがいい方
pipenvを使用 (リンク先に使用法)

3. お金($14.95/月)を払ってでも機械学習や数値計算の処理速度を求める方
Anaconda Commercial Editionを使用

仮想環境とは?

上記選択を判断するうえで、「仮想環境って何やねん?」と思われる方もいるかもしれないので、解説します。

まず、ここでいう「環境」とは、「Python本体+ライブラリ」の事を指します。
WindowsにPythonやそのライブラリをインストールすると、これが「デフォルト環境」になるわけですが、
この「デフォルト環境」とは別に、プロジェクトごとに環境を作る事を「仮想環境」といいます

こちらに記載されているように、トラブル特定や管理の容易さからインストールされるライブラリは必要最低限が望ましく、これを実現するためにはプロジェクトごとに仮想環境を持つ事があるべき姿と言えます。
一方で仮想環境の管理には手間も学習コストも掛かり、個人のちょっとした開発では割に合わないとも言えるので、結局のところ

仮想環境: チーム開発や、リリースに近い開発向け
デフォルト環境をそのまま使用: 個人での分析向け
といった感じで使い分けるのが良いのではと、個人的には思います。

「ようわからん!」という方へのおすすめは、まずは「1.venv+pip」の仮想環境不要の方法で開発環境を作り、チームで仮想環境を共有したくなったら「2.pipenv」に移る、という方法です



以下で各代替法の詳細な比較をします
※「パッケージ」と「ライブラリ」で表記ゆれがありますが、同じ意味と捉えて頂いて構いません

何が困る?

Anacondaの機能には、大きく分けて下記2種類があります。

①仮想環境の管理 : Pythonやライブラリのバージョンをプロジェクトに合わせ切り替える
②パッケージインストール : ライブラリをインストールする (例:conda install pandas)

どちらもPythonを使用する上で必須ともいえる機能 (特に後者)なので、代替手段を準備する必要があります。

また、上記以外にもAnacondaは③行列やFFTの計算速度向上という機能を持っています。
AIや数値計算にPythonを利用されている方には重要機能だと思うので、こちらについても詳説します。

対処のための前提条件

まず最初に、この記事を開いた時点で「従業員数200名以上の企業」で「Anacondaを業務使用」している方が多いかと思います。

明確な根拠があるわけではありませんが (参考)、規模の大きな企業ではWindowsの普及率が極めて高いと思われ、またAnaconda自体もWindowsでのシェアが高いツールなので、

Windows環境
を前提として記事を書かせて頂きます。

また、エディタとしてVisual Studio Code (VSCode)を使用している前提で記事を書きますが、
Pycharm等他のエディタでもほぼ同様の手順で対処できるかと思います。

どう対処する?

こちらの記事でまとめられているように、Pythonの主要パッケージ管理ツールを一覧表記します

※なお、前述「①仮想環境の管理」は、「インタプリタ切替」「パッケージ切替」の2つに分けております
・インタプリタ切替:Python自体のバージョン切替(例:Python3.8.6)
・パッケージ切替:ライブラリのバージョン切替(例:Pandas1.1.4)

①インタプリタ切替 ①パッケージ切替 ②パッケージインストール ③行列計算高速化
pyenv
virtualenv
pip
venv
virtualenvwrapper
pipenv
Anaconda

選択肢が多く見えますが、見ての通りAnacondaの機能をカバーするためには複数ツールを組み合わせる必要があり、
私が調べた限り (こちらが非常に参考になりました)、Windows環境において
Anaconda無償版の代替には、大きく下記3つの選択肢がありそうです

1. venv + pip
2. pipenv
3. Anaconda有償版

※pyenvが含まれていない事が気になる方もいるかもしれませんが、
こちらの記事のフローチャートを見る限り、pyenvが必要となるのは
「Linux or Mac環境」かつ「インタプリタ切替が必要となる」
場合なので、Windows環境を前提とした本記事では、pyenvに関しての言及は省略します

※Dockerとpipを併用する方法もありますが、Linux、Docker、VSCodeからのアクセス等学習コストが他の方法と比べ非常に大きいため、本記事のコンセプトから外れると判断して割愛します (希望が多ければ記事にしようと思います)

各手法の比較

トップに記載したように、上記3種類の代替手段のメリットデメリットをまとめました

venv+pip pipenv Anaconda有償版
A.導入・学習の手軽さ 〇 (インストール不要) × (pipfileの学習が必要) × (一見手軽だが実は複雑)
B.ライブラリ依存関係の解決 ×
C.行列計算速度 × ×
D.コスト 〇 (無料) 〇 (無料) × ($14.95/月)

詳細を下記します

A. 導入・学習の手軽さ

導入や学習コストの低さに関して

venv+pip

全てPythonにデフォルトで組み込まれた機能のため、新たなインストールが不要です
こちらに記載されているように、Python3.3以降はvirtualenvもvenvに取り込まれたようなので、pipと組み合わせてひととおりの機能を実現できるようになりました。

基本的な動作はシェルコマンドのみで実現できるため、学習コストが低いです。
使用するたびに仮想環境切替のコマンドを打つ必要があるのが手間ですが、後述のようにプロジェクトフォルダ直下に環境を作れば、VSCodeが勝手に読み込んでくれます。

pipenv

インストールが必要ですが、pip install pipenvの1コマンドで完了するので手間ではありません。
仮想環境の作成も簡単なコマンドで実現できますが、
pipenvコマンドの種類がやや多く、ある程度の学習コストが必要なことや、
パッケージ本体がどこにあるかが分かりにくい (こちらに記載されているようにC:\Users\[ユーザ名]\.virtualenvs\にあるようです)ことから、
venv+pipほどのシンプルさはないように感じました

Anaconda

Anacondaはよく使うパッケージを一括でインストールしてくれ、Navigator等のGUI管理ツールも充実しており一見使いやすく見えます。

ただし、ライブラリの依存関係管理がブラックボックス気味であったり他のツールのパスを隠蔽してしまったり、pipと併用するとパッケージが壊れることがある、という問題があるようです。

私もついやってしまいがちですが、Qiitaの記事を見ながらとっさにpip installと打ってしまってpipとの併用状態となり、パッケージが壊れて1からインストールし直し、という事もありました。

上記のような現象を完全にコントロールするためには相当な知識が必要とされるため、いざトラブルが起きた時の対処の難しさを考えると、「学習コストが低い」とは言い難いのではと思います。

B.依存関係の管理

ライブラリ間の依存関係を管理する機能についてです。

venv-pip

後述するように、インストールしたライブラリの一覧をrequirements.txtというファイルでエクスポートすることができます。
ただ、requirements.txtにはライブラリのバージョンしか記載されておらず、依存関係の記載が不完全ですし、requirements.txtの出力自体も手動で行う必要があります。
というわけで、依存関係の管理機能に関しては、下記のpipenvやAnacondaの方が充実していると言えそうです。

pipenv

使用法は後述しますが、pipenvの仮想環境ではこちらに記載されているように、「Pipfile」と「Pipfile.lock」という2つのファイルで依存関係を管理できます。
PipfileやPipfile.lockはライブラリインストール時に自動で更新されるため、requirement.txtのように出力し忘れて更新が反映されない心配もありません。

Anaconda

こちらに記載されているように、Anaconda社のリポジトリ管理の一環としてライブラリの依存関係もチェックしているため (ややブラックボックスではありますが)、ある程度の依存関係管理が保証されていると思われます。

また、conda info [パッケージ名]で依存関係を表示させることもできます。

C.行列計算速度

こちらに実際の比較結果が記載されている通り、Anacondaは公式のPythonと比べて行列の計算速度が高くなる機能が実装されています。

具体的には
・NumPy
・NumExpr
・SciPy
・Scikit-Learn
・Tensorflow

などのライブラリ (どれも速度が求められるツールですね‥)において、Anacondaバージョンのみ
Intel MKL
と呼ばれる数値計算を高速化するためのライブラリが採用されているそうです。
(ちなみに、PytorchはAnacondaのみならずpipでも、Intel MKL採用バージョンがインストールされるようです)

これらの数値計算・機械学習系ライブラリで速度を求める場合、有償版Anacondaの採用を検討しても良いかもしれません

実際どれくらい早くなるか気になるところではあるので、よく使う機械学習処理で比較して別途記事を作成しようと思います。

【追記】
手間は掛かりますが、どうやらpipでもIntel MKLバージョンのNumpy/SciPyが使えるようです。調べた限りはLinuxの記事ばかりでWindowsでも実現できるかは分かりませんが、上記速度比較記事を書く際に試してみようと思います。

D.コスト

これは見ての通り、Anaconda有償版が$14.95掛かるのに対し、venv + pipやpipenvは無償で利用できます。


ここまで3種類の代替手段を比較しましたが、ここからは実際の使用法を記載します

1. venv + pipの使用法

VSCodeのシェルを使って作業する前提で、使用手順を解説します。

【参考にさせて頂いた記事】
https://qiita.com/fiftystorm36/items/b2fd47cf32c7694adc2e
https://qiita.com/po3rin/items/5c853bc172e5de330148
https://www.python.jp/install/windows/venv.html

Anacondaのアンインストール

こちらの「方法2」で、完全にアンインストールしてください

PythonとVSCodeのインストール

こちらを参考にPythonとVisual Studio Codeをインストールしてください
PyCharm等を使用したい場合は、VSCodeのインストールは飛ばして頂いて構いません

Pythonインストール時は、下図赤枠Add Python3.* to Pathをチェックすると、デフォルトでのPython環境に設定できます
image.png

また、VSCodeインストール時は下記2箇所をチェックするとエクスプローラの右クリックで開けて便利です
image.png

ここからのフローですが、仮想環境を使用するかしないかでフローが分かれます

・仮想環境を使用し、プロジェクトによりライブラリのバージョンを分けたい方
そのまま読み進めてください
(Pythonのバージョンも環境ごとに変えたい場合はこちら参照)

・仮想環境は使用せず、デフォルト環境でパッケージインストールできれば良い方
ライブラリのインストール(仮想環境不使用時)に飛んでください
(ようわからん!という方はこちら推奨)

1-1. 仮想環境使用時の手順

下記のように
仮想環境の作成 → アクティベート → ライブラリのインストール → VSCodeの設定
の順で進めてください

venv仮想環境の作成

下図のような感じで環境構築したいディレクトリ(プロジェクトディレクトリ)を準備しておきます
image.png

ディレクトリ内に新環境を作成します

> cd [環境構築したいディレクトリ]
> python -m venv [新しい環境名]

下図のように環境が生成されます(下例では環境名「newvenv」で作成)
image.png

中を見るとexeファイルもろもろが生成されていることが分かります
image.png
なお、環境の管理方法としては、
プロジェクトディレクトリ(≒VSCodeのようなエディタで開くディレクトリ。上例では「venv_test」ディレクトリ)直下に環境を作ることが一般的なようです。

アクティベート

コマンドプロンプトの場合、下記のコマンドで仮想環境をアクティベートできます
(カレントフォルダの直下に環境作成した場合のコマンドなので、適宜パスは読み替えてください)

> .\[環境名]\Scripts\activate

下記のようにアクティベートした環境名がシェルの最初にカッコ表示されれば成功です

(newvenv) [カレントフォルダ名]>

deactivateコマンドで、仮想環境を終了(デアクティベート)できます

(newvenv) > deactivate



※ちなみにPowerShellでアクティベートコマンドを実行すると、「このシステムではスクリプトの実行が無効になっているため、‥を読み込むことができません」というようなエラーが出ることがあります。
こちらのリンクのように、 PowerShell のスクリプトの実行が実行ポリシーによって許可されていないことが原因なので、最初の1回のみ

> Set-ExecutionPolicy RemoteSigned -Scope CurrentUser -Force

としてあげれば、スクリプトの実行が有効となり

> .\[環境名]\Scripts\activate.ps1

としてあげれば、アクティベートされます

ライブラリのインストール

下記のようにpip install [ライブラリ名]で、ライブラリをインストールできます

(newvenv) >pip install pandas

pip listで、インストールされているライブラリのバージョンが確認できます

(newvenv) >pip list
Package         Version
--------------- -------
numpy           1.19.4 
pandas          1.1.4  
pip             20.2.1 
python-dateutil 2.8.1  
pytz            2020.4 
setuptools      49.2.1 
six             1.15.0 

VSCodeの設定

実行のたびにアクティベートコマンドを打つのは面倒なので、VSCodeでデバッグ実行時は勝手に仮想環境を読み込んでくれるようにしたいです。

こちらのリンクを参考に仮想環境の読込先を指定してください

→なんとプロジェクトフォルダ(VSCodeで開くフォルダ)直下に仮想環境を作っておけば、仮想環境を勝手に読み込んでくれるようです!VSCodeさん気が利きすぎですね!
image.png
なお、上の赤枠に現在デバッグに使用されているPythonのバージョンが表示されており、クリックするとデバッグ環境を切り替えることもできます


以上で仮想環境の使用ができるようになりました。以下は必要に応じて実行してください

仮想環境のコピー

下記コマンドで、パッケージ構成をrequirements.txtに出力できます

(newvenv) >python -m pip freeze > requirements.txt

このテキストファイル内に、インストールされているパッケージ一覧が記載されています

requirement.txt
numpy==1.19.4
pandas==1.1.4
python-dateutil==2.8.1
pytz==2020.4
six==1.15.0

このテキストファイルを他のプロジェクトフォルダに移し、仮想環境をアクティベートしたのち下記コマンドを打つと、該当するパッケージ全てをまとめてインストールしてくれます

(newvenv2) >python -m pip install -r requirements.txt

チーム内でパッケージを共通化したい場合に便利ですね!

ちなみに、仮想環境に入らずに同様のコマンドを実行すると、デフォルト環境にrequirements.txtに記載されたパッケージをまとめてインストールしてくれます。

バージョンを指定してライブラリをインストール

pipenv install numpy==1.19.3

仮想環境の初期化

デアクティベートした状態で、下記コマンドで初期化できます

python -m venv --clear [環境名]

仮想環境の削除

デアクティベートした状態で、仮想環境のフォルダを丸ごと削除するだけでOKです。
エクスプローラで仮想環境のフォルダ(上の例では「newvenv」フォルダ)を削除するか、
PowerShellで環境の存在するディレクトリに移動して下記コマンドを実行します

Remove-Item [環境名] -Recurse -Force

1-2. 仮想環境でインタープリタを変えたい場合

インタープリタ(Pythonのバージョン)をデフォルトから変えたい場合、
複数のバージョンのPythonをWindowsにインストールすることで対応できます。

Windowsでpyenvが不要と言われる所以はここにあります!
(逆に複数の3.0系PythonがインストールできないLinuxでは、pyenvが必要となることもあります)

手順を下記します

デフォルト以外バージョンのPythonインストール

Python公式から好きなバージョンのWindows x86-64 executable installerをダウンロードし、実行してインストールします。
デフォルト以外のPythonインストール時は、下図赤枠Add Python3.* to Pathをチェックしないでください
image.png

下図のように、複数のバージョンのPythonがインストールされていれば成功です
(通常はC:\Users\[ユーザ名]\AppData\Local\Programs\Pythonにインストールされるようです)
image.png

インタプリタがデフォルトと異なる仮想環境の構築

以降は、デフォルトのPythonバージョンが3.8.6、仮想環境のPythonバージョンを3.9.0である前提で解説します。

下記コマンドで、Pythonバージョンを変えた仮想環境を作成できます
(デフォルトインタプリタではpython -mとしていた部分を、py -[仮想環境に適用したいPythonのバージョン] -mに変えます)

> cd [環境構築したいディレクトリ]
> py -3.9 -m venv [新しい環境名]

デフォルト状態(仮想環境アクティベート前)のPythonバージョン確認

> python -V
Python 3.8.6

作成した仮想環境のPythonバージョン確認(環境名は「py39env」とした)

> .\py39env\Scripts\activate.ps1
(py39env) > python -V
Python 3.9.0

仮想環境においてPythonインタプリタのバージョンを変えられていることが分かります

1.3 仮想環境不使用時のライブラリのインストール

仮想環境不使用(デフォルトのPython環境使用)でライブラリをインストールしたい場合は、
アクティベートされていない状態(デフォルト状態)でpip install [ライブラリ名]と打つだけでOKです

> pip install pandas

pip listで、インストールされているライブラリのバージョンが確認できます

>pip list
Package         Version
--------------- -------
numpy           1.19.4 
pandas          1.1.4  
pip             20.2.1 
python-dateutil 2.8.1  
pytz            2020.4 
setuptools      49.2.1 
six             1.15.0 

2. pipenvの使用法

VSCodeのシェルを使って作業する前提で、使用手順を解説します。

【参考にさせて頂いた記事】
https://qiita.com/subarunari/items/dcbdad68ba1809b721b9
https://qiita.com/y-tsutsu/items/54c10e0b2c6b565c887a
https://pipenv-ja.readthedocs.io/ja/translate-ja/

Anacondaのアンインストール

こちらの「方法2」で、完全にアンインストールしてください

PythonとVSCodeのインストール

こちらを参考にPythonとVisual Studio Codeをインストールしてください
PyCharm等を使用したい場合は、VSCodeのインストールは飛ばして頂いて構いません

Pythonインストール時は、下図赤枠Add Python3.* to Pathをチェックすると、デフォルトでのPython環境に設定できます
image.png

また、VSCodeインストール時は下記2箇所をチェックするとエクスプローラの右クリックで開けて便利です
image.png

pipenvのインストール

下記コマンドで、pipenvをインストールします

> pip install pipenv

仮想環境の構築と使用法

下記のように
仮想環境の作成 → ライブラリのインストール → VSCodeの設定
の順で進めてください

pipenv仮想環境の作成

下図のような感じで環境構築したいディレクトリ(プロジェクトディレクトリ)を準備しておきます
image.png

上記の「プロジェクトディレクトリ」は、VSCodeのようなエディタで開くディレクトリ(上例では「pipenv_test」)であり、1. venv+pipのときと同様、この直下に仮想環境を作ることが一般的なようです。

ディレクトリ内に新環境を作成します

> cd [環境構築したいディレクトリ]
> pipenv install

※注意したいのが、OS内に複数のPythonをインストールしている場合です
ここで使用されるPythonバージョンはパスを通している(インストール時にAdd Python * to Pathをチェックした)デフォルトのバージョンではなく、インストールされている中で最新のバージョンとなります。
例えばPython3.8.6と3.9.0をインストールしており、3.8.6の方にパスを通していたとしても、pipenvではPython3.9.0が仮想環境に採用されてしまいます。

バージョンを指定して新環境を構築したい場合pipenv --python [Pythonバージョン]のように打ちます

> pipenv --python 3.8

指定したバージョンがOS内に存在する場合はそのPythonが使用され、存在しない場合には新たにそのバージョンのPythonがインストールされます

環境作成が成功すれば、PipfileおよびPipfile.lockというファイルが生成するはずです
(Pipfile.lockは生成しない事もありますが、その場合も後でライブラリインストール時に生成するので気にしないでください)
image.png

Pipfileの中身は下記のようになっています

[[source]]
url = "https://pypi.org/simple"
verify_ssl = true
name = "pypi"

[packages]

[dev-packages]

[requires]
python_version = "3.8"

ライブラリのインストール

pipenv install [ライブラリ名]コマンドで、ライブラリをインストールできます。

>pipenv install pandas

インストール後のPipfileを見ると、packagesにpandasが加わっている事が分かります。
このようにPipfileには、インストールしたパッケージの履歴が記録されます

[[source]]
url = "https://pypi.org/simple"
verify_ssl = true
name = "pypi"

[packages]
pandas = "*"

[dev-packages]

[requires]
python_version = "3.8"

Pipfile.lockを見ると、下記のように何やら色々と書かれていますが、ここにライブラリ間の依存関係等を記録しているようです。

Pipfile.lock(一部抜粋)
        },
        "pandas": {
            "hashes": [
                "sha256:09e0503758ad61afe81c9069505f8cb8c1e36ea8cc1e6826a95823ef5b327daf",
                "sha256:0a11a6290ef3667575cbd4785a1b62d658c25a2fd70a5adedba32e156a8f1773",
                "sha256:0d9a38a59242a2f6298fff45d09768b78b6eb0c52af5919ea9e45965d7ba56d9",
                "sha256:112c5ba0f9ea0f60b2cc38c25f87ca1d5ca10f71efbee8e0f1bee9cf584ed5d5",
                "sha256:185cf8c8f38b169dbf7001e1a88c511f653fbb9dfa3e048f5e19c38049e991dc",
                "sha256:3aa8e10768c730cc1b610aca688f588831fa70b65a26cb549fbb9f35049a05e0",
                "sha256:41746d520f2b50409dffdba29a15c42caa7babae15616bcf80800d8cfcae3d3e",
                "sha256:43cea38cbcadb900829858884f49745eb1f42f92609d368cabcc674b03e90efc",
                "sha256:5378f58172bd63d8c16dd5d008d7dcdd55bf803fcdbe7da2dcb65dbbf322f05b",
                "sha256:54404abb1cd3f89d01f1fb5350607815326790efb4789be60508f458cdd5ccbf",
                "sha256:5dac3aeaac5feb1016e94bde851eb2012d1733a222b8afa788202b836c97dad5",
                "sha256:5fdb2a61e477ce58d3f1fdf2470ee142d9f0dde4969032edaf0b8f1a9dafeaa2",
                "sha256:6613c7815ee0b20222178ad32ec144061cb07e6a746970c9160af1ebe3ad43b4",
                "sha256:6d2b5b58e7df46b2c010ec78d7fb9ab20abf1d306d0614d3432e7478993fbdb0",
                "sha256:8a5d7e57b9df2c0a9a202840b2881bb1f7a648eba12dd2d919ac07a33a36a97f",
                "sha256:8b4c2055ebd6e497e5ecc06efa5b8aa76f59d15233356eb10dad22a03b757805",
                "sha256:a15653480e5b92ee376f8458197a58cca89a6e95d12cccb4c2d933df5cecc63f",
                "sha256:a7d2547b601ecc9a53fd41561de49a43d2231728ad65c7713d6b616cd02ddbed",
                "sha256:a979d0404b135c63954dea79e6246c45dd45371a88631cdbb4877d844e6de3b6",
                "sha256:b1f8111635700de7ac350b639e7e452b06fc541a328cf6193cf8fc638804bab8",
                "sha256:c5a3597880a7a29a31ebd39b73b2c824316ae63a05c3c8a5ce2aea3fc68afe35",
                "sha256:c681e8fcc47a767bf868341d8f0d76923733cbdcabd6ec3a3560695c69f14a1e",
                "sha256:cf135a08f306ebbcfea6da8bf775217613917be23e5074c69215b91e180caab4",
                "sha256:e2b8557fe6d0a18db4d61c028c6af61bfed44ef90e419ed6fadbdc079eba141e"
            ],
            "index": "pypi",
            "version": "==1.1.4"
        },
        "python-dateutil": {
            "hashes": [
                "sha256:73ebfe9dbf22e832286dafa60473e4cd239f8592f699aa5adaf10050e6e1823c",
                "sha256:75bb3f31ea686f1197762692a9ee6a7550b59fc6ca3a1f4b5d7e32fb98e2da2a"
            ],
            "markers": "python_version >= '2.7' and python_version not in '3.0, 3.1, 3.2'",
            "version": "==2.8.1"
        },

venv+pipでは手動でrequirements.txtを作る必要があったのが、pipenvではPipfileとPipfile.lockに自動でライブラリインストールの履歴を残せるので、バックアップやコピー時に便利です。
Gitで管理すればチームでの共有も容易になりますね!

VSCodeの設定

VSCodeでのデバッグ実行時に勝手に仮想環境を読み込んでくれるようにします。
最初の1回のみ下図の手順で仮想環境を指定してあげれば、次回以降は勝手に設定が読み込まれます
image.png


以上で仮想環境の使用ができるようになりました。以下は必要に応じて実行してください

pipenv各種操作

pipenvでよく使う操作を記載します
多くの操作は、下記コマンドで仮想環境内に入って実行します

> pipenv shell

仮想環境でのスクリプト実行

仮想環境に入ったのち、python [スクリプト名].pyで実行できます

> pipenv shell
> python script1.py

環境のコピー

ここがvenv+pipと比べた際のpipenv最大のメリットですが、
環境のコピーを容易に実施することができます。

・ただ同じパッケージをインストールしたい場合は、Pipfileを同フォルダに置いて
pipenv install

>pipenv install

・パッケージのバージョンや履歴も含め詳細に再現したい場合は、Pipfile.lockを同フォルダに置いてpipenv sync

>pipenv sync

また、少し進んだ使い方として、開発環境のみのパッケージを管理することもできるようです

requirement.txtからのインストール

1.venv+pipで使用していたrequirement.txtからまとめてパッケージインストールすることもできます

> pipenv shell
> pipenv install -r ./requirements.txt

インストールされているライブラリの確認

仮想環境内に入ってpip listで仮想環境内にインストールされているライブラリ一覧を確認できます

> pipenv shell
> pip list

バージョンを指定してライブラリをインストール

pipenv install numpy==1.19.3

仮想環境の削除

仮想環境に入らずにpipenv --rmで削除できます

pipenv --rm

仮想環境の初期化

削除して再作成すれば初期化できます

pipenv --rm
pipenv install

3. Anaconda有償版

会社ごとに導入フローが異なると思うので、会社のIT部門に相談して下さい
使用法に関しては、基本的には無償版と同じはずです (仮想環境の構築法はこちらの記事が非常にわかりやすいです)

参考:有償版の価格設定

おわりに

・仮想環境は結構容量を食う(pandasやscikit-learnがデカい)ので、ご注意ください

・本記事とは直接関係ありませんが、2020/12現在Numpyをpipでインストールすると下記の問題が起こるので、バージョンを下げて対応してください。
https://qiita.com/bear_montblanc/items/b4b75dfd77da98076da5
(Anacondaはこの手の問題が解決した比較的Stableなバージョンをインストールしてくれるので、その意味での導入メリットもあるかもしれませんね)

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

すぐに使える!様々なデータへのクラスタリングとその評価方法〜行列・グラフ・時系列〜

はじめに

はじめまして、NTTドコモ R&D Advent Calendar 2020の10日目の記事を担当いたします、廣川です!
業務では,マーケティング分野を中心にデータ分析に取り組んでいます。私は、クラスタリングに縁があるようで、気がつくと様々なデータに対するクラスタリング手法を学んでいました。そこで本記事では,それらのクラスタリング手法についてご紹介したいと思います!

本記事の目的

  • 行列・グラフ・時系列データに対するクラスタリング手法をご紹介します。
  • クラスタ数の決め方や評価指標などについても簡単にご紹介します。

ご紹介する手法

ソースコードあり

  • 行列:Latent dirichlet allocation (LDA).
  • グラフ:Spectral clustering, Clauset-Newman-Moore greedy modularity maximization, Louvain method.
  • 時系列:Hidden Markov model (HMM).

概要のみ

  • 行列:Non-negative matrix factorization (NMF).
  • グラフ:A structural clustering algorithm for networks (SCAN).
  • 時系列:Dynamic time warping (DTW), Toeplitz Inverse Covariance-Based Clustering of Multivariate Time Series Data(TICC).

評価指標(ソースコードあり)

  • ラベルあり:Adjusted mutual information (AMI)
  • 自然言語処理:Perplexity, Coherence.
  • グラフ:Modularity, Conductance.
  • 確率モデル:Bayesian information criterion (BIC), log likelihood.

クラスタリングとは?

クラスタリングは、教師なし学習と呼ばれる手法の1つで、大規模なデータを少数のグループにまとめる手法です。主な用途としては、解釈性の向上やそれに伴うインサイトの獲得(データマイニング)。また、計算リソースの削減や計算モデルが扱いやすくすること(次元削減・特徴量抽出)等があげられます。

最近は、次元削減や特徴量抽出の用途としては、分散表現が用いられることも多いと思いますが、それでも、クイックな分析が求められるケースや、手元にある計算リソースの関係などで、クラスタリング手法が用いられるケースもあります。

教師あり学習と比較すると、教師ラベルがなくても分析が可能という大きなメリットがあります。一方で、正解データがないことから、何をもって「よいクラスタ」や「よい手法」とするのか?という評価が難しい点がデメリットとして考えられます。このような背景から、よく使われるクラスタリングの評価指標についても簡単にご紹介していきます。

行列データのクラスタリング

それでは、まず行列データのクラスタリングから見ていきましょう!
行列データのクラスタリングというと、トピックモデルと、行列因子分解(Matrix Factorization)が有名です。両者とも同様の問題設定となっており、大きな行列を2つの小さな行列に分解する手法と考えることができます1

MF.png

トピックモデルは自然言語処理の分野で用いられてきた手法で、単語と文章の共起関係に基づきクラスタリングを行います。Pythonですぐに使えるトピックモデルとしては、Latent Dirichlet Allocation (LDA)や、LDAのトピック生成部分を階層化することで、クラスタ数も同時に学習できるようなになったHDP−LDAがあります。

もう一方の行列因子分解は、協調フィルタリングというレコメンドシステムの用途として有名になった手法です。また、次元削減後に、元の行列を復元することで、欠損値補間の用途として利用されることもあります。こちらもPythonの利用と、マーケティング分野などでの応用を考えると、非負制約を追加したNon-negative Matrix Factorizationが扱いやすいと思います。

このように、異なるバックグラウンドから生まれた2つの手法ですが、実はこの2つの手法はかなり似ているようです。厳密には、LDAのベイズ化前であるpLSAと一般化KLダイバージェンス基準のNMFは等価です。

今回は、scikit-learnにある20 newsgroups datasetを用いて、LDAのデモンストレーションを行いたいと思います。

LDAによるクラスタリング

import pandas as pd
import gensim
from gensim.utils import simple_preprocess
from gensim.parsing.preprocessing import STOPWORDS
from gensim import corpora
from gensim.corpora.dictionary import Dictionary
from gensim.models.ldamodel import LdaModel

from nltk.stem import WordNetLemmatizer, SnowballStemmer
from nltk.stem.porter import *
import numpy as np
import nltk

from sklearn.datasets import fetch_20newsgroups
from sklearn.metrics import adjusted_mutual_info_score

np.random.seed(2020)
nltk.download('wordnet')
stemmer = SnowballStemmer('english')

データの読み込み

dataset = fetch_20newsgroups(shuffle=True, random_state=32, remove=('headers', 'footers', 'qutes'))

documents = pd.DataFrame({'News': dataset.data, 'Target': dataset.target})

print(documents.shape)
documents.head()
# (11314, 2)

前処理

前処理は、こちらが非常にわかりやすかったため、ほぼそのまま利用させてもらっています。
行った前処理は、以下の内容です。

  • ドキュメントのトークン化:具体的には,テキストを文単位に,文を単語単位にそれぞれ分割します.大文字の単語を小文字に修正し,句読点を取り除きます.
  • 3文字未満の単語の除去
  • ストップワードの除去:ストップワードとは,頻出あるいは一般的すぎて,利用価値が低いと考えられる単語のことをいいます.ex. the, a, for
  • 単語を見出し語化 (lemmatize) :具体的には,三人称の単語は一人称に,過去時制及び未来時制の動詞は現在時制に変更されます.
  • 単語の語幹抽出 (stemming)
def lemmatize_stemming(text):
    return stemmer.stem(WordNetLemmatizer().lemmatize(text, pos='v'))

def preprocess(text):
    result = []
    for token in gensim.utils.simple_preprocess(text):
        if token not in gensim.parsing.preprocessing.STOPWORDS and len(token) > 3:
            result.append(lemmatize_stemming(token))
    return result


# 前処理のテスト
# 適当な文で前処理の前後を比較
doc_sample = documents.iloc[0].News

words = []
for word in doc_sample.split(' '):
    words.append(word)

print('original document: ')
print(words)
print('\n\ntokenized and lemmatized document: ')
print(preprocess(doc_sample))

# original document: 
#['Please', 'could', 'someone', 'in', 'the', 'US', 'give', 'me', 'the', 'current', 'street', '\nprices', 'on', 'the', 'following,', 'with', 'and', 'without', 'any', 'relevant', 'taxes:\n\n', '8', 'Mb', '72', 'pin', 'SIMM\n16', 'Mb', '72', 'pin', 'SIMM', '(both', 'for', 'Mac', 'LC', 'III)\n\nAre', 'any', 'tax', 'refunds', 'possible', 'if', 'they', 'are', 'to', 'be', 'exported\nto', 'the', 'UK?', 'Can', 'you', 'recommend', 'a', 'reliable', 'supplier?']

#tokenized and lemmatized document: 
#['current', 'street', 'price', 'follow', 'relev', 'tax', 'simm', 'simm', 'refund', 'possibl', 'export', 'recommend', 'reliabl', 'supplier']

# すべてのデータに適用
processed_documents = documents['News'].map(preprocess)
processed_documents[:10]

前処理後の文章が、シンプルな単語(語幹)の集合になっていることがわかります。最後に、全データに前処理を適用しています。

次に、データを学習用とテスト用に分けます。

# データサイズの確認
len(processed_documents)

# 11314

# 学習用データとテスト用データに分ける
# すでにシャッフル済みのため、そのまま利用する
train_set = processed_documents[:10000]
test_set = processed_documents[10000:]
print(len(train_set), len(test_set))

# 10000 1314

# LDA用にフォーマットを整理する
dictionary = gensim.corpora.Dictionary(processed_documents)
dictionary.filter_extremes(no_below=15, no_above=0.5, keep_n=100000)

train_corpus = [dictionary.doc2bow(doc) for doc in train_set]
test_corpus = [dictionary.doc2bow(doc) for doc in test_set]

クラスタ数の決定

クラスタ数の決定は、応用上よく遭遇する問題だと思います。あくまで私の経験ですが、以下のように決めることが多いです。

  1. 今回の問題のように事前知識からクラスタ数の目安がある
  2. いくつかのクラスタ数で実験し定性評価で決定
  3. クラスタ数まで最適化するモデルを利用する(HDP-LDA等)
  4. 目的が別にある場合(予測タスクなど)、それに基づくCross Validation
  5. 確率モデルならBIC等の情報量規準を利用する

クイックな案件やデータの概要をつかみたいような探索的なモチベーションの場合は、1.や2.を使うこともあります。一方で、論文を書く場合など客観性が求められる場合は、3.〜5.のうち適用可能なものを選択しています。

今回は、正しいクラスタ数(20)が与えられていることから、こちらを利用したいと思います。参考として、クラスタ数が10、30のときの振る舞いも確認します。
それでは、クラスタ数を10、20、30でLDAを実行します。

lda_10 = gensim.models.LdaMulticore(train_corpus, num_topics=10, id2word=dictionary, passes=2, workers=2)
lda_20 = gensim.models.LdaMulticore(train_corpus, num_topics=20, id2word=dictionary, passes=2, workers=2)
lda_30 = gensim.models.LdaMulticore(train_corpus, num_topics=30, id2word=dictionary, passes=2, workers=2)

AMIとPerplexity/Coherenceによる定量評価

クラスタリング手法の評価をする場合、教師ラベルがないことから、データとの類似度を指標とすることが一般的です。また、クラスタの利用目的がはっきりしている場合には、その目的やタスクに合わせた独自の指標が定義されることもあります。具体例としては、トピックモデルをトピック品質の観点から評価するCoherence等が該当します。

今回は、正解データ(20のラベルデータ)が存在することとから、Adjusted Mutual Information(AMI)によって評価します。また参考として、トピックモデルの評価指標としてよく利用されるPerplexityとCoherenceも見てみます。
このAMIという指標は、簡単に言うと教師あり学習の評価で利用するAccuracyのようなものです2。 ただし、クラスタリングの場合には、教師あり学習とは異なり、同じクラスタ割当をしていたとしてもラベルが異なることがあり得ます。例えば、以下の2つの割当を考えます。

  • Aさん/Bさん/Cさんをクラスタ1、Dさん/Eさんをクラスタ2とする割当
  • Aさん/Bさん/Cさんをクラスタ2、Dさん/Eさんをクラスタ1とする割当

この2つは、クラスタリング結果としては等価ですが、ラベルが異なっています。そのため、通常のAccuracyで評価すると、性能が本来よりも低く見積もられてしまいます。そこで、ラベルを無視し、純粋に割当のみを評価する指標がAMIです。具体的には、異なる割当の時ほど、値が0に近く、この例のように全く等価の場合、1を示します。
それでは、LDAの結果を見てみましょう。

# AMIによる定量評価
dist = lda_20.get_document_topics(test_corpus)
clusters = [max(lis,key=lambda item:item[1])[0] for lis in dist]
print(f"{adjusted_mutual_info_score(clusters, documents[10000:].Target):.3f}")

# 0.377

AMIが高くないことから、正解ラベルの意図するようなクラスタ割当になっていない可能性が高いです。今回のように求める割当がはっきりしているケースで、教師データや特徴量も潤沢にあるのであれば、教師あり学習を用いる方がよいかもしれません。一方で、教師データが少ない場合には、モデルの改良(ハイパーパラメータの調整や制約の追加など)によって、目的のラベルに近づけることができるかもしれません。
今回の目的は、あくまでデータに対する探索的なアプローチなので、これ以降は教師データがない状況を想定して分析を続けます。

次に、クラスタ数を10、20、30と変化させた際のPerplexityとCoherenceの変化を見てみます。Coherenceには、複数の種類がありますが、今回は、対象データのコーパスのみから算出可能なU-mass Coherenceを用いています。

# perplexity
perplexity = np.exp2(-lda_20.log_perplexity(test_corpus))
print(f'Perplexity score: {perplexity:.2f}')

# coherence(u-mass)
cm = gensim.models.CoherenceModel(lda_20,corpus=test_corpus,dictionary=dictionary,coherence='u_mass')
print(f'u_mass coherence score: {cm.get_coherence():.2f}')

# perplexity score: 590.84
# u_mass coherence score: -2.61
#clusters Perplexity Coherence
10 352.53 -1.99
20 590.83 -2.62
30 1256.94 -2.25

あくまで1回の実験結果ですが、PerplexityやCoherenceの観点からは、クラスタ数:10が良さそうです。

トピックの上位単語による定性評価

最後に定性的な分析を行います。定性的な分析(解釈)には、それなりの時間と労力がかかることから、私は何かしらの定量的な指標でモデルの性能を評価してから、定性評価(解釈)することにしています。その方が、手戻りが少なくなることに加え、その結果が「新発見なのか?」「何かのミスなのか?」というところを悩むコストが低減されるためです。

それでは、各トピックの上位単語を可視化してみましょう!
※クラスタ数は10のものを可視化します。

for idx, topic in lda_10.print_topics(-1, num_words=5):
    print('Topic {}: {}'.format(idx, topic))

# Cluster 0: 0.008*"know" + 0.008*"peopl" + 0.006*"like" + 0.005*"think" + 0.005*"time"
# Cluster 1: 0.011*"drive" + 0.011*"articl" + 0.010*"like" + 0.008*"know" + 0.006*"think"
# Cluster 2: 0.011*"space" + 0.007*"articl" + 0.007*"launch" + 0.006*"nasa" + 0.006*"like"
# Cluster 3: 0.012*"game" + 0.008*"articl" + 0.008*"team" + 0.007*"window" + 0.007*"year"
# Cluster 4: 0.009*"peopl" + 0.007*"right" + 0.007*"christian" + 0.007*"believ" + 0.007*"articl"
# Cluster 5: 0.011*"file" + 0.010*"window" + 0.007*"program" + 0.006*"articl" + 0.006*"like"
# Cluster 6: 0.007*"like" + 0.007*"articl" + 0.007*"year" + 0.007*"good" + 0.006*"know"
# Cluster 7: 0.010*"think" + 0.009*"articl" + 0.008*"time" + 0.007*"know" + 0.007*"say"
# Cluster 8: 0.013*"articl" + 0.006*"time" + 0.005*"orbit" + 0.005*"think" + 0.005*"like"
# Cluster 9: 0.011*"peopl" + 0.010*"armenian" + 0.007*"say" + 0.006*"know" + 0.006*"think"

少しだけ中身を見てみます。今回のクラスタ2はspaceやlaunch、nasa等の宇宙系のクラスタのように見えます。また、クラスタ5は、file, window, programとありますので、コンピュータ系のクラスタのようです。このように、教師データがなくても、クラスタ数や各クラスタの特徴など、ある程度データの特徴が判明しました。

グラフデータのクラスタリング

それでは続いて、最近のトレンドであるグラフデータについて、クラスタリングの観点から見てみましょう!
まず、グラフデータのクラスタリングを、以下のように整理してみます。どちらも出力されるのはNodeのコミュニティ(集合)であり、Linkの集合ではありません。

  • Nodeベースのクラスタリング
  • Linkベースのクラスタリング

Nodeベースのクラスタリングでは、A Structural Clustering Algorithm for Networks(SCAN)と呼ばれるアルゴリズムがあります。こちらは2つのNode間でどれだけ同じNodeに隣接しているかという指標(Structural Similarity)に基づきクラスタリングを行います。

Linkベースのクラスタリングは、スペクトラルクラスタリング(Normalized cut等)や、コミュニティ検出(モジュラリティベースの手法など)等が存在します。特に、コミュニティ検出と呼ばれるアルゴリズムでは、自動でクラスタ数が定まるという特徴があります。 例えば、モジュラリティに基づくクラスタリング手法では、以下のモジュラリティという量を最大化する問題を考えます。

Q=\sum_{l \in C}(e_{ll}-a_l^2)\\
e_{lm} = \frac{1}{2M}\sum_{i \in V_l}\sum_{j \in V_m} A(i,j)\\
a_l = \sum_{m \in C} e_{lm} = \frac{1}{2M}\sum_{i \in V_l}\sum_{j \in V}A(i,j)

ここで、$C$は検出したいコミュニティの集合であり、$i,j$はNode、$M$はLinkの総数、$V_l, V_m$はコミュニティ$l, m$に含まれるNode集合、$A(i,j)$は対象グラフの隣接行列を表します。また、$e_{lm}$はコミュニティ$l, m$間に存在するLinkの割合、$a_l$はコミュニティ$l$内のNodeに接続するLink数の割合を表します。

この指標は定性的には、同じクラスタに所属するすべてのLinkに対して、ランダムグラフで期待されるものより、どれだけ多くのLinkが存在しているかを評価しています。先にも述べましたが、この手法では、クラスタ集合Cを含めて最適化されるため、ユーザがクラスタ数を指定する必要はありません。

Modularityによるコミュニティ検出

それでは、Modularityによるコミュニティ検出の実験を行います。具体的には、Louvain methodと呼ばれる手法と、Clauset-Newman-Moore greedy modularity maximizationという手法を用いてコミュニティ検出を行います。

Louvain methodの方が高速・高精度な手法として有名です。Greedy modularity maximizationの方は、今回のグラフのサイズでは問題なく動作しますが、もう少し大きくなると計算負荷が高くなる印象です。
今回は、Stochastic Block modelという手法で人工的にデータを生成します3

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from collections import defaultdict

import networkx as nx
from networkx.algorithms.community import greedy_modularity_communities
from sklearn.cluster import SpectralClustering
from cdlib import algorithms
from networkx.algorithms.cuts import conductance

可視化用の関数を定義

def show_graph(g, node2community=None, k=0.1):
    """\
    可視化用の関数

    Parameters
    ----------
    g : Graph
        対象のグラフ
    node2community : dict
        Node IDとコミュニティIDの辞書
    k : float
        可視化パラメータ
    """
    plt.figure(figsize=(10, 8))
    pos = nx.spring_layout(g, k=k)

    # コミュニティを指定したとき
    if node2community is not None:
        cmap = cm.get_cmap('viridis', max(node2community.values()) + 1)
        nx.draw_networkx_nodes(g, pos, node2community.keys(), node_size=40, 
                               cmap=cmap, node_color=list(node2community.values()))
    else:
        nx.draw_networkx_nodes(g, pos, node_size=40)

    nx.draw_networkx_edges(g, pos, alpha=0.5)

    plt.axis('off')
    plt.show()

Stochastic Block modelによるグラフ生成と可視化

# パラメータ設定
## 対角成分(クラスタ内)を密に、他を疎にしています。
sizes = [75, 100, 300]
probs = [
    [0.1, 0.001, 0.002], 
    [0.001, 0.2, 0.004],
    [0.002, 0.004, 0.5]
]
g = nx.stochastic_block_model(sizes, probs, seed=0)

# 可視化
show_graph(g, k=0.1)

toy_graph.png

Modularityによるコミュニティ検出を実行します(Greedy modularity maximization)。

community = list(greedy_modularity_communities(g))
print("コミュニティ数: ", len(community))
print("コミュニティラベル: ", list(range(len(community))))

# コミュニティ数:  3
# コミュニティラベル:  [0, 1, 2]

コミュニティ数を設定せずに、きちんと3つのクラスタが検出されたことがわかります!
続いて、Louvain methodによるコミュニティ検出を実行します。

com = algorithms.louvain(g, weight='weight', resolution=1., randomize=True)
node2community = com.to_node_community_map()
node2community = {k:v[0] for k, v in node2community.items()}

community = defaultdict(list)
for i in range(len(node2community)):
    community[node2community[i]].append(i)
community = community.values()

print("コミュニティ数: ", len(community))
print("コミュニティラベル: ", list(range(len(community))))

# コミュニティ数:  5
# コミュニティラベル:  [0, 1, 2, 3, 4]

今度は、5つのコミュニティが検出されてしまいました。しかし、この結果は次のModularityを確認してみると、意外とよい結果であることがわかります。

ConductanceとModularityによる定量評価

先ほどのモジュラリティだけでなく、グラフクラスタリングの評価指標としてよく使われるコンダクタンスも算出します。コンダクタンスは、Ncutに似た指標ですが、物理学をバックグラウンドにもつ指標です。具体的には、あるクラスタの頂点集合$C_i$のコンダクタンス$\phi(C_i)$は、以下のように定義されます。

\phi(C_i) = \frac{\sum_{i \in C_i}\sum_{j \notin C_i}A(i,j)}{min(\kappa_i, \overline{\kappa_i})}

ここで、$A(i,j)$は隣接行列、$\kappa_i, \overline{\kappa_i}$は、それぞれ、頂点集合$C_i$内、それ以外の頂点集合$\overline{C_i}$内の次数の和です。

そしてグラフ全体のコンダクタンス$\phi(G)$は、すべてのクラスタ(頂点集合)のコンダクタンスの最小値として定義されます。

\phi(G) = min_{C_i}(\phi(C_i))

Pythonコードとしては、下記のようにすると算出できます。

for c in list(community.values()):
    print(conductance(g, c))

print(f"modularity : ", modularity(g, community))

# conductance for 0: 0.06059465562664659 #<-最小値
# conductance for 1: 0.06469135802469136
# conductance for 2: 0.07278481012658228

# modularity : 0.09895468017762296

次に、比較対象として、Ncutのクラスタ数:2の場合と5の場合を作成します。コードとしては、クラスタ2のケースのみ記載します。

adjacency_matrix= nx.to_numpy_matrix(g)

num_community = 2
sc = SpectralClustering(num_community, affinity='precomputed', n_init=100, assign_labels='kmeans')
node2community = sc.fit_predict(adjacency_matrix)

# 可視化用
vis_node2community = {v : node2community[i] for i, v in enumerate(g.nodes)}

community = defaultdict(list)
for i in range(len(node2community)):
    community[node2community[i]].append(i)

for c in list(community.values()):
    print(conductance(g, c))

print(f"modularity : ", modularity(g, community))

# conductance for 0: 0.06059465562664659 <-最小値
# conductance for 1: 0.06059465562664659

# modularity : 0.09817090642685321
Greedy Louvain Ncut:2 Ncut:5
Conductance 0.0606 0.0647 0.0606 0.0647
Modularity 0.0990 0.1044 0.0982 0.0935

グラフ全体のコンダクタンスとして、最小値を採用しているため、真のクラスタ数:3とは異なるクラスタ数:2のNcutでも、モジュラリティベースの手法と同じコンダクタンスの値となっています。この現象は、こちらでも指摘されており、評価指標としてコンダクタンスを利用する場合は、ご自身の目的に適合するかを確認するとよいかもしれません。

他方のModularityを見てみると、Louvain法が最大となっていることがわかります。今回は、3つのクラスタになるようにタスクを設定したつもりだったため、トイデータの生成方法や指標が厳密な意味で合わなかったと考えられます。一方で、先ほどの図を見ると、定性的には3つのクラスタがあるようにも見えると思います。あくまで個人的な考えですが、このような場合は、インサイトを得るためと割り切り(クラスタ数:3を採用し)、施策設計などでカバーできればよいのではないかと思っています。

可視化による定性評価

今回検出されたコミュニティ(GreedyとLouvainの結果)をNodeの色を変えて可視化してみます。

## フォーマットを整えて、再度クラスタを色分けして可視化
node2community = {}
for i in range(len(community)):
    for n in community[i]:
        node2community[n] = i

# 可視化
show_graph(g, node2community, k=0.1)

toy_graph_cluster.png

Louvain.png

上図のGreedyの方の結果は、想定通り、3つの集団が検出されていることが定性的にも確認できました!下図のLouvainの方は、1つと想定していた真ん中のクラスタが複数に分割されていることがわかります。

時系列データのクラスタリング

最後に時系列データのクラスタリングについてご紹介いたします!
時系列データのクラスタリングは、大きく2つに分けられます。

  1. 時系列データ全体を1点として扱い、別の時系列との類似度を評価する whole time series clustering.
  2. 時系列データそのものを複数のセグメントに分割し、そのセグメントをクラスタリングするsubsequence time series (STS) clustering.

Whole time series (WTS) clusteringの代表的なものは、Dynamic Time Warping (DTW)と呼ばれる手法です。DTWは、まず、2つの時系列の各点間の距離を総当たりで計算します。次に、時系列のインデックスを表す2次元平面上で、最初のインデックス(0,0)から、それぞれの系列の最後のインデックス(N,M)までの経路を探索します。その際、重複を許容した上で距離の和が最小となるようにマッチングをとっていきます。

そのため、系列の長さが異なる場合でも適用可能であり、また、位相がずれていても、形状が似ている場合の類似度が高くなるという特徴があります。後者に関しては、時刻が重要な場合など、問題によって向き不向きがあると思いますので、適用の際は注意が必要です。

もう一つのSTS clusteringは、ユーザのセンサーデータを例に説明すると、WTS clusteringが、複数人の1日の行動や1週間の行動などに基づきユーザをクラスタリングするのに対し、STS clusteringは、ユーザ個人4の1日や1週間に基づき、時間帯をクラスタリングします。STS Clusteringのイメージ図を載せておきます(実際の結果ではありませんので値とクラスタラベルは無関係です)。

sts.png

出力される結果に関しても、WTSでは、よく走るユーザ集合や、階段の上り降りが多いユーザ集合などであるのに対して、STSでは、朝と夜は走っている、お昼は階段が多い等の時間帯集合です。このようにSTS clusteringでは、時間帯という概念が残るため、解釈性を保った次元削減手法としても利用できます(朝と夜に走りやすく、お昼は階段が多いユーザ等)。

STS clusteringの有名な手法として、KDD2017のRunner upとなったToeplitz Inverse Covariance-Based Clustering of Multivariate Time Series Data(TICC)があげられます。
この手法は、まず、各データのクラスタラベルを初期化します。次に、ラベルごとにToeplitz graphical lassoという手法でスパースなMRFの構造を推定します。そして、時系列的な近傍で変化が起こりにくい(temporal consistency)という仮定のもと、尤度とこのペナルティ項に基づき、近傍の時系列のノイズを除去しながら、各時点のデータにクラスタ(MRF)を割り当てていきます。この2つの処理を収束するまで繰り返すことで、最終的にクラスタシークエンスと、各クラスタの構造(スパースなMRF)が推定されます。

ここで、MRFのPrecision matrix (共分散行列の逆行列)を、各センサーをNodeとするグラフとして考えてみます。このとき、もし2つのNode間にLinkが張られていないならば、この2つのNodeは条件付き独立となる性質があります。これは、他のすべてのセンサーの値で条件付けした際に、この2つセンサー値は互いに影響しないことを意味します。そのため、スパースなMRFが得られたということは、そのクラスタ内のセンサーの依存関係が判明したということになり、クラスタ自身を解釈する手助けになります。

今回は、このTICCを使うことも考えましたが、実験のしやすさの観点から、Hidden Markov model(HMM)によるSTS clusteringを行うことにします。具体的には、HMMから生成した人工データに対する復元のデモンストレーションを行います

HMMによるSTS clustering

時系列的なデータとしての性質を強調するため、HMMの遷移行列を工夫しています。具体的には、HMMが少し有利になるように、一部の遷移が頻出しやすいように遷移確率に偏りを持たせました。

import numpy as np
from hmmlearn import hmm
from sklearn.mixture import GaussianMixture as GMM
from sklearn.metrics import adjusted_mutual_info_score
import statistics
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(2020)

人工データの生成

num_clusters = 5

model = hmm.GaussianHMM(n_components=num_clusters, covariance_type="full")
model.startprob_ = np.array([0.2, 0.3, 0.1,0.2,0.2])
model.transmat_ = np.array([
    [0.6, 0.01, 0.01,0.35,0.03],  # 0->0, 0->3
    [0.01, 0.5, 0.45, 0.01,0.03], # 1->1, 1->2 
    [0.45, 0.03, 0.5, 0.01,0.01], # 2->2, 2->0
    [0.03, 0.02, 0.2, 0.5,0.25],  # 3->3, 3->2, 3->4
    [0.01, 0.25, 0.03,0.01,0.7]  # 4->4, 4->1
])

model.means_ = np.array([[0.0, 0.0, 0.0],  [1.0, 2.0,3.0], [5.0,-10.0,0.0],[-5.0, 10.0,5.0], [10,5.0,0.0]])
model.covars_ = np.tile(np.identity(3), (5, 1, 1))

X, Z = model.sample(200)
train_X = X[:100]
test_X = X[100:]

train_Z = Z[:100]
test_Z = Z[100:]

print(train_Z)
# array([4, 1, 2, 0, 0, 0, 0, 0, 3, 4, 4, 4, 4, 4, 4, 1, 1, 2, 4, 4, 4, 4, 2, 2, 2, 2, 2, 0, 3, 2, 2, 0, 0, 3, 3, 3, 3, 4, 1, 1, 1, 2, 2, 1, 1, 2, 2, 0, 3, 2, 2, 2, 0, 0, 0, 0, 3, 3, 3, 3, 0, 0, 0, 0, 1, 2, 2, 0, 4, 4, 1, 1, 1, 1, 2, 0, 0, 3, 3, 3, 4, 4, 4, 1, 1, 2, 2, 0, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4])

クラスタ数の決定

今回もトイデータのため、クラスタ数はわかっていますが、実際のケースを意識して、BICを見てみたいと思います。

bic_clusters = []
# cluster数を2〜10でBICを確認する
for n_components in range(2,11):
    bics = []

    # 初期値を変えて10回試行する
    for i in range(10):
        hmm_model = hmm.GaussianHMM(n_components=n_components, covariance_type="full", n_iter=100)
        hmm_model.fit(train_X)
        score = hmm_model.score(train_X)

        n_features = hmm_model.n_features
        mean = hmm_model.means_
        covs = hmm_model.covars_

        # 共分散を対角行列にしているので、非ゼロ成分のみを見ています。
        num_parameters = len(np.where(mean*mean >= 0.1)[0]) +  len(np.where(covs*covs >= 0.01)[0]) + n_components*(n_components-1) + (n_components-1)

        cur_bic =np.log(len(train_X)) * num_parameters - 2 * score
        bics.append(cur_bic)

    bic_clusters.append(statistics.mean(bics))

plt.plot(list(range(2,11)), bic_clusters)

BIC.png

図のように、こちらが設定したクラスタ数:5でBICが最小となることが確認できました。

AMIと対数尤度による定量評価

こちらも、行列の場合と正解ラベルが存在するため、AMIを算出できます。また、通常の場合を想定して、データとの類似度としてテストセットに対する対数尤度も計算します。ベースライン手法として、時系列を考慮していないGaussian mixture model (GMM)を設定しました。
まずは、AMIを確認します。

ami_hmm = []
ami_gmm = []
for i in range(10):
    hmm_model = hmm.GaussianHMM(n_components=5, covariance_type="full", n_iter=100)
    hmm_model.fit(train_X)
    Z2 = hmm_model.predict(test_X)

    gmm_model = GMM(n_components=5).fit(train_X)
    Z3 = gmm_model.predict(test_X)

    ami_hmm.append(adjusted_mutual_info_score(test_Z,Z2))
    ami_gmm.append(adjusted_mutual_info_score(test_Z,Z3))

print(f"AMI achieved by HMM: {statistics.mean(ami_hmm):.3f}")
print(f"AMI achieved by GMM: {statistics.mean(ami_gmm):.3f}")

# AMI achieved by HMM: 0.972
# AMI achieved by GMM: 0.933

次に、テストセットに対する対数尤度を見てみます。

ll_hmm = []
ll_gmm = []
for i in range(10):
    hmm_model = hmm.GaussianHMM(n_components=5, covariance_type="full", n_iter=100)
    hmm_model.fit(train_X)
    score_HMM = hmm_model.score(test_X)

    gmm_model = GMM(n_components=5).fit(train_X)
    score_GMM = gmm_model.score(test_X)*len(test_X)

    ll_hmm.append(score_HMM)
    ll_gmm.append(score_GMM)

print(f"Log likelihood achieved by HMM: {statistics.mean(ll_hmm):.3f}")
print(f"Log likelihood achieved by GMM: {statistics.mean(ll_gmm):.3f}")

# Log likelihood achieved by HMM: -610.684
# Log likelihood achieved by GMM: -625.196
AMI Log likelihood
HMM 0.972 -610.684
GMM 0.933 -625.196

データ生成時の遷移行列を工夫したこともあり、HMMの方が高いAMIと対数尤度を達成していることがわかります。一方で、今回のトイデータを作成する際に平均値を大きく異なる値にしたことから、今回の人工データのクラスタ間の距離は大きいものと考えられます。そのため、時系列的な振る舞いを考慮していないGMMでも高精度に復元できているものと思われます。

HMMの遷移行列の可視化による定性評価

定性評価として、HMMの遷移行列を確認してみます。クラスタラベルが元の遷移行列と異なるため、容易には比較できませんが、値の傾向を見る限り、良さそうなフィッティングをしていそうなことが確認できました。

このようにデータのみに基づき、クラスタ数やクラスタの遷移傾向までつかむことができました。実際のセンサーデータの場合は、クラスタを表現するガウシアンのパラメータや発生している時刻なども参考に、クラスタ自体を解釈していきます。

sns.heatmap(hmm_model.transmat_, linewidths=.5, cmap="Reds", annot=True, fmt="1.1f")

heatmap.png

さいごに

最後までお読みいただきありがとうございます!
今回は、行列、グラフ、時系列データ等、様々なデータに対するクラスタリングについてご紹介しました。
それぞれの手法について定量的な側面と定性的な側面の2方向からアプローチし、データのみから分析可能な特徴についてご紹介しました。
3種類のデータ構造のうち、いずれか1つでも皆さまのビジネス・研究課題のお役に立てれば幸いです!

参考文献


  1. コア行列を追加し3つとすることもあると思いますが、一般的にはこちらのように思います。 

  2. 厳密には、偶然性を考慮して補正しています。 

  3. Stochastic block modelによるクラスタリングも可能ですが、今回は生成のみの目的で利用します。 

  4. 複数人を考慮することも可能です。ただし、手法によっては、ユーザ間でラベルを共有したい場合に、少し工夫する必要があります。 

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

FastAPIのStaticファイルの設定(WebSocket Chatの例を直す)

アドベントカレンダー8日目です。カレンダー空いてました!やばい!
昨日は@sh05_sh05vimでPDCAを回すでした!

FastAPIの静的ファイルの置き方法についてです。
扱うコードは公式ドキュメントのWebSocketの例を利用します。

アドベントカレンダーにしては狭いトピックですがFastAPIの紹介とWebSocketによるchatとドキュメントが充実しているよという意味合いに取ってください。

記事として冗長なので解決したくてアクセスされた場合は本編に飛んでください

FastAPIとは

Pythonの軽量ウェブフレームワークです。素晴らしさなどは他のページ、サイトで十分説かれているのでこの記事では公式より

  • Fast: Very high performance, on par with NodeJS and Go (thanks to Starlette and Pydantic). One of the fastest Python frameworks available.
  • Fast to code: Increase the speed to develop features by about 200% to 300%. *
  • Fewer bugs: Reduce about 40% of human (developer) induced errors. *
  • Intuitive: Great editor support. Completion everywhere. Less time debugging.
  • Easy: Designed to be easy to use and learn. Less time reading docs.
  • Short: Minimize code duplication. Multiple features from each parameter declaration. Fewer bugs.
  • Robust: Get production-ready code. With automatic interactive documentation.
  • Standards-based: Based on (and fully compatible with) the open standards for APIs: OpenAPI (previously known as Swagger) and JSON Schema.

FastAPI

と引用しておきます。
ドキュメントが自動で生成されたり、公式ドキュメントが充実していてとても好きです。

扱うコード

このページの最下部からです。
WebSockets - FastAPI

websocket_chat.py
from typing import List

from fastapi import FastAPI, WebSocket, WebSocketDisconnect
from fastapi.responses import HTMLResponse

app = FastAPI()

### 割愛 ###

@app.websocket("/ws/{client_id}")
async def websocket_endpoint(websocket: WebSocket, client_id: int):
    await manager.connect(websocket)
    try:
        while True:
            data = await websocket.receive_text()
            await manager.send_personal_message(f"You wrote: {data}", websocket)
            await manager.broadcast(f"Client #{client_id} says: {data}")
    except WebSocketDisconnect:
        manager.disconnect(websocket)
        await manager.broadcast(f"Client #{client_id} left the chat")

そのまま動かせばWebSocketを利用したChatアプリができます。すごいですね。

その中にこんな2箇所があります。

websocket_chat.py
html = """
<!DOCTYPE html>
<html>
    <head>
        <title>Chat</title>
    </head>
    <body>
        <h1>WebSocket Chat</h1>
        <h2>Your ID: <span id="ws-id"></span></h2>
        <form action="" onsubmit="sendMessage(event)">
            <input type="text" id="messageText" autocomplete="off"/>
            <button>Send</button>
        </form>
        <ul id='messages'>
        </ul>
        <script>
            var client_id = Date.now()
            document.querySelector("#ws-id").textContent = client_id;
            var ws = new WebSocket(`ws://localhost:8000/ws/${client_id}`);
            ws.onmessage = function(event) {
                var messages = document.getElementById('messages')
                var message = document.createElement('li')
                var content = document.createTextNode(event.data)
                message.appendChild(content)
                messages.appendChild(message)
            };
            function sendMessage(event) {
                var input = document.getElementById("messageText")
                ws.send(input.value)
                input.value = ''
                event.preventDefault()
            }
        </script>
    </body>
</html>
"""


@app.get("/")
async def get():
    return HTMLResponse(html)

明示的にHTMLと宣言して文字列を返しているようです。
もちろん注意書きで本番環境には適さない、フロントとちゃんと連携しろなどと書いてあります。

少しいじって試そうと思ったときに、(設定してなければ)シンタックスハイライトが効きづらいし、(折り畳みが無効であれば)ファイルが長くなります。
(特に強い理由はないですが練習としてやっていきます)

なので今回はFastAPIの静的ファイル、特にHTMLの配置をやります。

読むページ

Static Files - FastAPI
このページを見てできそうですが、ベースになっているStarLetteのドキュメントを読むとHTML用のオプションもあったのでそちら(下記)を読みます。
余談ですがresponderもStarletteがベースになっているそうです。

Static Files - Starlette
こちらを読んでみると

html - Run in HTML mode. Automatically loads index.html for directories if such file exist.

とあります。単なる静的ファイルだけでなくindex.htmlは読み込んでくれるそうです。

本編

構成

app
├── html
│   └── index.html
└── websocket_chat.py

aiofilesのインストール

pip install aiofiles

コード

websocket_chat.py
from typing import List

from fastapi import FastAPI, WebSocket
from fastapi.staticfiles import StaticFiles

app = FastAPI()


app.mount("/chat", StaticFiles(directory="/app/html", html=True), name="html")

# 以降は変更なし

app.mountがルーティングとなっています。
そして先述の通りindex.htmlを読み込むように指定したので@app.getの追加などをすることなく/chatにアクセスするとできます。
スクリーンショット 2020-12-08 17.17.15.png

おわり

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

【第5回】Python で競馬予想してみる  ~ LightGBM ~

前回までに、scikit-learn のロジスティック回帰を使って競馬予測のモデルを作ってみた。
説明変数などを工夫してみても馬券回収率(単勝)は 80%超程度までで上昇は見込めず

この回からは、流行りの決定木アルゴリズムに基づいた勾配ブースティング(Gradient Boosting)の機械学習フレームワーク LightGBM を使って競馬予想してみるのだ

データの取得

ターゲットフロンティアのレース検索データを出力して使用

前回までと同様にターゲットフロンティアのレース検索で 2000年からのレースをすべて検索
前走読込み ボタンを押して前走のレースID等を取得して、今回はすべてのデータを CSV で出力した
全データを出力すると出力時間や pandas での読込みに時間がかかるのだが、使いたいデータがないときの出力し直しが面倒なのでとにかく全データを出力したのだ

参照:【第1回】Python で競馬予想してみる ~ 環境設定 とデータの取得 ~

データの前処理

pandasでcsvファイル読み込み

使用する 特徴量 を選択

予想するレース(直近レース)で使用する特徴量は下記の通り

Index(['レースID(新)', '前走レースID(新)', '性別', '年齢', '父タイプ名', '母父タイプ名', '場所', 'クラス名',
       '距離', '芝・ダ', 'コーナー', '馬場状態', '出走頭数', '枠番', '馬番', '騎手コード', '斤量', 'キャリア',
       '間隔'],
      dtype='object')

上記の直近レースの 前走レースID(新) に紐づけで結合するための過去走で使用する特徴量は

Index(['レースID(新)', '前走レースID(新)', '場所', 'クラス名', '距離', '芝・ダ', 'コーナー', '馬場状態',
       '出走頭数', '枠番', '馬番', '騎手コード', '斤量', 'キャリア', '間隔', '確定着順', '斤量体重比',
       '単勝オッズ', '人気', '着差', '上3F地点差', '脚質', 'Ave-3F', '上り3F', '上り3F順', '馬体重増減',
       '補正', '補9'],
      dtype='object')

ダミー変数を作成

build.ipynb
# ダミー変数を生成
rd_last = pd.get_dummies(rd_last, columns=['場所', 'クラス名', '芝・ダ', '馬場状態', '性別', '父タイプ名', '母父タイプ名']) # 直近レース用
rd_past = pd.get_dummies(rd_past, columns=['場所', 'クラス名', '芝・ダ', '馬場状態', '脚質']) # 過去レース用

直近レースデータに加工走のデータをマージ

今回は 5走前まで連結

欠損値処理はしない

LightGBM は欠損値があっても学習・予測できるらしい
欠損値処理はいつでもできるので、とりあえず残しておくのだ

目的変数の処理

目的変数は、出走馬の着順を 3着以内なら 1、4着以下なら 0 に変更した

build.ipynb
# 着順が3着以内かどうかのカラムを追加する
f_ranking = lambda x: 1 if x in [1, 2, 3] else 0
df['3着以内'] = df['確定着順'].map(f_ranking)

参照:【第2回】Python で競馬予想してみる ~ データの前処理(CSVデータの読み込み~各種処理) ~

モデルの学習

アンダーサンプリングは必要 ?

アンダーさんプリンをした方が良いのか分からないので、今回はしない

build.ipynb
import lightgbm as lgb #LightGBM
# 説明変数をdataXに格納
dataX = df

# 目的変数をdataYに格納
dataY = df['3着以内']

# データの分割を行う(学習用データ 0.8 評価用データ 0.2) stratify:層化抽出
X_train, X_test, y_train, y_test = train_test_split(dataX, dataY, test_size=0.2, stratify=dataY)

# LightGBM用のデータセットに入れる
lgb_train = lgb.Dataset(X_train.values, y_train)
lgb_eval = lgb.Dataset(X_test.values, y_test)

# LightGBM parameters
params = {
    'task': 'train', 
    'objective': 'binary', 
    'metric': 'binary_logloss'
}

# モデルの学習
model = lgb.train(params, lgb_train, valid_sets=lgb_eval,
                  verbose_eval=50,  # 50イテレーション毎に学習結果出力
                  num_boost_round=1500,  # 最大イテレーション回数指定
                  early_stopping_rounds=100)

# テストデータの予測 (クラス1の予測確率(クラス1である確率)を返す)
y_pred_proba_train = model.predict(X_train)
y_pred_proba = model.predict(X_test)
accuracy_train = roc_auc_score(y_train, y_pred_proba_train)
accuracy_test = roc_auc_score(y_test,y_pred_proba)

print("  Accuracy_train = {}".format(accuracy_train))
print("  Accuracy_test = {}".format(accuracy_test))

結果は
予想外に時間がかかるのだ

Training until validation scores don't improve for 100 rounds
[50]    valid_0's binary_logloss: 0.440026
[100]   valid_0's binary_logloss: 0.434186
[150]   valid_0's binary_logloss: 0.431671
[200]   valid_0's binary_logloss: 0.430362
[250]   valid_0's binary_logloss: 0.429819
[300]   valid_0's binary_logloss: 0.42937
[350]   valid_0's binary_logloss: 0.429038
[400]   valid_0's binary_logloss: 0.428755
[450]   valid_0's binary_logloss: 0.428568
[500]   valid_0's binary_logloss: 0.428517
[550]   valid_0's binary_logloss: 0.428273
[600]   valid_0's binary_logloss: 0.428092
[650]   valid_0's binary_logloss: 0.427977
[700]   valid_0's binary_logloss: 0.427929
[750]   valid_0's binary_logloss: 0.427825
[800]   valid_0's binary_logloss: 0.427821
[850]   valid_0's binary_logloss: 0.427842
[900]   valid_0's binary_logloss: 0.427784
[950]   valid_0's binary_logloss: 0.427763
[1000]  valid_0's binary_logloss: 0.427756
[1050]  valid_0's binary_logloss: 0.42771
[1100]  valid_0's binary_logloss: 0.427678
[1150]  valid_0's binary_logloss: 0.427665
[1200]  valid_0's binary_logloss: 0.427658
Early stopping, best iteration is:
[1111]  valid_0's binary_logloss: 0.42764
  Accuracy_train = 0.8222143026699238
  Accuracy_test = 0.7791752907403984

この結果についてググって調べてみるのだ
今回はここまで:wave:

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

[Python]forとifの処理を1行で書く方法[リスト内包表記]

概要

  • リスト内包表記を使うことでfor・ifの処理を1行で書くことができます!
  • (ifのみの場合は、リスト内包表記を使えません)

コード(リスト内包表記を使う前)

  • リストの宣言やfor・ifを使うと4行書くことになります。
  • リスト内包表記を使うと…???
before_list = [1,2,3,4,5,6,7,8,9,10] # 1から10のリスト
# ----------リスト内包表記を使う前----------------------------------------
after_list = [] # 偶数だけを入れるリスト
for i in before_list:
    if i % 2 == 0: # 2で割り切れる(偶数)ときにリストへ追加する
        after_list.append(i)
# -----------------------------------------------------------------------
print(after_list)
  • 実行結果
[2, 4, 6, 8, 10]

コード(リスト内包表記を使った後)

  • リスト内包表記を使うと1行で簡潔に書くことができます!
before_list = [1,2,3,4,5,6,7,8,9,10] # 1から10のリスト
# ----------リスト内包表記を使った後---------------------------------------
after_list = [i for i in before_list if i % 2 == 0]
# -----------------------------------------------------------------------        
print(after_list)
  • 実行結果
[2, 4, 6, 8, 10]

まとめ

  • リスト内包表記を使うことでコードを1行で書くことができる。
  • リスト内包表記はいいぞ
  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

atmacup 08 初心者向け講座 from youtube

atmacupに初参戦することとなった
実質稼働日数は3日 & ド素人だがとりあえず出てみることにした
親切なことに主催者の方々が初心者向けの講座をyoutubeに挙げていたのでそれをまとめて理解しようと思う

まず初めにどういうデータを取り扱って何を予測してその結果がどのように評価されるのかについての理解を深める

一番初めにやることはコードを書くことではなくデータのスキーマと実際にどういうデータが入っているのかを確認すること。自分が使える情報は何か、それはどのような構造を持っているか実際にどういう値が入っているかをしっかりと把握する

現象を理解する データがどのように作られているかについて考える

実際に眺める エクセルで開いてデータを確認する。
例えばnullがめっちゃ起こっている→特定の場合には二つの情報が同時に欠落するようなオペレーションになっているかも?テキストのカラムでの表記ゆれなど。

pandas profilingというライブラリがありこれを使うとデータフレームの統計量を一気に確認することができる。簡単な使い方はProfileReport classに可視化したいデータフレームを渡してreportを作りto_fileを呼び出してhtmlファイルとして保存する方法です。
保存したhtmlをブラウザで開くと以下のような統計量をカラムごとに確認出来て便利
1,ユニークな値の数 2,頻度の高い値 3,意味のないカラムやNullが多いカラムかどうか...etc

pandas_profiling.ProfileReport(train_df)

学習用データと予測するデータの差分を知る
学習用データと予測データは別のデータなので、差分が強すぎるとよくない
差分をvenn図で見定める

c="Year_of_Release"
fig,ax=plt.subplots(figsize=(5,3))
venn2(
    subsets=(set(train_df[c].unique()),set(test_df[c].unique())),
    set_labels=("Train","Test"),
    ax=ax
)
ax.set_title(c)

year_of_Release以外のデータも見てみる

columns = test_df.columns
n_figs = len(columns)
n_cols = 4
n_rows = n_figs // n_cols + 1

fig, axes = plt.subplots(figsize=(n_cols * 3, n_rows * 3), ncols=n_cols, nrows=n_rows)

for c, ax in zip(columns, axes.ravel()):
    venn2(
        subsets=(set(train_df[c].unique()), set(test_df[c].unique())),
        set_labels=('Train', 'Test'),
        ax=ax
    )
    ax.set_title(c)

fig.tight_layout()

モデリングのための特徴量作成、アルゴリズムが理解しやすい特徴量を作ることを特徴量エンジニアリングという

連続変数の特徴量 連続変数の場合そのままモデルが解釈できるため基本的にはそのまま用いられる

欠損値の処理 通常欠損値があると欠損のことを考慮する必要があります。その場合平均値や中央値などで欠損を埋める必要がありますが、lightgbmなどのgbdtは欠損をそのまま自然に扱えるので基本的に気にする必要がない

連続変数とみなす場合にはその大小関係に意味があると思っていることと同値であることに注意する。大小関係に意味がない値の場合には連続変数ではなくてカテゴリ変数とみなして処理を行ったほうが適切な場合がある。勾配ブースティングは連続値とみなしてもある程度うまくやってくれるが線形モデルのように値の大小と予測結果が線形に比例するようなモデルでは予測性能が悪化する場合がある

def preprocess(input_df: pd.DataFrame) -> pd.DataFrame:
    output_df = input_df.copy()
    output_df['User_Score'] = input_df['User_Score'].replace('tbd', None).astype(float)
    return output_df

def create_continuous_features(input_df):
    input_df = preprocess(input_df)
    use_columns = [
        # 連続変数
        'Critic_Score',
        'Critic_Count', 
        'User_Score', 
        'User_Count',
        'Year_of_Release'
    ]
    return input_df[use_columns].copy()
assert len(create_continuous_features(train_df)) == len(train_df)
assert create_continuous_features(train_df.head()).equals(create_continuous_features(train_df.head()))

booleanの特徴量
特定の条件を満たしているかどうかを表す特徴量。例えばテキストカラムに特定の文字列が存在しているかどうかなどが該当する

def create_boolean_feature(input_df):
    output_df=pd.DataFrame()
    texts=[
        "japan","nintendo"
    ]
    for t in texts:
        output_df[f'Developer_has_{t}']=input_df["Developer"].fillna('').str.lower().str.contains(t).astype(int)
    return output_df
# https://github.com/nyk510/vivid/blob/master/vivid/utils.py
from contextlib import contextmanager
from time import time

@contextmanager
def timer(logger=None, format_str='{:.3f}[s]', prefix=None, suffix=None):
    if prefix: format_str = str(prefix) + format_str
    if suffix: format_str = format_str + str(suffix)
    start = time()
    yield
    d = time() - start
    out_str = format_str.format(d)
    if logger:
        logger.info(out_str)
    else:
        print(out_str)
from tqdm import tqdm

def to_feature(input_df):
    processors = [
        create_continuous_features,
        create_boolean_feature
    ]

    out_df = pd.DataFrame()

    for func in tqdm(processors, total=len(processors)):
        with timer(prefix='create ' + func.__name__ + ' '):
            _df = func(input_df)

        assert len(_df) == len(input_df), func.__name__
        out_df = pd.concat([out_df, _df], axis=1)

    return out_df
train_feat_df = to_feature(train_df)
test_feat_df = to_feature(test_df)

lightgbmによる学習
特徴量を作成したので次にlightgbmによる学習を行っていく交差検証という考え方が大事

cross validationとは 学習用データセットを複数に分割してそれぞれの分割で学習検証のデータセットを作りモデルの性能を見積もる枠組み。学習データの中で今の枠組みの性能を評価したい。

一番ナイーブな方法はkfold これは何も考えずにとにかくランダムに学習データを分割する。
その他にターゲット(目的変数)の分布が同じになるように分割するStratifiedと呼ばれる方法もある。groupkfolfは訓練データと検証データに同じグループが現れないように分ける方法

明日また更新します

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

ネットワーク分析2 実践編:karate club

Aidemy 2020/12/8

はじめに

 こんにちは、んがょぺです!バリバリの文系ですが、AIの可能性に興味を持ったのがきっかけで、AI特化型スクール「Aidemy」に通い、勉強しています。ここで得られた知識を皆さんと共有したいと思い、Qiitaでまとめています。以前のまとめ記事も多くの方に読んでいただけてとても嬉しいです。ありがとうございます!
 今回は、ネットワーク分析の2つ目の投稿になります。どうぞよろしくお願いします。

*本記事は「Aidemy」での学習内容を「自分の言葉で」まとめたものになります。表現の間違いや勘違いを含む可能性があります。ご了承ください。

今回学ぶこと
・karate clubについて
・中心性について
・グラフの分割

ネットワーク分析の実践

karate club

・このデータセットは、米国のある大学の空手クラブの交友関係を示したものである。これを使って、中心人物の特定コミュニティの抽出などの様々なネットワーク分析が行える。

スクリーンショット 2020-12-04 12.58.41.png

・Chapter1で扱ったnetworkxでは、「nx_karate_club_graph()」とするだけで、このデータセットを呼び出すことができる。また、この描画も今までにやった方法で行える。
・また、この頂点は「属性」を有している。ここでいう属性は、大学の2つのクラブ(club)を表している。例えば「G.node[18]['club']」とすると「Officer」と出力され、他の頂点についてもみると「Mr. Hi」と出力されるものもある。以下のコードでは、この2つの属性ごとに頂点を色分けしたものである。

スクリーンショット 2020-12-04 13.37.42.png

中心人物を見つける

・直感的な理解の通り、グラフの辺は、頂点同士に繋がりがあることを示す。社会ネットワークの例で言えば、人間関係のグラフにおいて、各頂点がその組織を構成する人間であるとすると、頂点AとBを結ぶ辺は、AさんとBさんの間に何らかのつながりがあることを示す。
・各頂点がネットワークの構造上どれぐらい有利な状態であると言えるかという指標を数値化したものが「中心性」である。人間関係の例で言えば、より多くの人と繋がりを持った方が有利であるといえ、また、多くのつながりを持った人と繋がることも有利であると言える。
・中心性の一つとして、「固有ベクトル中心性」がある。グラフGにおける固有中心性は「nx.eigenvector_centrality_numpy(G)」で取得できる。この値が大きいほど重要度が高いということになる。以下では、karate clubにおける固有ベクトル中心性を算出している。

・コードスクリーンショット 2020-12-04 14.14.10.png

・結果スクリーンショット 2020-12-04 14.14.21.png

媒介中心性

・ここでは、固有ベクトル中心性とは別の中心性を見ていく。また、ここでいう「経路」は全て最短経路を示すこととする。
・人間関係の例で、AさんとBさんは直接のつながりはないが、他の人物を経由することでつながりが見出せる場合を考える。すなわち、頂点AからBへの経路が存在する場合である。この場合、経路上に存在する頂点は、AさんからBさんへの情報伝達の役割を担っているといえ、この「情報伝達の役割を持った頂点」を有利な状態と考えて定義された中心性が「媒介中心性」と呼ばれるものである。この値が高いほど、より多くの情報伝達の役割を担っていると言える。
・グラフGにおける媒介中心性は「nx.communicability_betweenness_centrality(G)」で求められる。

・コードスクリーンショット 2020-12-04 14.29.51.png

・結果スクリーンショット 2020-12-04 14.30.00.png

最も中心性が高い頂点を求める

・karate clubで最も中心性の高い頂点を求める。求め方は、固有ベクトル中心性、あるいは媒介中心性をリストと見て、その中の最大値を取れば良い。具体的には、以下の通り。

スクリーンショット 2020-12-04 14.44.26.png

中心性の大きさに応じて頂点の大きさも変える

頂点の大きさを、中心性の大きさに応じて大きくするということを行う。頂点の大きさは、「nx.draw_networkx()」の引数で「node_size」を指定することで変更できる。以下では、それぞれの頂点の媒介中心性の値(size)に10000をかけることで、そのままnode_sizeとしている。

・コードスクリーンショット 2020-12-04 14.56.56.png

・結果スクリーンショット 2020-12-04 14.57.06.png

グラフの分割

連結成分

グラフをいくつかのグループに分けることを考える。この時のそれぞれのグループのことを「サブグループ」という。このようにグラフが分割される時、サブグループを連結成分ともいう。「nx.connected_component_subgraphs(G)」とすることで、それぞれの連結成分が返される。

クリーク

・前項でおこなったグラフの分割は、連結成分があることが前提であったため、これがない場合、すなわち全結合グラフである時は、グラフの密度に注目して分割していく。頂点の個数をn、辺の個数をmとすると、有向グラフの時は$\frac{m}{n(n-1)}$、無向グラフの時は$\frac{2m}{n(n-1)}$で算出することができる。
頂点が3つ以上あり、かつ密度が1(最大)であるサブグラフのことを「クリーク」という。以下のグラフは、頂点が3つ以上あり、全ての頂点が、他の全ての頂点と辺で繋がっているので、クリークであると言える。

スクリーンショット 2020-12-04 15.31.49.png

・また、以下では、グラフGと、そのグラフの任意の3つの頂点のリストを引数とし、その3つの点がクリークを形成するかを判定する関数「is_clique()」である。
・まずは3つの頂点をそれぞれペアにして、リストに格納する。渡された3点のうち2つをiとjとし、それぞれ同じものでなければ(i,j)としてnode_pair_listに格納する。
・そらにこのそれぞれのペアについて、(頂点A,頂点B)あるいは(頂点B,頂点A)のどちらかのペアが、グラフの辺「G.edges()」の中にあれば、「AとBは繋がっている」ということができ、全てのペアについて「繋がっている」判定がなされた時は、クリークが形成されるとしてTrueを返す。

スクリーンショット 2020-12-04 15.44.26.png

n-クリーク

・前項でクリークについて学んだが、実際のグラフでクリークが形成されることはあまりない。クリークの条件は、言い換えれば「サブグラフに含まれる全ての頂点間の距離」が「1」であるものを指すので、この条件を緩くして、「距離の最大値がn以下である」ものを「n-クリーク」と呼ぶ。
・以下は、頂点のリストと「n」、最短距離が与えられた時、与えられた頂点がn-クリークを形成するかを判定する関数「is_n_clique()」である。判定するには、距離の最大値を求めれば良いので、各頂点i,jについての距離dist[i][j]の最大値がn以下だったらTrueを返すようにすれば良い。

スクリーンショット 2020-12-04 15.59.12.png

クラスタリング

・ネットワーク分析におけるクラスタリングは、グラフ理論を用いた分析手法を指す。具体的には、連結されたグラフをいくつかのサブグラフに分割する操作のことクラスタリングという。また、この時のサブグラフをコミュニティという。
・ネットワークにおけるクラスタリングのメリットとして、大規模なネットワーク構造でも、自動でノードの分類ができることが挙げられる。また、具体的な活用例として、SNSの個々のアカウントをノード(頂点)としてクラスタリングを行うことで、個人に最も適した広告やレコメンドなどを行うことが可能になる。
・networkxでは、communityモジュールのbest_partition()を利用することでクラスタリングが行える。以下ではそれを実行している。各ノードと、何番目のコミュニティに属するかが辞書型で渡されている。

・コードスクリーンショット 2020-12-04 23.44.11.png

コミュニティごとに色分けする

・前項で算出した「partition」の値は、属するコミュニティを示すので、partitionの値が同じもの同士で色分けすることで、属するコミュニティをわかりやすく可視化する。以下ではこれをおこなっているが、各色について、引数の「cmap」では「plt.cm.RdYlBu」で指定している。これは、ランダムに色を変化させるもので、振り分けられる色は発散的になっている。今回は頂点番号が1桁のものと2桁のもので色分けしている。その場合、node_colorにはこれを満たす条件式を指定する。

・コードスクリーンショット 2020-12-05 11.43.40.png

・結果スクリーンショット 2020-12-05 11.43.58.png

モジュラリティ指標

・ここでは、「良い分割」を評価する手法を紹介する。「良い分割」とは、各コミュニティ内の結びつきが強くコミュニティ同士の結びつきは弱いものとなるような分割である。この指標となるのが「モジュラリティ指標」である。これは「community.modularity(partition,G)」で求められる。この値が高いほど良い分割であると言える。
・以下では、ここまでの分割手法「community.best_partitionによる分割」と「頂点番号による分割」のモジュラリティ指標を算出している。

スクリーンショット 2020-12-05 11.53.07.png

まとめ

・データセット「karateclub」を呼び出すには「nx.karate_club_graph()」とすれば良い。
・ネットワーク上で頂点がどれくらい有利であるかを示すのが中心性である。人間関係では、より多くのつながりを持っていたり、つながりを持った人とつながることで中心性が高いとみなされる。中心性は「固有ベクトル中心性」「媒介中心性」がある。
・グラフをいくつかのグループに分けるとき、分けられたそれぞれのグループをサブグループ(連結成分)という。また、全ての頂点がつながっているサブグラフをクリークという。また、このサブグラフに分割することをクラスタリングということもあり、この時のサブグラフをコミュニティという。また、このコミュニティの分割を評価する指標となるのが「モジュラリティ指標」である。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

[python] リストを数字に変換する

今回はPythonで競技プログラミング用に備忘録として書いていこうと思います〜

int型のリストをもらった時にint型に変換したい時のために書きました!

[1,2,3] ====>>>> 123

実際のコード

sample.py
sum_number = int("".join([str(n) for n in digits]))

説明

"".join()

でList型をStr型に変換します

int()

最後にint型に変換します

誰かのためになれば幸いです!

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

AIにクラッピーを見せると何を思うのか(Azureの画像分析APIです)

はじめに

#クラッピーチャレンジ Advent Calendar 2020 8日目の記事です。
クラッピーが手元にないので、過去に撮影したクラッピーの画像をAzure Computer Visionの画像分析APIに投げてみました。画像分析APIを用いると、画像の説明文を生成してくれます。
HoloLens2でもこんな記事「HoloLens2 × Azure Cognitive Services(画像分析APIで画像説明文生成)」を書いているので見てみてください。

開発環境

  • Windows 10
  • Python 3.6
  • Azure

導入

画像分析APIに投げるPythonプログラムは下記になります。AzureポータルからComputer Vision APIを作成し、エンドポイントとキーをコピペしてください。

analyze.py
import os
import sys
import requests
# If you are using a Jupyter notebook, uncomment the following line.
# %matplotlib inline
import matplotlib.pyplot as plt
from PIL import Image
from io import BytesIO
import os

subscription_key = "<Insert Your API Key>"
endpoint = "https://<Insert Your Resource Name>.cognitiveservices.azure.com/"
analyze_url = endpoint + "vision/v3.1/analyze"

# Set image_path to the local path of an image that you want to analyze.
# Sample images are here, if needed:
# https://github.com/Azure-Samples/cognitive-services-sample-data-files/tree/master/ComputerVision/Images

dirname = "clappy/"
files = os.listdir(dirname)

for file in files:
    # Read the image into a byte array
    image_data = open(dirname + file, "rb").read()
    headers = {'Ocp-Apim-Subscription-Key': subscription_key,
            'Content-Type': 'application/octet-stream'}
    params = {'visualFeatures': 'Categories,Description,Color'}
    response = requests.post(
        analyze_url, headers=headers, params=params, data=image_data)
    response.raise_for_status()

    # The 'analysis' object contains various fields that describe the image. The most
    # relevant caption for the image is obtained from the 'description' property.
    analysis = response.json()
    print(analysis)
    image_caption = analysis["description"]["captions"][0]["text"].capitalize()

    # Display the image and overlay it with the caption.
    image = Image.open(BytesIO(image_data))
    fig = plt.figure()
    plt.imshow(image)
    plt.axis("off")
    _ = plt.title(image_caption, size="x-large", y=-0.1)
    # plt.show()
    fig.savefig("result/"+file)

実行

19030338_1491492647538744_773272452102692181_n.jpg
19055112_1491486954205980_1297198113259920897_o.jpg
19264650_1499046780116664_3424860175366470380_o.jpg
23318983_1635601143127893_871019205316474673_n.jpg
23319335_1635695543118453_5621071620397277174_n.jpg
23376365_1636484413039566_2707645763496136760_n.jpg
23376647_1635572793130728_923685704748732178_n.jpg
23473031_1635696433118364_7171265576451749495_n.jpg
25446047_1673387339349273_34401910827914340_n.jpg
30710917_1788011027886903_6576225328295837696_n.jpg
37179064_1897943366893668_7397617731798827008_n.jpg
43250600_2016093018412035_3469528603807449088_n.jpg
16665649_1374984332522910_4747950487955808095_o.jpg
16797175_1386987047989305_7491875210929899417_o.jpg
17349923_1413753251979351_164391718457275837_o.jpg
17358723_1413887895299220_7035955493849861248_o.jpg
17545527_1421129511241725_2620203238592672616_o.jpg
18952605_1491492500872092_9196702837233785555_n.jpg

まとめ

割と正確にtoyだった。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Python 入門「その1」:Hello World!

前書き

現在,Python を学習しています.
私は学習した知識をまとめて誰かに教えるのが好きで,そうするとよく覚えられるので,私の学習のためにも,何回かに分けて Python の入門を Qiita で投稿していきたいと思います.プログラミング未経験者の方にも,できるだけわかりやすいよう心がけていきます.とはいえ,私も Python に関しては初心者ですので,間違いもあるかもしれません.その時は,指摘していただけると幸いです.

開発環境としては,Windows 10 で,Jupyter Notebook を使って学習しています.Notebook でレポートのようにまとめているのですが,Qiita には,Notebook にまとめた内容をさらにわかりやすくし書いていきたいと思います.Python 3 となりますので,ご注意ください.

なお,秀和システム社の「Python プログラミングパーフェクトマスター」という本で学習しています.当たり前ですが,本の内容をパクったりしているわけではなく,本で勉強し,その応用として Qiita に独自にまとめています.そのため,このシリーズでは,簡潔に解説しています.詳しい解説や説明,開発環境の用意や使い方についてはその本をご覧ください.

Hello World!

さて,さっそく Hello World してみましょう.
プログラミングをやったことない方に説明すると,Hello World とは初めてプログラミングするときなどの恒例行事(?)で,「Hello World!」と表示するものです.

以下のコードを実行してみてください.この時,コピペするのではなく,自分で手打ちで写すと覚えられると思います.

print ('Hello World!');

「Hello World!」と表示されたと思います.

  • print は,カッコ内の文字列などを表示する関数(次回説明)です.
  • シングルクォーテーション(')で囲むと,囲まれた範囲内が文字列として認識されます.
  • ダブルクォーテーション(")でも,同様の動作をします.
  • 文末のセミコロンは,なくても構いません.Python では不要ですが,ほとんどの言語でセミコロンをつける必要があるので,書く癖をつけておくといいと思います.Python ではあってもなくても同じです.
  • print とカッコの間にスペースがありますが,なくても構いません.見やすさのためにこうしていますが,見やすいほうで構いません.

このシリーズでは,セミコロンあり・空白ありで解説していきます.

では,次はこれを実行してみましょう.

print (Hello World!);

この例では,シングルクォーテーションやダブルクォーテーションで囲っていないため,エラーが出ます.文字列を扱うときは,必ず囲むようにしましょう

print ("Hello World!");

先ほどはシングルクォーテーションを使いましたが,この例ではダブルクォーテーションを使っています.シングルクォーテーションのときと変わらず,「Hello World!」と出力されたはずです.

シングルクォーテーションとダブルクォーテーション

シングルクォーテーションとダブルクォーテーション,どちらを使っても同じと解説しました.では,どちらを使うべきなのでしょうか.基本的には,好みで構いません.シングルクォーテーションのほうがすっきりしてて見やすいか,ダブルクォーテーションのほうがわかりやすいか,どちらでも大丈夫です.しかし,以下の例を試してみてください.

print ('Airi's computer');
print ("She said, "programming is fun".");

英語の勉強ではないので,文が間違えていても気にしないでください()
上は「あいりのコンピュータ」下は「彼女は「プログラミングは楽しい」と言いました」という意味です.

どちらも,エラーが出たと思います.
よく見てみると,文字列中にシングルクォーテーション・ダブルクォーテーションが入っているのがわかるでしょうか.この時,その地点で文字列の終わりと判断されてしまいます.このようなときに,シングルクォーテーションとダブルクォーテーションを使い分けるといいと思います.

print ("Airi's computer");
print ('She said, "programming is fun".');

また,これ以外に,\を使う方法もあります.詳しくはまた今度説明しますが,\'\" とすることで,エラーは出ず,正しく出力されるようになります.

print ('Airi\'s computer');
print ("She said, \"programming is fun\".");

トリプルクォーテーション

シングル,ダブルときたらトリプルです.シングルクォーテーションもしくはダブルクォーテーションを 3 つ使うことで,トリプルクォーテーションになります.これは,文字列中に改行があっても,文字列の続きとして扱われるようになります.

print ('''りんご
いちご
みかん''');

トリプルクォーテーションを使わなくても,\n と入れることで,改行することができます.これも詳しくはまた今度説明します.

print ('りんご\nいちご\nみかん');

終わりに

Hello World から,文字を表示する方法までわかったでしょうか.
Qiita で本格的に記事を書くのは初めてなので,間違いや,見づらいところがあるかもしれません.
毎日勉強を進めているので,毎日更新できたらなと思いますが,それは厳しいかもしれません.
それではまた次回.

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

Qiitaに投稿した記事をWordPressにも自動投稿してみた

はじめに

普段はpythonを使ってエンジニアしたり、PHPでWeb系のお仕事したりと広く浅くプログラミングやってます。

pythonでこんなことできるよみたいな紹介ができればと思います。

スクレイピングとかも面白いかなと思ったのですが、そのあたりの記事はたくさんあるしなぁとか悩んだ結果、せっかくだしQiitaAPIを触って何かしてみようと思いました。

なので今回は、QiitaAPIを使ってQiitaに自分が投稿した記事を自分のサイト(WordPress)にも投稿できるよというのを紹介したいと思います。

python始めた方向けにpythonこんなことできるよという雰囲気を味わってもらえたらいいなと思っているので、細かい部分は省きながら説明していきます。

準備が必要なもの

今回の自動化で事前に準備が必要なのは、

1、QiitaAPIのアクセストークン
2、自分のWordPressのAPI用アカウント情報

それぞれの準備については今回の目的と外れますので参考になるサイトをそれぞれ張っておきます。

1、こちら
2、こちら

サンプルコード

今回は解説ではなくてこんなことできるよの紹介記事になるのでサンプルコードと実際にこの記事を自分のサイトに自動投稿してみた結果をご紹介します。

import requests
import json


def main():
    qiita_token = "QiitaAPIのアクセストークン"
    headers = {
            'Content-Type': 'application/json',
            'Charset': 'utf-8',
            'Authorization': 'Bearer {}'.format(qiita_token)
        }
    pf = PostFunctions()
    r_get = requests.get("https://qiita.com/api/v2/users/QiitaのユーザーID/items/", headers=headers)
    results = r_get.json()
    for result in results:
        article_id = result["id"]
        if pf.postcheck(article_id):
            pass
        else:
            title = result["title"]
            content = result["rendered_body"]
            pf.postadd(article_id)
            Article2WP(title, content)


class PostFunctions():
    # 既にWPに投稿されているか確認
    def postcheck(self, article_id):
        with open("article_id.txt", "r", encoding="utf8") as f:
            id_list = [s.strip() for s in f.readlines()]
            if article_id in id_list:
                return True
            else:
                return False
    # WPに投稿した記事IDを追加
    def postadd(self, article_id):
        with open("article_id.txt", "a", encoding="utf8") as f:
            print(article_id, file=f)

def Article2WP(title, content):
    user_id = "WP投稿用アカウントID"
    password = "WP投稿用アカウントパスワード"
    end_point_url ="https://自分のサイトドメイン/wp-json/wp/v2/posts/"
    payload = {
        'title': title,
        'content': content,
        'status': "draft",
        'comment_status': 'closed',
        'featured_media': アイキャッチ画像ID,
        'categories': カテゴリーID,
    }
    headers = {'content-type': "Application/json"}
    r = requests.post(end_point_url, data=json.dumps(
        payload), headers=headers, auth=(user_id, password))

基本的なコードはこれだけです。

定期実行するように作ってありますので既に自動投稿した記事をテキストファイルにメモしておくようにしていますが、そのあたり不要な方は削ってしまってもきちんと動きます。

最後に

個人的にはもう少し時間かかるかなと思ったのですが、かなり便利なAPIが双方に用意されていたのでただAPI叩いただけになってしまいました。

この記事を実際に自分のサイトに自動投稿してみましたのでどんな感じに投稿されたか気になる方は見てみてください。

こちら

Qiitaと自分のサイトに投稿したい需要があるのか分からないですが、誰かの参考になったら嬉しいです。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

ネットワーク分析1 基礎知識

Aidemy 2020/12/8

はじめに

 こんにちは、んがょぺです!バリバリの文系ですが、AIの可能性に興味を持ったのがきっかけで、AI特化型スクール「Aidemy」に通い、勉強しています。ここで得られた知識を皆さんと共有したいと思い、Qiitaでまとめています。以前のまとめ記事も多くの方に読んでいただけてとても嬉しいです。ありがとうございます!
 今回は、ネットワーク分析の一つ目の投稿になります。どうぞよろしくお願いします。

*本記事は「Aidemy」での学習内容を「自分の言葉で」まとめたものになります。表現の間違いや勘違いを含む可能性があります。ご了承ください。

今回学ぶこと
・ネットワーク分析について
・基本用語
・グラフの作成

ネットワーク分析について

ネットワーク分析とは

ネットワーク分析とは、グラフを用いてさまざまな事象と事象の関係について調べるという分析手法である。ここでいう「ネットワーク」とは、いわゆるコンピューターネットワークに限らず、広く「複数の点がなんらかの関係で繋がっている(線になっている)」物を指す。例えば、人と人をそれぞれ点とすれば、それを繋げば「人間関係」となる。ちなみに、友人関係、血縁関係といった人間関係は「社会ネットワーク」と呼ばれる。このネットワークを可視化したものが「グラフ」となる。このグラフも、一般的な意味とは必ずしも一致しないため、注意する。
・ネットワーク分析の応用例として、ネットワークの中心性を分析して影響力のある人物を特定する「ソーシャルネットワーク」、企業間の取引データから業種の分類などを行う「企業間ネットワーク」などがある。
・今回は、networkxというネットワーク分析用ライブラリを使い、ネットワークの「クラスタリング」「中心的なノードの発見法」を学んでいく。Chapter1ではこのライブラリの関数の説明を行う。Chapter2ではkarateclubというデータセットを使って、実際にデータ分析を行う。

グラフの描画例

・Chapter1では、networkxライブラリの使い方の中でも、特にグラフの記述についてを学んでいく。グラフは点と線からなる関係構造を可視化したもの、ということができる。
・グラフの作成はこのnetworkxで行うが、表示については今までと同様に「plt.show()」で行う。以下は実際にグラフを作成し、表示したものである。(コードについては後述する。)

スクリーンショット 2020-12-03 20.22.29.png

グラフ理論の基礎

基本用語

・前項のような、グラフを数学的に扱う分野「グラフ理論」という。
・基本用語について、グラフを構成する各点を「頂点(ノード)」という。そして、この頂点同士を結んだ線を「辺(エッジ)」という。
・ベクトルのように「矢印の向き」がある辺を含むグラフのことを「有向グラフ」といい、全ての辺に向きがないグラフを「無向グラフ」という。
・そして、この点と辺を連続で辿ったものを、文字通り「道(path)」という。有向グラフの場合は、辺の向きが進む方向と一致していなければ、道とはならない。また、一般的に「道」というときは、同じ頂点は二度と通らない物を指す。
スクリーンショット 2020-12-03 20.43.05.png

その他の用語(応用)

・「辺」に数字が書いてあるものを、「重み付きグラフ」といい、その数字のことを「重み」という。ニューラルネットワークでよく見られる。また、ある頂点から頂点に移動する道を考えるとき、通過する重みの合計が最小になるもの「最短経路」という。
・道のうち、循環になっている(一回りすると最初の頂点に戻る)物を「閉路」という。
スクリーンショット 2020-12-03 20.53.31.png

・また、閉路を持たないグラフのことを「木」という。そして、木の一番上の頂点を「根」、一番下の頂点を「葉」という。自分より上位にある頂点を「親」、下位にある頂点を「子」という。
・木の特徴として、辺の数の合計が「頂点の数-1」で求められることが挙げられる。
スクリーンショット 2020-12-03 21.00.14.png

・他のグラフとしては、「二部グラフ」というものがある。これは、頂点を2グループに分けて、各グループ同士では辺をつくらず、違うグループ同士でのみ辺を作るものを指す。以下の図では、赤いグループと青いグループで分け、同じグループ内では辺が作られていないことがわかる。
スクリーンショット 2020-12-03 22.02.40.png

用語(発展)

・有向グラフにおいて、方向を無視しないでどの頂点との間にも道が存在するものを「強連結」という。一方で、方向を無視した場合(無効グラフとして見た場合)に任意の二点間に道が存在するものを「弱連結」という。

スクリーンショット 2020-12-04 0.20.06.png

・次に、隣接行列について見ていく。以下は、上記グラフのうちの一つ目(強連結の有向グラフ)の隣接行列である。スクリーンショット 2020-12-04 0.25.38.png

・それぞれの行についてみてみる。まず一行目については2列目と5列目が「1」となっている。これは、頂点「1」から頂点「2」と「5」への辺が存在していることを表す。逆に「0」の部分は頂点1からの辺が存在していないことを表す。また、自分自身(頂点「1」)への辺は存在していないものとして扱う。それ以降の行についても同様に考えることができる。
・そして、以下がこのグラフの隣接リストと呼ばれるものである。スクリーンショット 2020-12-04 0.29.53.png

・これは、i番目のリストには、頂点iから辺が出ている頂点が格納されている。基本的には隣接行列から得られる情報と変わりはないが、それぞれ以下のような特徴がある。
 ・隣接行列の方が計算で扱いやすい。
 ・隣接リストの方が消費メモリが少ない。
 ・頂点iからどの頂点への辺が出ているかを確認するとき、隣接リストの方が参照数は少ない。

・最後に、ある頂点につながっている辺の数を「次数」という。特に有向グラフについて、その頂点に入ってくるような辺の数を「入次数」、出ていくような辺の数を「出次数」という。

networkxについて

グラフの描画の実装

・以下は、最初に作成したグラフの描画のコード例である。このコードについて、詳しく説明していく。
スクリーンショット 2020-12-04 0.40.18.png

・まず、「G=nx.path_graph(5)」について、これは「nx.path_graph()」によって「G」という無向グラフを作成していることを示す。引数で、頂点の数を指定する。頂点は後から追加できるので、この値は厳密に決める必要はない。今回の例で言えば、5と指定しているので、それぞれラベルが「0~4」の5つの頂点が一繋ぎの辺となって作成されている。ちなみに、「nx.Graph()」とすることで、まっさらなグラフを作成できる。
「G.add_edge()」で、指定した頂点を結ぶ辺を追加している。この段階で存在していない頂点が指定された場合は自動的に追加される
「pos = nx.spring_layout(G)」で、Fruchterman-Reingold force-directed algorithmという力学モデルの描画アルゴリズムで、頂点を描画する位置を決めている。この部分はなくても構わない。
「nx.draw_networkx(G)」で、作成したグラフGを描画する。これを「plt.show()」とすることで、表示できる

有向グラフの描画

有向グラフについては、グラフの作成時に「nx.DiGraph()」とすれば良い。あとは基本的には同じ。
・有向グラフに限らず、まとめて辺を作成したいときは「G.add_edges_from()」に、辺のリストを渡せば良い。有向グラフの場合、辺の指定は(入力点,出力点)とする。例えば(0,1)なら頂点0から頂点1に向かう辺が作成され、(1,0)なら頂点1から頂点0に向かう辺が作成される。このように指定順で意味が異なるので注意する。
また、頂点をまとめて追加したい時は「G.add_nodes_from()」に、同じように頂点のリストをまとめれば良い。

スクリーンショット 2020-12-04 1.13.35.png

頂点と辺の描画

「draw_networkx_nodes()」で頂点を、「draw_networkx_edges()」で辺を、それぞれ色や大きさを指定して描画することができる。
・各引数について、「G」「pos」はこれまでと同様に指定する。ちなみに、ここでのposは必須となる。
「nodelist」で、描画したい頂点をリストで渡して指定する。「node_color」で頂点の色を指定でき、「node_size」で大きさを指定できる。この3つについては「node」の部分を「edge」とすれば辺についても同じように指定できる。
・頂点について、「alpha」で透過度を指定できる。0〜1の値を取る。また、辺については「width」で太さを指定できる。

・コードスクリーンショット 2020-12-04 1.28.45.png

・結果スクリーンショット 2020-12-04 1.29.02.png

グラフに重みを加える

・グラフGに、「add_weighted_edges_from()」というメソッドを使うと、重みを付け加えることができる。引数には「(頂点1,頂点2,重み)」をまとめたリストを渡す。
スクリーンショット 2020-12-04 10.04.48.png

・これを使って作成した以上のグラフについて、networkxのメソッドを試していく。
「nx.info(G)」とすると、グラフGの基本情報が出力される。
スクリーンショット 2020-12-04 10.06.35.png

Typeはグラフの種類、Number of nodes(edges)は頂点(辺)の数、Average degreeは全ての頂点から出ている辺の平均数である。
「G.nodes()」というようにすると、頂点をリスト型で取得できる。同様に、「G.edge()」とすると辺を構成する頂点をリスト型で取得できる。
・また、「G.degree()」とすると、すべての点における次数(隣接する点の数)が取得できる。
「nx.adjacency_matrix(G)」とすることで、Gの隣接行列が取得できる。ただし、返り値はリストではないので、そのままでは扱えない。以下では、重み付きグラフGについて、隣接行列を取得し、その各要素(頂点)について、i行j列目(頂点iとj)の辺の重みが0でなければ道が存在するということになるので、その部分の辺について、重みを取得して、隣接リストadj_listに格納している。

・コード(Gの作成〜隣接行列の取得まで)スクリーンショット 2020-12-04 11.25.21.png

・結果スクリーンショット 2020-12-04 11.25.34.png

・これを見ると、「adj_matrix」には「頂点1,頂点2,重み」というような隣接行列が作成されていることがわかる。

・コード(隣接リストの作成と表示)スクリーンショット 2020-12-04 11.28.42.png

・結果スクリーンショット 2020-12-04 11.28.53.png

「adj_matrix[(i, j)]」で、頂点(i, j)間の辺の重みが取得できる。辺が存在しなければ0を返す、それ以外の値であれば、その値を重みとする辺が存在するということになる。よって、隣接リスト「adj_list」の「i」番目について、(j, adj_matrix[(i,j)])として格納することで、例えばadj_list[0]の値が(1,5)であれば、「頂点0から頂点1の辺の重みは5である」という意味を持つリストが作成できる。

networkxのその他のメソッド(応用)

・他のメソッドとしては、グラフを連結成分で分割する「nx.connected_component_subgraphs(G)」がある。「連結成分」については、次のChapterで説明する。以下では、「0~4の線」「5~6の線」を連結し、これで分割したものをリスト化して、前者はG[0]、後者はG[1]と、分けて描画している。
スクリーンショット 2020-12-04 11.55.08.png

・結果スクリーンショット 2020-12-04 11.55.22.png

・次のメソッドは「nx.all_pairs_dijksta_path(G)」で、これは、「ダイクストラ法」というもので、頂点間の最短経路を求めることができる。最短経路は「重みが最も小さくなるような道」のことである。よって、前に行った方法で重み付きグラフを作成し、それをGとしてこのメソッドを使い、辞書型で取得すると、{i:{j:[最短経路の時に通る頂点のリスト]}}という辞書で、iからjまでの最短経路を取得できる。以下では、この最短経路の距離(重みの合計)も算出している。
・距離の求め方について、辺の重みは「edge_labels[i,j]」で求められる。これを応用して、距離を足し合わせている。

スクリーンショット 2020-12-04 12.18.14.png

・結果(一部のみ)スクリーンショット 2020-12-04 12.19.01.png

まとめ

ネットワーク分析とは、グラフを使ってさまざまな事象の関係について調べる分析手法である。この時のグラフ上でのそれぞれの事象を頂点(node)、繋がりを辺(edge)という。
・ネットワーク分析では、networkxというライブラリを使用する。グラフもこれを使って作成する。無向グラフ「nx.path_graph()」で作成できる。有向グラフ「nx.DiGraph()」で作成でき、「nx.Graph()」は条件のないグラフが作成される。
・このグラフに頂点や辺を追加していくことでグラフを作成していく。作成したグラフの描画は「nx.draw_networkx(G)」のように行う。
・グラフの辺に重みを付け加える時は「add_weighted_edges_from()」を使うことで行える。また、この重みについて、隣接行列隣接リストに変換することで、計算などにも使えるようになる。

今回は以上です。最後まで読んでいただき、ありがとうございました。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む

pyenv-virtualenv 環境で jupyter notebook を開きたい

概要

とあるOSSライブラリを使ってみようとしていたとき、
jupyter notebookの形でexample コードが実装されていた。

まずはいつも自分が使っているpyenv-virtualenvの環境から
単純にjupyter notebook を起動。

jupyter notebook

としてブラウザを開いて、

import ****

とすると、modulue not found

と怒られた。

pyenv-virtualenv 環境から jupyter notebook を開く

この状況では、
https://qiita.com/shierote/items/421af3c741c0f9956b09
ここを参考にすることで解決した。ありがとうございます。

  • このエントリーをはてなブックマークに追加
  • Qiitaで続きを読む