PythonのWaveモジュールを使ってwavファイルを編集するの記事の以下のコードをコピペしたら、エラーが出てしまいました。

pywave.py

import wave
import struct
from scipy import fromstring, int16

# ファイルを読み出し
wavf = '/data/input/test.wav'
wr = wave.open(wavf, 'r')

# waveファイルが持つ性質を取得
ch = wr.getnchannels()
width = wr.getsampwidth()
fr = wr.getframerate()
fn = wr.getnframes()

print("Channel: ", ch)
print("Sample width: ", width)
print("Frame Rate: ", fr)
print("Frame num: ", fn)
print("Params: ", wr.getparams())
print("Total time: ", 1.0 * fn / fr)

# waveの実データを取得し、数値化
data = wr.readframes(wr.getnframes())
wr.close()
X = fromstring(data, dtype=int16)

threshold.py:25: DeprecationWarning: The binary mode of fromstring is deprecated, as it behaves surprisingly on unicode inputs. Use frombuffer instead
  X = fromstring(data, dtype=int16)

fromstringはunicode入力のときに何かすごい挙動するからやめて、みたいな事のようです。そして、frombufferを使え、と。
調べてみると、frombufferはnumpyのメソッド(?)のようなので、importにnumpyを入れて、fromstringはnp.frombufferにします。

pywave.py

import numpy as np
import wave
import struct
from scipy import fromstring, int16

# ファイルを読み出し
wavf = 'test.wav'
wr = wave.open(wavf, 'r')

# waveファイルが持つ性質を取得
ch = wr.getnchannels()
width = wr.getsampwidth()
fr = wr.getframerate()
fn = wr.getnframes()

print("Channel: ", ch)
print("Sample width: ", width)
print("Frame Rate: ", fr)
print("Frame num: ", fn)
print("Params: ", wr.getparams())
print("Total time: ", 1.0 * fn / fr)

# waveの実データを取得し、数値化
data = wr.readframes(wr.getnframes())
wr.close()
X = np.frombuffer(data, dtype=int16)

これで無事動きました。

こんにちは。ひろちょんです。

機械学習を始めたばかりの初心者が、様々な記事を参考にしつつアヤメ分類をしてみます。
説明を加えながら書いていくので、是非参考にしてみてください！

↓記事を読む前に筆者について知りたい方へ↓
>>>詳しいプロフィールはコチラ

～目次です～
1. 対象読者とか書いてみる
2. 《Google colaboratory》で進めていきます
3. 色々とimportしていく
　　・オートパイロット状態でimportする奴
　　・scikit-learnの様々な機能をimport!!
4. アヤメのデータセットをじっくり見ていく
　　・インスタンスを生成
　　・アヤメデータはどうなってるの？
5. Pandasを用いてデータを分かりやすくする
　　・配列をDataFrameに変換しよう！
　　・実際にDataFrameを見ていこう！
6. データセットを分割する
　　・説明変数と目的変数とは？
　　・学習用とテスト用のデータはどうするの？
　　・train_test_split関数を使う！
7. データをmatplotlibで可視化する
　　・データを可視化する意味は？《特徴量選択》
　　・散布図をplt.scatterで描いていく！
8. 機械学習アルゴリズムを使っていこう！
　　・特徴量を選択しよう！
　　・モデルを構築⇒学習⇒予測させる
　　・モデルが予想したデータの答え合わせ
　　・2つの結果の違いについて詳しく見る

1. 対象読者とか書いてみる

今回は機械学習でアヤメの分類を行っていくことが主題なので、レベルとしてはこんな感じ。

• 何かしらのNotebook形式でプログラムを実行できる人
• Python初心者以上
• 機械学習の流れが何となくわかる人

2. Google Colaboratoryで進めていきます

今回はバージョンとか、仮想環境とか…諸々のしがらみに囚われたくなかったので、Google Colaboratoryにてアヤメ分類をしていきたいと思います。

始め方はすごく簡単で、Google Colaboratoryのサイト左上の《ファイル》から《Python3の新しいノートブック》をクリックしてください。

すると新しいノートブックが生成されるので、左上のUntitled0.ipynbから名前を変更してあげてください！

※僕は名前を『初めてのアヤメ分類』にしておきました。笑

初学者の方はページを見つつ、一から作成することをオススメしますが、一応完成形も↓のGitHubにてノートブックを公開しています。
>>>GitHubへはコチラから

3. 色々importしていく

オートパイロット状態でimportする奴

とりあえず以下をインポートしていきますね。(矢印先は僕のライブラリへの感想です。)

• Numpy　⇒　計算機能凄めライブラリ
• Pandas　⇒　表を扱いやすいライブラリ
• Matplotlib　⇒　グラフ可視化ライブラリ
• warnings　⇒　たまに邪魔な警告を消すライブラリ

そしてプログラム化するとこんな感じ。

importするライブラリ.ipynb

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')

scikit-learnの色々な機能をimport!!

機械学習ライブラリのscikit-learnには分類を行う以外にも様々な機能が存在しています。今回使用するものをサラっと紹介しますね。(後で深堀りしていきます)

1つ目はscikit-learnに用意されているデータセットです。
scikit-learnのデータセット集一つにアヤメのデータセットがあります。
>>>詳しくはコチラ(UCI Machine Learning Repository: Iris Data Set)

２つ目が学習用とテスト用にデータセットを分ける機能です。
分ける理由としては、《学習用のデータで予測するモデルを作成する》⇒《テスト用データで作成したモデルは良い物か判定する》という機械学習では鉄板の2つの処理を行いたいからです。

3つ目が分類を行うアルゴリズムになります。機械学習には様々な分類のアルゴリズムが存在して、以下のプログラムではアルゴリズムとして『線形のSVM』を用いるので、LinearSVCをインポートしています。
>>>SVM(サポートベクターマシーン)について詳しく知りたい方はコチラ(Qiita)

では上記で紹介した3つをimportしていきます。

scikit-learnの様々な機能をimport.ipynb

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.svm import LinearSVC

4. アヤメのデータセットをじっくり見ていく

まず今回の目的をハッキリさせておきましょう！

どのサイトでも分類の部分を重要視しすぎて、結局何がしたかったのか見失いがちです。

今回の目的は…

『アヤメのがく片や花びら幅や長さの数値』を用いて『アヤメ属(花)の種類を分類する』こと

そこで既にアヤメの種類によっての幅や長さのデータを数値を集めてくれているデータセット(先ほどimportした奴)を使って、機械学習により予測していきます。

フィッシャーのアヤメのデータから散布図行列を描く(https://teenaka.at.webry.info/201803/article_12.html )から引用

ちなみに上の図がアヤメになります。この花の長さや幅を測定して数値化したデータセットをimportしました(/・ω・)/

インスタンスを生成

今はアヤメのデータセット機能をインポートしただけなので、まずインスタンスを生成していきます。

アヤメインスタンスを生成する.ipynb

heacet = load_iris()

インスタンスの名前に特に深い意味はないです。単に自サイトの宣伝をしているだけです！笑

もし嫌だという方はコピーするだけではなく、自分で書いていくと理解も深くなるので、ここのインスタンス名を変えて実行してみてください！

アヤメデータはどうなってるの？

さっそくデータセットの中身を見ていきましょう。
↓コチラを実行してください↓

アヤメデータの中身を確認する.ipynb

print("与えられたデータ")
print(heacet.data)
print(heacet.data.shape)
print("-----------------")
print("予測するデータ")
print(heacet.target)
print(heacet.target.shape)
print(heacet.target_names)

実行して出力を見ると、どうやら《数値が書かれた150×4の2次元配列》と《0,1,2と書かれた150×1の1次元配列》が得られたみたいですね。

また0,1,2はそれぞれ『setosa』,『versicolor』,『virginica』対応していることも分かります。

よって与えられた《数値が書かれた150×4の2次元配列》を《データセットを学習用と予想用に振り分け》て、《どれが0,1,2に対応するのか》を《機械学習アルゴリズム》に通して、《予想して》いく流れが見えてきます。

さて、これらのデータをpandasを用いて見やすくしていきましょう(/・ω・)/

5. Pandasを用いてデータを分かりやすくしよう！

今のままでは配列が出力されただけで、何のデータなのかよくわからない状況になっています。

そんな時に表を列ごとに名前を付けて見やすくできたり、平均値や標準偏差などを自動で出してくれるというPandasライブラリを使っていきます。

配列をDataFrameに変換しよう！

今からやっていくことはこんな感じ↓

1. DataFrameの第一引数にデータセット、第二引数にカラムの名前を与える。
2. DataFrameの第一引数に目的変数、第二引数にカラムの名前を与える。
3. 1と2のDataFrameを横に結合したDataFrameを作る

↓なのでプログラムはこんな感じになります↓

Pandasの定義する.ipynb

heacet_data = pd.DataFrame(heacet.data, columns=["がく片の長さ","がく片の幅","花びらの長さ","花びらの幅"])
heacet_target = pd.DataFrame(heacet.target, columns=["花の種類"])
heacet_all = pd.concat([heacet_data,heacet_target], axis=1)

恐らく1行目と2行目のプログラムは大体何をしているか分かると思いますが、3行目に『concatメソッド』を用いて、引数に『1行目のデータと2行目のデータを選択』して、『axis=1』と定義しています。

axis=1ってなに？？

『axis』は恐らくよく見かけているのではないでしょうか。今後も使っていくと思うので、しっかりと確認しておきましょう！

【　axisとは軸を指定する引数　】

図で簡単に表すとこんな感じ↓

axis=0だと縦の行を示して、axis=1だと横の列を示していることになります。

そこで今回は、DataFrameであるheacet_dataのcolumnsに《花の種類》というcolumnを横に追加したかったので、axis=1と指定したわけですね。

実際にDataFrameを見ていこう！

さてさて、PandasのDataFrameを使ったことで見やすくなったはずなので、見にいきましょう～

headメソッドを使って最初の10行を見てみます。

とてもまとまっていて、2次元配列を単に出力した時とは見やすさが大違いですね。笑

またPandasでは平均値などを出してくれるdescribeメソッドもあるので、是非実行してみてください！

6. データセットを分割する

そもそも論なんですが、予測するデータはすでに答えが存在してしまっているので、自分で学習させるデータと予測するデータを創り出さなければなりません。

冒頭で訓練用とテスト用に分ける機能をインポートしましたよね。やっとここで使います。
ですが少し待ってください。ある用語を説明していなかったので、ここで挟みます。

説明変数と目的変数とは？

ゴリゴリペイントで書いた感のある図で説明していきます！

上の図のように今まで《がく片の長さ、がく片の幅、花びらの長さ、花びらの幅》といったカラムを付けていた部分を説明変数、《花びらの種類》とカラムを付けていた部分を目的変数と呼びます。
これから多用していくので、是非覚えておいてください！！

そして単純に説明変数と目的変数を使って学習をさせていきたいのですが、まだ学習用とテスト用のデータを分割していませんでした！

え、そもそも学習用とテスト用のデータってなに？

↑上図のように学習用データを使って学習させてモデルを作成して、テスト用データを使ってモデルが正確に動作しているか確かめるという流れが機械学習にはあります。

この一連の流れをこなす為に《学習用データ》と《テスト用データ》が必要となってきます。

実はこの『学習用とテスト用(詳しくはvalidation用)のデータを分ける』という段階は深堀りすれば、Kaggler達で様々な議論が行われているような分野みたいです…

この記事はそんな話題には触れず、サラッと進めていきます。笑

てか『データなんてPythonでシンプルに分ければ良いじゃないか！』という話なのですが、先ほどの《花の種類》の目的変数を見て頂けたら分かる通り、0⇒1⇒2の順番で目的変数が並んでいます。

ここでもしスライスを使ってデータを半分に分けたりすると、上部分の目的変数が0と1のデータしかモデルは識別できないので、目的変数が2であるデータに対して有効的なモデルを作ることができません。

ペイント間満載の図を使っていくと…

↑上では『0のデータ』と『1のデータ』しか学習していません。

↓結果的にこのような事が起こります。

図にしてしまえば単純な話ですね！

そこで、
『データを学習用とテスト用に分ける且つシャッフルしてくれるような機能があれば良いなぁ～』
なんて思う訳です。

『あります。train_test_split関数があります！！！』

散々言葉の定義について書いてきましたが、やっとプログラムを書いていきます。笑
ここで冒頭にてimportしたtrain_test_split関数を使います。

train_test_split関数を使う

今回は学習用の説明変数と目的変数、テスト用の説明変数と目的変数を以下のように定義します。

• 学習用の説明変数　⇒　setsumei_train
• 学習用の目的変数　⇒　mokuteki_train
• テスト用の説明変数　⇒　setsumei_test
• テスト用の目的変数　⇒　mokuteki_test

これらをプログラム化するとこうです↓

学習用とテスト用の説明変数と目的変数を定義する.ipynb

setsumei_train,setsumei_test,mokuteki_train,mokuteki_test = train_test_split(heacet_data, heacet_target, test_size=0.33)

引数としてtest_sizeというものがあります。これはデータセットをどれだけテスト用に使うかを割合で設定する引数です。

今回は『データセット50個分 = 全体の1/3 = 0.33』をテスト用に使うので、引数を『0.33』としています。

しっかりと分けれているか見ていきます。僕はこうなりました↓

良い感じにシャッフルされて、100個50個で分かれていますね。またID番号も同じなので良い感じです( *´艸｀)

7. データをmatplotlibで可視化する

次は効果的な特徴量を見つけるために、matplotlibでデータを可視化していきます。

データを可視化する意味は？《特徴量選択》

特徴量とは説明変数のことを指していて、『特徴量　≠　説明変数』であることは覚えておいて欲しいのですが、今は同じものとしておきます。

さてデータを可視化することで、効果的な特徴量選択が行えるということを例えを出して説明していきます。

グラフ上にプロットしている赤と緑と青の点についてそれぞれで分類したいとします。

このとき↓下の２つの図ではどちらの方が正しく分類できるでしょうか？

↑これが1つ目の図

↑これが2つ目の図

どう見ても2つ目の図のほうが正しく分類できそうですよね。実際にそうです。笑

少し極端な例でしたが、上下の図のように選択する特徴量によって、分類のしやすさが変化することがあります。

今回のアヤメの分類で言えば、《がく片の長さ》《がく片の幅》《花びらの長さ》《花びらの幅》で2つ特徴量を選ぶとすれば、組み合わせによって分類のしやすさが変わってくるということです。

散布図をplt.scatterで描いていく！

matplotlibのscatterメソッドは引数として以下を指定します。

1. 横軸にしたいデータ
2. 縦軸にしたいデータ
3. label=凡例
4. cmap=カラーマップの種類

カラーマップを用いて色分けしていくのですが、様々な色分けの方法が存在します。
>>>指定できるカラーマップの一覧はコチラ(matplotlib公式リファレンス)

早速matplotlibを使ったプログラムを書いていきたいのですが、花の種類によって視覚的に分類の正確さを測るために目的変数別に色を変えていきます。

目的変数をフィルタリングする方法が僕は少し躓きました。詳しい方がいたらコメントを頂けたら嬉しいですm(__)m

縦軸と横軸データの選択時に上下どちらの引数を使っても大丈夫だったので、実際に左右で使っています。

• setsumei_train[(mokuteki_train == 0).values]["がく片の長さ"]
• setsumei_train[mokuteki_train["花の種類"] == 0]["がく片の長さ"]

がく片の長さと幅を使って、プロットしていきます。

pltでプロットする(がく片).ipynb

plt.scatter(setsumei_train[mokuteki_train["花の種類"]==0]["がく片の長さ"],setsumei_train[(mokuteki_train == 0).values]["がく片の幅"],label="setosa",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==1]["がく片の長さ"],setsumei_train[(mokuteki_train == 1).values]["がく片の幅"],label="versicolor",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==2]["がく片の長さ"],setsumei_train[(mokuteki_train == 2).values]["がく片の幅"],label="virginica",cmap="rgb")

## X軸の範囲を指定
plt.xlim(3,9)
## Y軸の範囲を指定
plt.ylim(1,5)

## X軸の名前
plt.xlabel("Length of sepal")
## Y軸の名前
plt.ylabel("Width of sepal")

## グラフのタイトル
plt.title("Relation between length and width of sepal")
## 凡例を出力
plt.legend()

先ほど上下どちらでも良いと言いましたが、目的変数をフィルタリングする方法が一つは『Numpyで2次元でブールインデックス参照をしている？』のと、『Seriesのブール値をDataFrameに入れることによって、フィルタリングしている？』方法があります。
この考えが正しいのか証明できる文献がイマイチ探せなかったので、詳しい方がいたらコメントを頂けたら嬉しいですm(__)m

一応その結果は残しておきます。

zikken.py

>>> mokuteki_train == 0
## ブール値でDataFrameが返ってくる
>>> mokuteki_train["花の種類"] == 0
## ブール値でSeriesが返ってくる
>>> (mokuteki_train == 0).values
## np行列でブール値の2次元配列が返ってくる
>>> (mokuteki_train["花の種類"] == 0).values
## np行列でブール値の1次元配列が返ってくる

# この時setsumei_train[mokuteki_train["花の種類"]==0]と
# setsumei_train[(mokuteki_train == 0).values]が全く同じ値を返します。

また花びらの長さと幅でプロットする例も示しておきます。

pltでプロットする(花びら).ipynb

plt.scatter(setsumei_train[mokuteki_train["花の種類"]==0]["花びらの長さ"],setsumei_train[(mokuteki_train == 0).values]["花びらの幅"],label="setosa",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==1]["花びらの長さ"],setsumei_train[(mokuteki_train == 1).values]["花びらの幅"],label="versicolor",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==2]["花びらの長さ"],setsumei_train[(mokuteki_train == 2).values]["花びらの幅"],label="virginica",cmap="rgb")

plt.xlim(0,8)
plt.ylim(0,4)

plt.xlabel("Length of petal")
plt.ylabel("Width of petal")

plt.title("Relation between length and width of petal")
plt.legend()

そして↓下図ががく片の長さと幅を使ってプロットした結果

これは分類しにくそうですね。人間である僕からしても、どこに区別する線を引いたら良いのかわかりません。笑

そして↓下図が花びらの長さと幅を使ってプロットした結果

これはそれぞれ目的変数の集合ができているので、簡単に分類する境界線が引けそうですね！
このように選択する特徴量によって、分類のしやすさが変わってくることが視覚的に分かりました。

8. 機械学習アルゴリズムを使っていこう！

ついに機械学習のアルゴリズムに触れることができます。

といってもLinearSVCを用いることを決めてしまっているので、もうそこまでやることはないです。

『ハイパーパラメータを決める』といった段階もありますが、今回は『機械学習の一連の流れをつかむ』ことを重きに置いているので触れないでおきます。気になる方は下のリンクで参考文献を貼っておきます。
>>>LinearSVCのハイパーパラメータの詳しい説明はコチラ

特徴量を選択しよう！

では先ほどMatplotlibを使って可視化させた2つのパターン(がく片コンビと花びらコンビ)の特徴量を選択していきます。
新たに2つのDataFrameとして定義すると、《名前による参照メソッドloc》を使って、プログラムは以下になる。

2つの特徴量を作成する.ipynb

## がく片コンビのDataFrameを作成する。
gakuhen_train = setsumei_train.loc[:,["がく片の長さ","がく片の幅"]]

## 花びらコンビのDataFrameを作成する。
hanabira_train = setsumei_train.loc[:,["花びらの長さ","花びらの幅"]]

LinearSVCでモデルを構築⇒学習⇒予測させる

やっとモデル構築までたどり着きました。笑

もう一度ここ付近の話を説明した図を引っ張ってきますと、

今までは『与えられたデータセットがどうなっているか見たり』、『学習用とテスト用でデータを分けたり』と上図での上の方をずっとやっていました。

ですがこれが『機械学習という分野の一種の特徴』だそうで、機械学習は『前処理に時間がとられる』と聞いたことありませんか？？

まさにそれを具現化しちゃいましたね。記事の大半を占めてしまっています。

そして次にやる機械学習アルゴリズムを実装させる部分は『機械学習アルゴリズムを理解するのは難しい』けれど、『実装させるのは超簡単』という分野で、すぐに終わります。

よって冒頭で定義したLinearSVCを使っていきます。プログラムはこんな感じ↓

LinearSVCでモデル構築⇒学習⇒予測する.ipynb

## それぞれモデルを構築
## それぞれモデルを構築
gakuhen_model = LinearSVC()
hanabira_model = LinearSVC()

## それぞれのモデルに学習させる
gakuhen_model.fit(gakuhen_train,mokuteki_train)
hanabira_model.fit(hanabira_train,mokuteki_train)

## それぞれのモデルで予測させて、予測値を代入させる
### モデルが《がく片の長さと幅》を使って学習しているので、予測する時も《がく片の長さと幅》を渡す必要がある。
gakuhen_predict = gakuhen_model.predict(setsumei_test.loc[:,["がく片の長さ","がく片の幅"]])
### モデルが《花びらの長さと幅》を使って学習しているので、予測する時も《花びらの長さと幅》を渡す必要がある。
hanabira_predict = hanabira_model.predict(setsumei_test.loc[:,["花びらの長さ","花びらの幅"]])

ついに答え合わせです！ここで上で立てた仮説を検証できますね。

Matplotlibで可視化させた図では『《がく片の幅と長さ》は正確に分類できなさそう』でした！
>>>がく片の長さと幅を可視化したグラフはコチラ

逆に『《花びらの幅と長さ》は《がく片の幅と長さ》よりかは正確に分類できそう』でしたね！
>>>花びらの長さと幅を可視化したグラフはコチラ

モデルが予想したデータの答え合わせ

では答え合わせできるプログラムをインポートして、実行していきます。

accuracy_scoreを使う.ipynb

## sklearnライブラリからscore算出の関数をimport
from sklearn.metrics import accuracy_score

## gakuhen_scoreとhanabira_scoreにそれぞれに結果を代入
gakuhen_score = accuracy_score(mokuteki_test, gakuhen_predict)
hanabira_score = accuracy_score(mokuteki_test, hanabira_predict)

print('がく片の長さと幅コンビの正解率：{}'.format(gakuhen_score),'花びらの長さと幅コンビの正解率：{}'.format(hanabira_score), sep='\n')

↓僕は出力結果として以下が得られました！↓

とりあえずどちらも正答率が8割を超えているので、機械学習によってアヤメの種類を分類することはひとまず成功しましたね！

ではでは…

結果について考察していきます。

2つの結果の違いについて詳しく見る

人間の目で見ても《がく片の長さと幅》より《花びらの長さと幅》の方が正確に分類しやすそうでしたが、LinearSVCアルゴリズムにとっても同じように《花びらの長さと幅》の方が正確に分類できるようです。笑

ですが『LinearSVCアルゴリズムにとっても同じように《花びらの長さと幅》の方が正確に分類できる』というのは僕の推論でしかないのです。

実際に確かめるには、どこで境界を作っているのかをMatplotlibを使って可視化していくと良いですよね！！

《がく片の長さと幅》と《花びらの長さと幅》でそれぞれ境界線を見たいので、代入できる関数として定義していきます。

分類の境界を可視化する.ipynb

def heacet_border_check(H, M, model, param1, param2, resolution=0.01):
    H1_min, H1_max = H[param1].min()-0.5, H[param1].max()+0.5
    H2_min, H2_max = H[param2].min()-0.5, H[param2].max()+0.5
    H1, H2 = np.meshgrid(np.arange(H1_min, H1_max, resolution),
                           np.arange(H2_min, H2_max, resolution))
    n = np.array([H1.ravel(), H2.ravel()]).T
    Z = model.predict(n)
    Z = Z.reshape(H1.shape)
    plt.contourf(H1, H2, Z, alpha=0.5, cmap="Set2")
    plt.xlim(H1_min, H1_max)
    plt.ylim(H2_min, H2_max)
    plt.xlabel("Length")
    plt.ylabel("Width")
    plt.scatter(H[param1],H[param2], c=M["花の種類"], cmap="brg")

関数が用意できたので、コチラを使って出力すると…

↑まずは上図が《がく片の長さと幅》を使った時の境界図です。
青丸は上手く境界を持てているようですが、赤丸と緑丸ごちゃごちゃしていて、微妙なところに境界線が引かれていますね。

そりゃあスコア低くなりますよね…という感じ。

↑続いて上図が《花びらの長さと幅》を使った時の境界図です。
青丸は完璧ですね。赤丸と緑丸同士がほんの少し境界を越えているくらいで、ほぼ綺麗に境界線を引けていると思います。

がく片と花びらとで場合分けしてきましたが、『データの可視化によって特徴を選択することはとても重要なこと』は伝わったのではないでしょうか？

まとめ

最後まで見て頂きありがとうございました。

普段から技術記事は書いていないもので、少し口語が多かったかもしれません…笑

色々と調べながら初心者だからこその目線で《アヤメ分類》について一から説明してみました。

文書が変だったり、間違っている点などございましたら気軽にコメント頂けると嬉しいです。
もちろん感想やよかった点などでも気軽にコメントください(´艸｀*)

こんにちは。ひろちょんです。

機械学習を始めたばかりの初心者が、様々な記事を参考にしつつアヤメ分類をしてみます。
説明を加えながら書いていくので、是非参考にしてみてください！

↓記事を読む前に筆者について知りたい方へ↓
>>>詳しいプロフィールはコチラ

～目次です～
1. 対象読者とか書いてみる
2. 《Google colaboratory》で進めていきます
3. 色々とimportしていく
　　・オートパイロット状態でimportする奴
　　・scikit-learnの様々な機能をimport!!
4. アヤメのデータセットをじっくり見ていく
　　・インスタンスを生成
　　・アヤメデータはどうなってるの？
5. Pandasを用いてデータを分かりやすくする
　　・配列をDataFrameに変換しよう！
　　・実際にDataFrameを見ていこう！
6. データセットを分割する
　　・説明変数と目的変数とは？
　　・学習用とテスト用のデータはどうするの？
　　・train_test_split関数を使う！
7. データをmatplotlibで可視化する
　　・データを可視化する意味は？《特徴量選択》
　　・散布図をplt.scatterで描いていく！
8. 機械学習アルゴリズムを使っていこう！
　　・特徴量を選択しよう！
　　・モデルを構築⇒学習⇒予測させる
　　・モデルが予想したデータの答え合わせ
　　・2つの結果の違いについて詳しく見る

1. 対象読者とか書いてみる

今回は機械学習でアヤメの分類を行っていくことが主題なので、レベルとしてはこんな感じ。

• 何かしらのNotebook形式でプログラムを実行できる人
• Python初心者以上
• 機械学習の流れが何となくわかる人

2. Google Colaboratoryで進めていきます

今回はバージョンとか、仮想環境とか…諸々のしがらみに囚われたくなかったので、Google Colaboratoryにてアヤメ分類をしていきたいと思います。

始め方はすごく簡単で、Google Colaboratoryのサイト左上の《ファイル》から《Python3の新しいノートブック》をクリックしてください。

すると新しいノートブックが生成されるので、左上のUntitled0.ipynbから名前を変更してあげてください！

※僕は名前を『初めてのアヤメ分類』にしておきました。笑

初学者の方はページを見つつ、一から作成することをオススメしますが、一応完成形も↓のGitHubにてノートブックを公開しています。
>>>GitHubへはコチラから

3. 色々importしていく

オートパイロット状態でimportする奴

とりあえず以下をインポートしていきますね。(矢印先は僕のライブラリへの感想です。)

• Numpy　⇒　計算機能凄めライブラリ
• Pandas　⇒　表を扱いやすいライブラリ
• Matplotlib　⇒　グラフ可視化ライブラリ
• warnings　⇒　たまに邪魔な警告を消すライブラリ

そしてプログラム化するとこんな感じ。

importするライブラリ.ipynb

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')

scikit-learnの色々な機能をimport!!

機械学習ライブラリのscikit-learnには分類を行う以外にも様々な機能が存在しています。今回使用するものをサラっと紹介しますね。(後で深堀りしていきます)

1つ目はscikit-learnに用意されているデータセットです。
scikit-learnのデータセット集一つにアヤメのデータセットがあります。
>>>詳しくはコチラ(UCI Machine Learning Repository: Iris Data Set)

２つ目が学習用とテスト用にデータセットを分ける機能です。
分ける理由としては、《学習用のデータで予測するモデルを作成する》⇒《テスト用データで作成したモデルは良い物か判定する》という機械学習では鉄板の2つの処理を行いたいからです。

3つ目が分類を行うアルゴリズムになります。機械学習には様々な分類のアルゴリズムが存在して、以下のプログラムではアルゴリズムとして『線形のSVM』を用いるので、LinearSVCをインポートしています。
>>>SVM(サポートベクターマシーン)について詳しく知りたい方はコチラ(Qiita)

では上記で紹介した3つをimportしていきます。

scikit-learnの様々な機能をimport.ipynb

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.svm import LinearSVC

4. アヤメのデータセットをじっくり見ていく

まず今回の目的をハッキリさせておきましょう！

どのサイトでも分類の部分を重要視しすぎて、結局何がしたかったのか見失いがちです。

今回の目的は…

『アヤメのがく片や花びら幅や長さの数値』を用いて『アヤメ属(花)の種類を分類する』こと

そこで既にアヤメの種類によっての幅や長さのデータを数値を集めてくれているデータセット(先ほどimportした奴)を使って、機械学習により予測していきます。

フィッシャーのアヤメのデータから散布図行列を描く(https://teenaka.at.webry.info/201803/article_12.html )から引用

ちなみに上の図がアヤメになります。この花の長さや幅を測定して数値化したデータセットをimportしました(/・ω・)/

インスタンスを生成

今はアヤメのデータセット機能をインポートしただけなので、まずインスタンスを生成していきます。

アヤメインスタンスを生成する.ipynb

heacet = load_iris()

インスタンスの名前に特に深い意味はないです。単に自サイトの宣伝をしているだけです！笑

もし嫌だという方はコピーするだけではなく、自分で書いていくと理解も深くなるので、ここのインスタンス名を変えて実行してみてください！

アヤメデータはどうなってるの？

さっそくデータセットの中身を見ていきましょう。
↓コチラを実行してください↓

アヤメデータの中身を確認する.ipynb

print("与えられたデータ")
print(heacet.data)
print(heacet.data.shape)
print("-----------------")
print("予測するデータ")
print(heacet.target)
print(heacet.target.shape)
print(heacet.target_names)

実行して出力を見ると、どうやら《数値が書かれた150×4の2次元配列》と《0,1,2と書かれた150×1の1次元配列》が得られたみたいですね。

また0,1,2はそれぞれ『setosa』,『versicolor』,『virginica』対応していることも分かります。

よって与えられた《数値が書かれた150×4の2次元配列》を《データセットを学習用と予想用に振り分け》て、《どれが0,1,2に対応するのか》を《機械学習アルゴリズム》に通して、《予想して》いく流れが見えてきます。

さて、これらのデータをpandasを用いて見やすくしていきましょう(/・ω・)/

5. Pandasを用いてデータを分かりやすくしよう！

今のままでは配列が出力されただけで、何のデータなのかよくわからない状況になっています。

そんな時に表を列ごとに名前を付けて見やすくできたり、平均値や標準偏差などを自動で出してくれるというPandasライブラリを使っていきます。

配列をDataFrameに変換しよう！

今からやっていくことはこんな感じ↓

1. DataFrameの第一引数にデータセット、第二引数にカラムの名前を与える。
2. DataFrameの第一引数に目的変数、第二引数にカラムの名前を与える。
3. 1と2のDataFrameを横に結合したDataFrameを作る

↓なのでプログラムはこんな感じになります↓

Pandasの定義する.ipynb

heacet_data = pd.DataFrame(heacet.data, columns=["がく片の長さ","がく片の幅","花びらの長さ","花びらの幅"])
heacet_target = pd.DataFrame(heacet.target, columns=["花の種類"])
heacet_all = pd.concat([heacet_data,heacet_target], axis=1)

恐らく1行目と2行目のプログラムは大体何をしているか分かると思いますが、3行目に『concatメソッド』を用いて、引数に『1行目のデータと2行目のデータを選択』して、『axis=1』と定義しています。

axis=1ってなに？？

『axis』は恐らくよく見かけているのではないでしょうか。今後も使っていくと思うので、しっかりと確認しておきましょう！

【　axisとは軸を指定する引数　】

図で簡単に表すとこんな感じ↓

axis=0だと縦の行を示して、axis=1だと横の列を示していることになります。

そこで今回は、DataFrameであるheacet_dataのcolumnsに《花の種類》というcolumnを横に追加したかったので、axis=1と指定したわけですね。

実際にDataFrameを見ていこう！

さてさて、PandasのDataFrameを使ったことで見やすくなったはずなので、見にいきましょう～

headメソッドを使って最初の10行を見てみます。

とてもまとまっていて、2次元配列を単に出力した時とは見やすさが大違いですね。笑

またPandasでは平均値などを出してくれるdescribeメソッドもあるので、是非実行してみてください！

6. データセットを分割する

そもそも論なんですが、予測するデータはすでに答えが存在してしまっているので、自分で学習させるデータと予測するデータを創り出さなければなりません。

冒頭で訓練用とテスト用に分ける機能をインポートしましたよね。やっとここで使います。
ですが少し待ってください。ある用語を説明していなかったので、ここで挟みます。

説明変数と目的変数とは？

ゴリゴリペイントで書いた感のある図で説明していきます！

上の図のように今まで《がく片の長さ、がく片の幅、花びらの長さ、花びらの幅》といったカラムを付けていた部分を説明変数、《花びらの種類》とカラムを付けていた部分を目的変数と呼びます。
これから多用していくので、是非覚えておいてください！！

そして単純に説明変数と目的変数を使って学習をさせていきたいのですが、まだ学習用とテスト用のデータを分割していませんでした！

え、そもそも学習用とテスト用のデータってなに？

↑上図のように学習用データを使って学習させてモデルを作成して、テスト用データを使ってモデルが正確に動作しているか確かめるという流れが機械学習にはあります。

この一連の流れをこなす為に《学習用データ》と《テスト用データ》が必要となってきます。

実はこの『学習用とテスト用(詳しくはvalidation用)のデータを分ける』という段階は深堀りすれば、Kaggler達で様々な議論が行われているような分野みたいです…

この記事はそんな話題には触れず、サラッと進めていきます。笑

てか『データなんてPythonでシンプルに分ければ良いじゃないか！』という話なのですが、先ほどの《花の種類》の目的変数を見て頂けたら分かる通り、0⇒1⇒2の順番で目的変数が並んでいます。

ここでもしスライスを使ってデータを半分に分けたりすると、上部分の目的変数が0と1のデータしかモデルは識別できないので、目的変数が2であるデータに対して有効的なモデルを作ることができません。

ペイント間満載の図を使っていくと…

↑上では『0のデータ』と『1のデータ』しか学習していません。

↓結果的にこのような事が起こります。

図にしてしまえば単純な話ですね！

そこで、
『データを学習用とテスト用に分ける且つシャッフルしてくれるような機能があれば良いなぁ～』
なんて思う訳です。

『あります。train_test_split関数があります！！！』

散々言葉の定義について書いてきましたが、やっとプログラムを書いていきます。笑
ここで冒頭にてimportしたtrain_test_split関数を使います。

train_test_split関数を使う

今回は学習用の説明変数と目的変数、テスト用の説明変数と目的変数を以下のように定義します。

• 学習用の説明変数　⇒　setsumei_train
• 学習用の目的変数　⇒　mokuteki_train
• テスト用の説明変数　⇒　setsumei_test
• テスト用の目的変数　⇒　mokuteki_test

これらをプログラム化するとこうです↓

学習用とテスト用の説明変数と目的変数を定義する.ipynb

setsumei_train,setsumei_test,mokuteki_train,mokuteki_test = train_test_split(heacet_data, heacet_target, test_size=0.33)

引数としてtest_sizeというものがあります。これはデータセットをどれだけテスト用に使うかを割合で設定する引数です。

今回は『データセット50個分 = 全体の1/3 = 0.33』をテスト用に使うので、引数を『0.33』としています。

しっかりと分けれているか見ていきます。僕はこうなりました↓

良い感じにシャッフルされて、100個50個で分かれていますね。またID番号も同じなので良い感じです( *´艸｀)

7. データをmatplotlibで可視化する

次は効果的な特徴量を見つけるために、matplotlibでデータを可視化していきます。

データを可視化する意味は？《特徴量選択》

特徴量とは説明変数のことを指していて、『特徴量　≠　説明変数』であることは覚えておいて欲しいのですが、今は同じものとしておきます。

さてデータを可視化することで、効果的な特徴量選択が行えるということを例えを出して説明していきます。

グラフ上にプロットしている赤と緑と青の点についてそれぞれで分類したいとします。

このとき↓下の２つの図ではどちらの方が正しく分類できるでしょうか？

↑これが1つ目の図

↑これが2つ目の図

どう見ても2つ目の図のほうが正しく分類できそうですよね。実際にそうです。笑

少し極端な例でしたが、上下の図のように選択する特徴量によって、分類のしやすさが変化することがあります。

今回のアヤメの分類で言えば、《がく片の長さ》《がく片の幅》《花びらの長さ》《花びらの幅》で2つ特徴量を選ぶとすれば、組み合わせによって分類のしやすさが変わってくるということです。

散布図をplt.scatterで描いていく！

matplotlibのscatterメソッドは引数として以下を指定します。

1. 横軸にしたいデータ
2. 縦軸にしたいデータ
3. label=凡例
4. cmap=カラーマップの種類

カラーマップを用いて色分けしていくのですが、様々な色分けの方法が存在します。
>>>指定できるカラーマップの一覧はコチラ(matplotlib公式リファレンス)

早速matplotlibを使ったプログラムを書いていきたいのですが、花の種類によって視覚的に分類の正確さを測るために目的変数別に色を変えていきます。

目的変数をフィルタリングする方法が僕は少し躓きました。詳しい方がいたらコメントを頂けたら嬉しいですm(__)m

縦軸と横軸データの選択時に上下どちらの引数を使っても大丈夫だったので、実際に左右で使っています。

• setsumei_train[(mokuteki_train == 0).values]["がく片の長さ"]
• setsumei_train[mokuteki_train["花の種類"] == 0]["がく片の長さ"]

がく片の長さと幅を使って、プロットしていきます。

pltでプロットする(がく片).ipynb

plt.scatter(setsumei_train[mokuteki_train["花の種類"]==0]["がく片の長さ"],setsumei_train[(mokuteki_train == 0).values]["がく片の幅"],label="setosa",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==1]["がく片の長さ"],setsumei_train[(mokuteki_train == 1).values]["がく片の幅"],label="versicolor",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==2]["がく片の長さ"],setsumei_train[(mokuteki_train == 2).values]["がく片の幅"],label="virginica",cmap="rgb")

## X軸の範囲を指定
plt.xlim(3,9)
## Y軸の範囲を指定
plt.ylim(1,5)

## X軸の名前
plt.xlabel("Length of sepal")
## Y軸の名前
plt.ylabel("Width of sepal")

## グラフのタイトル
plt.title("Relation between length and width of sepal")
## 凡例を出力
plt.legend()

先ほど上下どちらでも良いと言いましたが、目的変数をフィルタリングする方法が一つは『Numpyで2次元でブールインデックス参照をしている？』のと、『Seriesのブール値をDataFrameに入れることによって、フィルタリングしている？』方法があります。
この考えが正しいのか証明できる文献がイマイチ探せなかったので、詳しい方がいたらコメントを頂けたら嬉しいですm(__)m

一応その結果は残しておきます。

zikken.py

>>> mokuteki_train == 0
## ブール値でDataFrameが返ってくる
>>> mokuteki_train["花の種類"] == 0
## ブール値でSeriesが返ってくる
>>> (mokuteki_train == 0).values
## np行列でブール値の2次元配列が返ってくる
>>> (mokuteki_train["花の種類"] == 0).values
## np行列でブール値の1次元配列が返ってくる

# この時setsumei_train[mokuteki_train["花の種類"]==0]と
# setsumei_train[(mokuteki_train == 0).values]が全く同じ値を返します。

また花びらの長さと幅でプロットする例も示しておきます。

pltでプロットする(花びら).ipynb

plt.scatter(setsumei_train[mokuteki_train["花の種類"]==0]["花びらの長さ"],setsumei_train[(mokuteki_train == 0).values]["花びらの幅"],label="setosa",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==1]["花びらの長さ"],setsumei_train[(mokuteki_train == 1).values]["花びらの幅"],label="versicolor",cmap="rgb")
plt.scatter(setsumei_train[mokuteki_train["花の種類"]==2]["花びらの長さ"],setsumei_train[(mokuteki_train == 2).values]["花びらの幅"],label="virginica",cmap="rgb")

plt.xlim(0,8)
plt.ylim(0,4)

plt.xlabel("Length of petal")
plt.ylabel("Width of petal")

plt.title("Relation between length and width of petal")
plt.legend()

そして↓下図ががく片の長さと幅を使ってプロットした結果

これは分類しにくそうですね。人間である僕からしても、どこに区別する線を引いたら良いのかわかりません。笑

そして↓下図が花びらの長さと幅を使ってプロットした結果

これはそれぞれ目的変数の集合ができているので、簡単に分類する境界線が引けそうですね！
このように選択する特徴量によって、分類のしやすさが変わってくることが視覚的に分かりました。

8. 機械学習アルゴリズムを使っていこう！

ついに機械学習のアルゴリズムに触れることができます。

といってもLinearSVCを用いることを決めてしまっているので、もうそこまでやることはないです。

『ハイパーパラメータを決める』といった段階もありますが、今回は『機械学習の一連の流れをつかむ』ことを重きに置いているので触れないでおきます。気になる方は下のリンクで参考文献を貼っておきます。
>>>LinearSVCのハイパーパラメータの詳しい説明はコチラ

特徴量を選択しよう！

では先ほどMatplotlibを使って可視化させた2つのパターン(がく片コンビと花びらコンビ)の特徴量を選択していきます。
新たに2つのDataFrameとして定義すると、《名前による参照メソッドloc》を使って、プログラムは以下になる。

2つの特徴量を作成する.ipynb

## がく片コンビのDataFrameを作成する。
gakuhen_train = setsumei_train.loc[:,["がく片の長さ","がく片の幅"]]

## 花びらコンビのDataFrameを作成する。
hanabira_train = setsumei_train.loc[:,["花びらの長さ","花びらの幅"]]

LinearSVCでモデルを構築⇒学習⇒予測させる

やっとモデル構築までたどり着きました。笑

もう一度ここ付近の話を説明した図を引っ張ってきますと、

今までは『与えられたデータセットがどうなっているか見たり』、『学習用とテスト用でデータを分けたり』と上図での上の方をずっとやっていました。

ですがこれが『機械学習という分野の一種の特徴』だそうで、機械学習は『前処理に時間がとられる』と聞いたことありませんか？？

まさにそれを具現化しちゃいましたね。記事の大半を占めてしまっています。

そして次にやる機械学習アルゴリズムを実装させる部分は『機械学習アルゴリズムを理解するのは難しい』けれど、『実装させるのは超簡単』という分野で、すぐに終わります。

よって冒頭で定義したLinearSVCを使っていきます。プログラムはこんな感じ↓

LinearSVCでモデル構築⇒学習⇒予測する.ipynb

## それぞれモデルを構築
## それぞれモデルを構築
gakuhen_model = LinearSVC()
hanabira_model = LinearSVC()

## それぞれのモデルに学習させる
gakuhen_model.fit(gakuhen_train,mokuteki_train)
hanabira_model.fit(hanabira_train,mokuteki_train)

## それぞれのモデルで予測させて、予測値を代入させる
### モデルが《がく片の長さと幅》を使って学習しているので、予測する時も《がく片の長さと幅》を渡す必要がある。
gakuhen_predict = gakuhen_model.predict(setsumei_test.loc[:,["がく片の長さ","がく片の幅"]])
### モデルが《花びらの長さと幅》を使って学習しているので、予測する時も《花びらの長さと幅》を渡す必要がある。
hanabira_predict = hanabira_model.predict(setsumei_test.loc[:,["花びらの長さ","花びらの幅"]])

ついに答え合わせです！ここで上で立てた仮説を検証できますね。

Matplotlibで可視化させた図では『《がく片の幅と長さ》は正確に分類できなさそう』でした！
>>>がく片の長さと幅を可視化したグラフはコチラ

逆に『《花びらの幅と長さ》は《がく片の幅と長さ》よりかは正確に分類できそう』でしたね！
>>>花びらの長さと幅を可視化したグラフはコチラ

モデルが予想したデータの答え合わせ

では答え合わせできるプログラムをインポートして、実行していきます。

accuracy_scoreを使う.ipynb

## sklearnライブラリからscore算出の関数をimport
from sklearn.metrics import accuracy_score

## gakuhen_scoreとhanabira_scoreにそれぞれに結果を代入
gakuhen_score = accuracy_score(mokuteki_test, gakuhen_predict)
hanabira_score = accuracy_score(mokuteki_test, hanabira_predict)

print('がく片の長さと幅コンビの正解率：{}'.format(gakuhen_score),'花びらの長さと幅コンビの正解率：{}'.format(hanabira_score), sep='\n')

↓僕は出力結果として以下が得られました！↓

とりあえずどちらも正答率が8割を超えているので、機械学習によってアヤメの種類を分類することはひとまず成功しましたね！

ではでは…

結果について考察していきます。

2つの結果の違いについて詳しく見る

人間の目で見ても《がく片の長さと幅》より《花びらの長さと幅》の方が正確に分類しやすそうでしたが、LinearSVCアルゴリズムにとっても同じように《花びらの長さと幅》の方が正確に分類できるようです。笑

ですが『LinearSVCアルゴリズムにとっても同じように《花びらの長さと幅》の方が正確に分類できる』というのは僕の推論でしかないのです。

実際に確かめるには、どこで境界を作っているのかをMatplotlibを使って可視化していくと良いですよね！！

《がく片の長さと幅》と《花びらの長さと幅》でそれぞれ境界線を見たいので、代入できる関数として定義していきます。

分類の境界を可視化する.ipynb

def heacet_border_check(H, M, model, param1, param2, resolution=0.01):
    H1_min, H1_max = H[param1].min()-0.5, H[param1].max()+0.5
    H2_min, H2_max = H[param2].min()-0.5, H[param2].max()+0.5
    H1, H2 = np.meshgrid(np.arange(H1_min, H1_max, resolution),
                           np.arange(H2_min, H2_max, resolution))
    n = np.array([H1.ravel(), H2.ravel()]).T
    Z = model.predict(n)
    Z = Z.reshape(H1.shape)
    plt.contourf(H1, H2, Z, alpha=0.5, cmap="Set2")
    plt.xlim(H1_min, H1_max)
    plt.ylim(H2_min, H2_max)
    plt.xlabel("Length")
    plt.ylabel("Width")
    plt.scatter(H[param1],H[param2], c=M["花の種類"], cmap="brg")

関数が用意できたので、コチラを使って出力すると…

↑まずは上図が《がく片の長さと幅》を使った時の境界図です。
青丸は上手く境界を持てているようですが、赤丸と緑丸ごちゃごちゃしていて、微妙なところに境界線が引かれていますね。

そりゃあスコア低くなりますよね…という感じ。

↑続いて上図が《花びらの長さと幅》を使った時の境界図です。
青丸は完璧ですね。赤丸と緑丸同士がほんの少し境界を越えているくらいで、ほぼ綺麗に境界線を引けていると思います。

がく片と花びらとで場合分けしてきましたが、『データの可視化によって特徴を選択することはとても重要なこと』は伝わったのではないでしょうか？

まとめ

最後まで見て頂きありがとうございました。

普段から技術記事は書いていないもので、少し口語が多かったかもしれません…笑

色々と調べながら初心者だからこその目線で《アヤメ分類》について一から説明してみました。

文書が変だったり、間違っている点などございましたら気軽にコメント頂けると嬉しいです。
もちろん感想やよかった点などでも気軽にコメントください(´艸｀*)

記事へのリンク

1. 理論編（←本記事）
2. 実装編（書いてます・・・）
3. 評価編（書いてます・・・）

本記事におけるお断り

本記事は筆者が自学自習のためにデータ取得から計算までを行ったものです。特定の団体や個人から許可許諾はとっていないため、もし内容に問題があるとお考えの場合にはご連絡いただければ幸いです。
また、タイトルにもある通り、単純に強い、上手いを判断するのが困難であるポーカーにおけるスキルの『推定』がやりたかったことです。
本記事を通じて少しでもポーカーに興味を持ってアミューズメントやオンラインでポーカーを始めたいと思う人が増えれば何よりです。

なんでこんなことしようと思ったのか

pandasを勉強する機会があり、継続学習のために興味がある領域でお勉強＆アウトプットしたかったらです。
あと、ポーカーというかなり運の要素が入るゲームにおいて実力を正確に計算する事が出来るのかという点は昔から興味がありました。
昔社交ダンスの試合結果100万件をELOレーティングで計算したりと、何かと順位を元にした実力推定は興味があったのです。

背景の説明

ポーカー

現在世界ではテキサスホールデムという種目が最も主流でありプレイ人口も多いです。
本記事では今後ポーカー＝ノーリミットテキサスホールデムという前提で記述します。
5枚引いて、好きな枚数チェンジして、勝負！ではなく、裏向きの自分しか使えない2枚＋全員が使える5枚のうち好きな5枚を使って勝負！
という感じのルールです。詳細は以下URLよりぜひご覧ください！
https://ja.wikipedia.org/wiki/%E3%83%86%E3%82%AD%E3%82%B5%E3%82%B9%E3%83%BB%E3%83%9B%E3%83%BC%E3%83%AB%E3%83%87%E3%83%A0

日本のポーカー事情

ご存知のとおり日本において賭博は違法です。そんな中どのように日本のプレイヤーはポーカーを楽しんでいるかというと、大きく3つに分けられます。

1. アミューズメント施設で遊ぶ
2. オンラインゲームで遊ぶ
3. 海外のカジノに行って遊ぶ

それぞれざっくり言うと以下のような特徴があります。

名称	特徴	代表的な場所	データの有無
1.アミューズメント施設	ゲームセンターと同じような仮想チップを使ってゲームを楽しむ。換金は出来ないが大きな大会ではスポンサーから商品が出る事もある。	アキバギルド、東京dePoker、カジノクエスト、パラハetc	全順位データ公開
2.オンラインゲーム	モバイルやPCから遊べる。一部サイトではお金を賭けて遊べてしまったりするが、合法か違法かはグレーなので自己責任。	PokerStars、Zynga、PokerPoker	一部開発者のみに限定公開
3. 海外のカジノ	カジノが有ってポーカーが遊べる国に行けば合法的にお金を賭けてポーカーが出来る。全てのカジノでポーカーが出来るわけではない（むしろ全体でみると少ない？）ので注意。	ラスベガス、マカオ、韓国、マニラ、バルセロナ、モナコetc	入賞者のみ公開

今回、「1.アミューズメント施設で遊ぶ」という試合結果をもとにして計算を行いました。
なぜかと言うと、データが公開された状態で十分な量存在していることと、結果に対する肌感があるため自分でも計算結果を一定評価出来るからです。

計算手法

今回は「総当たり式ELOレーティング」と「TrueSkill」の2つの方法で計算を行いました

総当たり式ELOレーティング

ELOレーティングはチェスの実力評価のために開発された（恐らく）最古の実力を定量化するためのアルゴリズムです。
https://ja.wikipedia.org/wiki/%E3%82%A4%E3%83%AD%E3%83%AC%E3%83%BC%E3%83%86%E3%82%A3%E3%83%B3%E3%82%B0

ただし、1vs1の試合結果しか計算出来ないため、将棋・チェス・ボクシング・剣道etcの競技では単純に適用可能なのですが、
ポーカー、ゴルフ、競馬、社交ダンスetcの様な複数のプレイヤーが同時に競技を行い順位がつく競技に適用するためには工夫が必要です。

その工夫が「総当たり式」であり、端的に言えば1試合の参加者を全2人の組み合わせに分解してそれぞれで勝敗を計算してしまおう、という考えです。
例：ABCDEの5人が100m走を走り、CABEDの順位だったとする
- AvsB：Aの勝ち、Bの負け
- AvsC：Cの勝ち、Aの負け
- （略）
- DvsE：Eの勝ち、Dの負け

といった感じでそれぞれの結果を計算し最終的に合算します
計算式については上記Wikipediaリンクに詳細がありますが、抜粋するとこんな感じです。

勝率の計算

Aさん（レーティング1700）とBさん（レーティング1500）が対戦する場合を考えます
まずAさんとBさんの対戦した時の勝率を計算します
Aさんから見た勝率は以下の式より約76%と計算できます

W_{ab} = \frac{1}{10^{(R_B-R_A)/400}+1} 

レーティング変動を計算

上記で算出した勝率をもとに、試合後のレーティング（R'）がどの程度変動するか計算します。

R'_A = R_A + K \times W_{BA}

R'_B = R_B - K \times W_{BA}

定数パラメーター

ELOレーティングの計算において、3つのパラメータを決定する必要があります。

• 初期レーティング：多くの場合1500と設定されていますが好きな数字でOK
• 初期レーティングに対する勝敗比：多くの場合400と置かれています。初期レーティングに対して乖離しすぎていないければOKの認識
• 変動率（K）：1回の試合結果でどの程度レーティングを変動させるのかのパラメータです。32が一般的と言われていますが、競技の特性によってここの値を適正に変える必要があります。試合頻度が高い、1対1ではなく順位式の競技、運の要素が大きい競技、などの場合にはKの値を小さく設定するのが良いです

計算してみた結果から言うと・・・

とてもではないが使いものになりませんでした。。。
理由はいくつか考えられますが、ポーカーのようなどんな実力者でもかなりの頻度で初心者に負けうる様な競技の場合
レーティングが1回の勝敗に大きく左右されすぎるため、人数が多い試合においては1回の負け/勝ちであまりに大きな変動となってしまう事がわかりました。
例として100人程度の規模である程度の実力者が下位20%くらいの順位だと一発で-2000くらいの変動を食らっていました。
では大きく変動しすぎないようにKを下げれば良いかと言われると、今度は少人数の試合で全く動かなくなるためこれもよろしくない。
では人数にKを反比例させる、等考えられますがそうすると1試合の価値が変わる事になりこれもどうかと。。。

ちなみにこちらは私の約80回あまりの試合ごとの変動をプロットしたものです。
見ての通り、概ね1500より高い位置に居ますが、最後の方で大きくマイナスを食らった結果、
残念ながら1500未満の結果となりました。（悲しい）

TrueSkill

TrueSkillはマイクロソフト謹製のスキル推定アルゴリズムです。
XBoxのゲーム内マッチングで利用する目的で開発され、マイクロソフトが特許を持っており商用利用には許諾が必要です
日本語記事だとQiita内では以下2つが大変参考になりました。
数学が苦にならない人は1つめの記事を苦手な方は2番目の記事を読むと大変わかりやすいと思います

誰でも分かるTrueSkill
TrueSkill「まだ Elo レーティングで消耗してるの?」

以下2つめの記事より引用した特徴です

• 収束が早い。 レーティングに初めて参加するプレイヤーの実力を推定するのに何度も何度も不適当なマッチングで対戦する必要がない。
• 複数人による対戦に対応している。 勝ちか負けかのみならず順位を定めるようなゲームやチーム戦1のゲームにも使用できる。
• ゲームへの参加に重み付けができる。 例えばチームメンバーのひとりが回線トラブルにより途中でゲームから抜けた場合でも適用できる。
• Microsoft の息がかかっている。 特許申請と商標登録がなされており、なんとなく手が出しづらい。 あと一部の狂信的 OSS 主義者が憤死する。

TrueSkillの概要

※ここから先の説明は数学素人の筆者が理解しやすさ重視で記載しているので誤っている箇所があるかもしれません
ELOレーティングは1回スキルの計算が終わると、「今回は良い結果だったからあなたのレーティングは1500⇒1600にアップしました」という様に具体的なスキルの数値を示します。
これを複数回繰り返すことで本来のスキルに収束してゆく、という考え方です。
一方TrueSkillは1回スキルの計算が終わると、「今回は良い結果だったからあなたのレーティングは、試合前は1200～1800の間のどこかで中央の1500である可能性が最も高かったけど、計算後は1400～1700の間のどこかで中央値の1550である可能性が最も高くなるとわかりました」という様にスキルを確率分布で示します。

数学に詳しい方はお気づきかと思いますが、これは正規分布で表現されているという事です。
TrueSkillで計算が進むとは、正規分布の中央の値（μ）が左右に動きその人のスキルに近づいてゆき、更に正規分布の範囲が狭まり（σ）その人のスキルの幅が確定してゆく。
といった過程を経て最終的なスキルが算出されます

こちらもELOレーティングと同様に私自身の計算過程である約80試合の計算過程をプロットしたものです。
低くて横に広いグラフが最初の状態、高くて横に狭いグラフが最終的な状態です。
見れば分かる通り私はグラフの中央が右方向に動いているため、一応計算上は強い人、ということになります。
（一応計算上は、です。。。）

次の記事

次の記事では、どのように上記理論を実装したのか簡単に記載したいと思います。
そして最後に計算結果から分析結果をまとめて記載したいと思います。
ちなみに実名（ポーカーの場合はポーカーネームというニックネーム）を公表すると色々問題になりそうなのですが、
どうしようかは現在まだ決めていません。。。

もし現時点で自分の値を知りたい方はTwitterで本人までDMいただければこっそり教えます！

ターミナルで実行する際にエラーが出た時、ブラウザを起動して検索するプロセスでとても億劫になっていたのでターミナルから検索したい文字を入力すると自動でブラウザが開いて検索してくれるsomethingを作りました。
これでエラーが出たときに気負う事なく検索ができるようになります。
小ネタ程度のものなので、コピペして必ず動くとは言えないほどガバなコードです。
フローだけ参考にして、よかったら自作してみてください。

詳細

コマンドライン引数に検索したいワードを入力し、実行することで入力したワードを検索します。
検索結果は、ブラウザが起動して表示されます。
使用したブラウザは、safariです。
数行で書けます。

コード

search.py

import webbrowser, sys

url = 'https://www.google.com/search?client=safari&rls=en&q={}&ie=UTF-8&oe=UTF-8'
search = ' '.join(sys.argv[1:])
webbrowser.open_new_tab(url.format(search))

実行する際の、python search.py hoge の python search.pyを環境変数にします。

terminal

$ echo export google=python\ search.pyのフルパス >> ~/.bash_profile
$ source .bash_profile

これで簡単にターミナルから検索ができるようになりました。

terminal

$ $google hoge

感想

大したものではないですがこの程度のもので効率が上がるならええやん。と思いました。

正確な発動条件が定かではないですが、ubuntu18.04にpipを手動インストールし、pyenvで環境構築した際にタイトルのエラーが発生したので、解決できた方法をメモします。

現象

上記条件で、pyenvで作成したpython環境でpipを実行しようとすると次のエラーが出て動かない。

>> pip
Traceback (most recent call last):
  File "/home/dev-user/.pyenv/versions/3.6.7/bin/pip", line 11, in <module>
    sys.exit(main())
TypeError: 'module' object is not callable

解決方法

pyenv環境で次のコマンドでpipインストールを再実行したところ解決。

curl -kL https://bootstrap.pypa.io/get-pip.py | python

あとがき

エラーメッセージで検索すると、aptでpipをアップデートしている記事や、pip install pip==18.0などでversionを戻している記事などはあったのですが、直接的な解法にたどり着くのに時間がかかったのでメモを兼ねて記事にしました。

何番煎じかわかりませんが、(タイムゾーン指定のない)シンプルな日時文字列から、入力と出力のタイムゾーンを指定して変換結果を得るコードのメモになります。

utcとの差分9時間を加減すれば処理できるのですが、ちゃんと?、'UTC'、'Asia/Tokyo'などとタイムゾーンを指定して変換したい場合のコードです。

日付関連の用語に慣れていないので、もし誤用などありましたらご指摘ください。

使っているライブラリ

import datetime
import pytz

JSTの日時文字列をUTCに変換する

# 文字列から(タイムゾーン情報を持たないnaiveな)datetimeオブジェクトを作成する
datetime_obj = datetime.datetime.strptime("2019-10-28 13:00:00", "%Y-%m-%d %H:%M:%S")

# 変換(JST -> UTC)
utc_datetime_obj = pytz.timezone('Asia/Tokyo').localize(datetime_obj).astimezone(pytz.timezone('UTC'))

# 結果を確認
print(utc_datetime_obj.strftime("%Y-%m-%d %H:%M:%S"))
# (9時間前の時刻が表示される)
# 2019-10-28 04:00:00

やっていること

datetimeオブジェクトを作成

日時文字列とその書式情報から、タイムゾーン情報を持たないnaiveなdatetimeオブジェクトを作成

datetime_obj = datetime.datetime.strptime("2019-10-28 13:00:00", "%Y-%m-%d %H:%M:%S")

datetimeオブジェクトにタイムゾーン情報を付与

作成したdatetimeオブジェクトに(JSTの)タイムゾーン情報を付与してawareなdatetimeオブジェクトを作成

jst_datetime_obj = pytz.timezone('Asia/Tokyo').localize(datetime_obj)

astimezoneでUTCに変換

utc_datetime_obj = jst_datetime_obj.astimezone(pytz.timezone('UTC'))

(逆) UTCの日時文字列をJSTに変換する

# この文字列をUTCに見立てて処理する
# 文字列から(タイムゾーン情報を持たないnaiveな)datetimeオブジェクトを作成する
datetime_obj = datetime.datetime.strptime("2019-10-28 13:00:00", "%Y-%m-%d %H:%M:%S")

# 変換(UTC -> JST)
jst_datetime_obj = pytz.timezone('UTC').localize(datetime_obj).astimezone(pytz.timezone('Asia/Tokyo'))

# 確認
print(jst_datetime_obj.strftime("%Y-%m-%d %H:%M:%S"))
# (9時間後の時刻が表示される)
# 2019-10-28 22:00:00

async/awaitを使う

背景

スクレイピングを効率的に行う方法としてmultiprocessingを使ってマルチプロセスを使用する方法がありますが、マルチプロセスでの実行は環境のコア数に依存してしまいます。
リクエストの結果を待ってる時にコアを掴んでおく必要はなく非同期で実行するようにすることで効率的に処理を進める方法を探していたら見つけたので実行方法について書きました。

grequestsとは

geventを使って非同期HTTPリクエストを簡単に実現することができるライブラリです。
geventとは非同期処理をベースとしたネットワーク処理のライブラリです。
bottleやflaskを使用してwebsocketを使う場合によく出てくるものと一緒です。

■ FlaskとWebSocketを使用してリアルタイム通信を行う

■ gevent/gevent

環境

$ uname -a
Darwin mbp01 19.0.0 Darwin Kernel Version 19.0.0: Wed Sep 25 20:18:50 PDT 2019; root:xnu-6153.11.26~2/RELEASE_X86_64 x86_64

$ python3 -V
Python 3.7.4

インストール

grequestsはpypiで公開されているのでpipでインストールできます。

$ pip install grequests

■ grequests 0.4.0

ソースコードはGitHub上で公開されています。

■ spyoungtech/grequests

使ってみる

非同期httpリクエストを実行する方法はとても簡単です。

app.py

import grequests

urls = [
    'http://www.heroku.com',
    'http://python-tablib.org',
    'http://httpbin.org',
    'http://python-requests.org',
    'http://fakedomain/',
    'http://kennethreitz.com'
]

# 非同期リクエスト用のオブジェクトを生成
rs = (grequests.get(u) for u in urls)

# httpリクエスト/結果を表示
print(grequests.map(rs))

サンプルではhttpのGETを実行していますがこれ以外にも基本的には対応しているようです。

grequests.py

# Shortcuts for creating AsyncRequest with appropriate HTTP method
get = partial(AsyncRequest, 'GET')
options = partial(AsyncRequest, 'OPTIONS')
head = partial(AsyncRequest, 'HEAD')
post = partial(AsyncRequest, 'POST')
put = partial(AsyncRequest, 'PUT')
patch = partial(AsyncRequest, 'PATCH')
delete = partial(AsyncRequest, 'DELETE')

さくっと検証

実際に自分の環境で検証してみたのですがスクレイピングプログラムの実行時間が早くなるのは実感できませんでした。
(そもそもそこまでリクエストを多数実行してるわけでもないので。。。)
せっかくなので非同期を体感できるような検証環境を作って試してみました。

サーバアプリケーション

flaskをuwsgiで多重化して実行。
4プロセスで起動するので4リクエストまでは(環境次第で)同時に処理します。
並列性を確認する簡単な方法として非同期sleepである標準モジュールのsleepを使用します。
返却値としてはリクエストの時間とレスポンスの時間をつめて返します。

ちなみにflaskではデフォルトで起動すると複数のリクエストを同時に処理することができず今回の検証ではuwsgiを使ってます。
ただ起動時にthreaded=Trueオプションを指定することでuwsgiを使わなくても検証することは可能です。

While lightweight and easy to use, Flask’s built-in server is not suitable 
for production as it doesn’t scale well and by default serves only one request at a time. 
Some of the options available for properly running Flask in production are documented here.

■ Deployment Options

(threadに関しては実際のシステムではWSGIなどを用いることが多いようなのであまり使われないオプションって認識)

サーバ

#!/usr/local/bin/python3
# coding: utf-8

import datetime
import os
import sys
import time
from flask import Flask, jsonify, request

app = Flask(__name__)


def get_date_formatting():
    return str(datetime.datetime.today())[:-7]


@app.route("/", methods=["GET"])
def hello_world():
    try:
        sec = int(request.args.get("sec"))
    except Exception:
        sec = 0

    req_time = get_date_formatting()
    time.sleep(sec)
    res_time = get_date_formatting()
    return jsonify({
        "pid": os.getpid(),
        "req-time": req_time,
        "res-time": res_time
    })


if __name__ == "__main__":
    app.run(debug=True, host="0.0.0.0", port=5000)

クライアント

同期的に実行
1リクエストごとに処理された時間が追加されていってるのが確認できます。

client.py

import requests

def main():
    urls = [
        "http://localhost:5000?sec=3",
        "http://localhost:5000?sec=5",
        "http://localhost:5000?sec=1",
        "http://localhost:5000"
    ]


    for i in urls:
        print((requests.get(i).text))


if __name__ == "__main__":
    main()

実行結果

{
  "pid": 21577,
  "req-time": "2019-10-28 20:55:22",
  "res-time": "2019-10-28 20:55:25"
}

{
  "pid": 21577,
  "req-time": "2019-10-28 20:55:25",
  "res-time": "2019-10-28 20:55:30"
}

{
  "pid": 21577,
  "req-time": "2019-10-28 20:55:30",
  "res-time": "2019-10-28 20:55:31"
}

{
  "pid": 21577,
  "req-time": "2019-10-28 20:55:31",
  "res-time": "2019-10-28 20:55:31"
}

grequestsを使って実行。
実行時間をみると並列に4リクエスト同時に実行してサーバ側で処理されていることがわかりました。

client_aio.py

import grequests

def main():
    urls = [
        "http://localhost:5000?sec=3",
        "http://localhost:5000?sec=5",
        "http://localhost:5000",
        "http://localhost:5000"
    ]

    rs = (grequests.get(u) for u in urls)

    for r in grequests.map(rs):
        if r is not None:
            print(r.text.rstrip())


if __name__ == "__main__":
    main()

実行結果

{
  "pid": 21577,
  "req-time": "2019-10-28 20:51:31",
  "res-time": "2019-10-28 20:51:34"
}
{
  "pid": 21577,
  "req-time": "2019-10-28 20:51:31",
  "res-time": "2019-10-28 20:51:36"
}
{
  "pid": 21577,
  "req-time": "2019-10-28 20:51:31",
  "res-time": "2019-10-28 20:51:31"
}
{
  "pid": 21577,
  "req-time": "2019-10-28 20:51:31",
  "res-time": "2019-10-28 20:51:31"
}

まとめ

ネットワークIOなどは非同期で処理をしやすいのは間違いないですが、
私が使ってるスクレイピングプログラムだとあまり効果が出ませんでした。大量にリクエストを送るプログラム(テストなど？)の並行実行には強いのではないでしょうか？

はじめに

DeepChemでScaffoldSplitterを普通に使ったところ、掲題のエラーが発生したので調べたときのメモ

環境

• python 3.6
• deepchem 2.2.1.dev54
• rdkit 2019.03.3.0
  

現象

splitterにScaffoldSplitterを指定して意気揚々とCross-Validationをやろうとしたら、deepchem.data.CSVLoaderでcsvファイルを読み込もうとした時に、以下エラーが発生。

  File "C:\Users\XXX\AppData\Local\conda\conda\envs\yyy\lib\site-packages\rdkit\Chem\Scaffolds\MurckoScaffold.py", line
108, in MurckoScaffoldSmiles
    raise ValueError('No molecule provided')
ValueError: No molecule provided

原因と解決策

エラーログを元に色々見たところ、deepchem.splits.splitters.pyの870行目に以下を発見。

 for ind, smiles in enumerate(dataset.ids):
      if ind % log_every_n == 0:
        log("Generating scaffold %d/%d" % (ind, data_len), self.verbose)
      scaffold = generate_scaffold(smiles)
      if scaffold not in scaffolds:
        scaffolds[scaffold] = [ind]
      else:
        scaffolds[scaffold].append(ind)

dataset.idsからsmilesとってんじゃん。ダメじゃん。そこはdataset.smilesからとるべきでしょうが。
といっても回避策はないので、ScaffoldSplitterを使う場合は、CSVLoaderのid_fieldsは指定しないようにしましょう。
そうすると、id_fieldにsmiles_fieldsが設定されるっぽいので。
id_fieldが何に使われてるか知らんけど、まあ今のところ大丈夫っぽい。

loader = dc.data.CSVLoader(tasks=[args.target_col],
#                              id_field=args.id_col, #　ここはコメントオフにしよう。
                               smiles_field=args.smiles_col,
                               featurizer=featurizer)

おわりに

OSSを使う場合は、適度にソースも見ながら、うまく付き合うことが大事ということで。。。

はじめに

DeepChemでScaffoldSplitterを普通に使ったところ、掲題のエラーが発生したので調べたときのメモ

環境

• python 3.6
• deepchem 2.2.1.dev54
• rdkit 2019.03.3.0
  

現象

splitterにScaffoldSplitterを指定して意気揚々とCross-Validationをやろうとしたら、deepchem.data.CSVLoaderでcsvファイルを読み込もうとした時に、以下エラーが発生。

  File "C:\Users\XXX\AppData\Local\conda\conda\envs\yyy\lib\site-packages\rdkit\Chem\Scaffolds\MurckoScaffold.py", line
108, in MurckoScaffoldSmiles
    raise ValueError('No molecule provided')
ValueError: No molecule provided

原因と解決策

エラーログを元に色々見たところ、deepchem.splits.splitters.pyの870行目に以下を発見。

 for ind, smiles in enumerate(dataset.ids):
      if ind % log_every_n == 0:
        log("Generating scaffold %d/%d" % (ind, data_len), self.verbose)
      scaffold = generate_scaffold(smiles)
      if scaffold not in scaffolds:
        scaffolds[scaffold] = [ind]
      else:
        scaffolds[scaffold].append(ind)

dataset.idsからsmilesとってんじゃん。CSVLoaderでid_fieldを明示的に指定した場合は、そこにはsmilesは格納されないわけだから、そのフィールドをsmilesと勝手に判断し、molに変換したらそりゃー、'No molecule provided'が発生するよ。
といっても、こんなファイル修正しようがないのでScaffoldSplitterを使う場合は、CSVLoaderのid_fieldsは指定しないようにしましょう。
そうすると、id_fieldにsmiles_fieldsの値が設定されるっぽいので、結果オーライです。
id_fieldが何に使われてるか知らんけど、まあ今のところ大丈夫っぽい。

loader = dc.data.CSVLoader(tasks=[args.target_col],
#                              id_field=args.id_col, #　ここはコメントオフにしよう。
                               smiles_field=args.smiles_col,
                               featurizer=featurizer)

おわりに

OSSを使う場合は、適度にソースも見ながら、うまく付き合うことが大事ということで。。。　
以上、取り急ぎご報告まで失礼します。

はじめに

この記事はゼロから作るディープラーニング 5章ニューラルネットワークの学習を自分なりに理解して分かりやすくアウトプットしたものです。
文系の自分でも理解することが出来たので、気持ちを楽にして読んでいただけたら幸いです。
また、本書を学習する際に参考にしていただけたらもっと嬉しいです。

学習機能を実装したニューラルネットワーク

今回は分類問題を解く専用の学習機能付きニューラルネットワークのクラスを実装して行きたいと思います。

今回の実装の流れを簡単に説明します。

1,　インスタンス変数でパラメータとハイパーパラメータを実装

2,　予測データを出力するニューラルネットワークの処理を行うクラスメソッドpredictの実装

3,　ニューラルネットワークから損失関数までの処理を行うクラスメソッドlossの実装

4,　予測データの正答率を出すメソッドaccuracyの実装

5,　学習を行うための勾配式を実装

1, インスタンス変数でパラメータとハイパーパラメータを実装

class MNST_net:# 中間層は今回一層

    def __init__(self, input_size, hiden_size, output_size, weight_init_std = 0.01):
        self.params = {}#パラメータの初期化
        self.params['W1'] = weight_init_std * np.random.randn(input_size, hiden_size)
        self.params['b1'] = np.zeros(hiden_size)
        self.params['W2'] = weight_init_std * np.random.randn(hiden_size, output_size)
        self.params['b2'] = np.zeros(output_size)

最初のクラスメソッドでハイパーパラメータとパラメータを実装します。

init(self):のカッコの中には、input_size(入力層のニューロン数),hiden_size(中間層のニューロン数),output_size(出力層のニューロン数)、weight_init_std(重みを小さなな値にして過学習を抑える)などの人間が設定しないといけないハイパーパラメータを記述します。

重みやバイアスなどのパラメータの実装(初期化)はparamsと言う辞書をインスタンス変数として作成し、その中に各層ごと各パラメータの種類ごとに分けて実装していきます。

重みは、np.random.randnを使って(入力信号✖︎ニューロン数)のランダムな整数が入った配列を作成し、weight_init_stdをかけることで実装(初期化)できます。

バイアスは初期値をゼロにしないといけないので、np.zerosでニューロン数のゼロが入った一次元配列を作成することで実装（初期化）できます。

2, 予測データを出力するニューラルネットワークの処理を行うクラスメソッドpredictの実装

class MNST_net:# 中間層は今回一層

    def __init__(self, input_size, hiden_size, output_size, weight_init_std = 0.01):
        self.params = {}#パラメータの初期化
        self.params['W1'] = weight_init_std * np.random.randn(input_size, hiden_size)
        self.params['b1'] = np.zeros(hiden_size)
        self.params['W2'] = weight_init_std * np.random.randn(hiden_size, output_size)
        self.params['b2'] = np.zeros(output_size)

    def predict(self, x):# ニューラルネットワークの処理
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']

        a1 = np.dot(x, W1) + b1# 入力層から隠れ層１段の処理
        z1 = sigmoid_function(a1)

        a2 = np.dot(z1, W2) + b2
        y = softmax_function_pro(a2)
        return y

初期化されたパラメータをそれぞれ変数に入れて、各層のニューラルネットワークのニューロンの処理を実装していきます。

3, ニューラルネットワークから損失関数までの処理を行うクラスメソッドlossの実装

class MNST_net:# 中間層は今回一層

    def __init__(self, input_size, hiden_size, output_size, weight_init_std = 0.01):
        self.params = {}#パラメータの初期化
        self.params['W1'] = weight_init_std * np.random.randn(input_size, hiden_size)
        self.params['b1'] = np.zeros(hiden_size)
        self.params['W2'] = weight_init_std * np.random.randn(hiden_size, output_size)
        self.params['b2'] = np.zeros(output_size)

    def predict(self, x):# ニューラルネットワークの処理
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']

        a1 = np.dot(x, W1) + b1# 入力層から隠れ層１段の処理
        z1 = sigmoid_function(a1)

        a2 = np.dot(z1, W2) + b2
        y = softmax_function_pro(a2)
        return y

    def loss(self, x, t):#  NNから損失関数までの処理
        y = self.predict(x)
        return cross_entropy_errors_label(t, y)

predictメソッドで予測データを取得して変数yに入れ、それを交差エントロピー誤差の引数に入れてreturnで返します。

4, 予測データの正答率を出すメソッドaccuracyの実装

class MNST_net:# 中間層は今回一層

    def __init__(self, input_size, hiden_size, output_size, weight_init_std = 0.01):
        self.params = {}#パラメータの初期化
        self.params['W1'] = weight_init_std * np.random.randn(input_size, hiden_size)
        self.params['b1'] = np.zeros(hiden_size)
        self.params['W2'] = weight_init_std * np.random.randn(hiden_size, output_size)
        self.params['b2'] = np.zeros(output_size)

    def predict(self, x):# ニューラルネットワークの処理
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']

        a1 = np.dot(x, W1) + b1# 入力層から隠れ層１段の処理
        z1 = sigmoid_function(a1)

        a2 = np.dot(z1, W2) + b2
        y = softmax_function_pro(a2)
        return y

    def loss(self, x, t):#  NNから損失関数までの処理
        y = self.predict(x)
        return cross_entropy_errors_label(t, y)

    def accuracy(self, x, t):# 正答率
        y = self.predict(x)
        y = np.argmax(y, axis = 1)# 一番確率の高いラベルのインデックスを出す
#         t = np.argmax(t, axis = 1)# 正解ラベルのインデックスを出す one-hot法

        accuracy = (y == t).sum() / float(x.shape[0])
        return accuracy

predictで予測データを出して、一番値が大きいindexをnp.argmaxで出し、それと正解データを比べて正解したものの合計 / 全体数　で正答率を出します。

5, 学習を行うための勾配式を実装

class MNST_net:# 中間層は今回一層

    def __init__(self, input_size, hiden_size, output_size, weight_init_std = 0.01):
        self.params = {}#パラメータの初期化
        self.params['W1'] = weight_init_std * np.random.randn(input_size, hiden_size)
        self.params['b1'] = np.zeros(hiden_size)
        self.params['W2'] = weight_init_std * np.random.randn(hiden_size, output_size)
        self.params['b2'] = np.zeros(output_size)

    def predict(self, x):# ニューラルネットワークの処理
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']

        a1 = np.dot(x, W1) + b1# 入力層から隠れ層１段の処理
        z1 = sigmoid_function(a1)

        a2 = np.dot(z1, W2) + b2
        y = softmax_function_pro(a2)
        return y

    def loss(self, x, t):#  NNから損失関数までの処理
        y = self.predict(x)
        return cross_entropy_errors_label(t, y)

    def accuracy(self, x, t):# 正答率
        y = self.predict(x)
        y = np.argmax(y, axis = 1)# 一番確率の高いラベルのインデックスを出す
#         t = np.argmax(t, axis = 1)# 正解ラベルのインデックスを出す one-hot法

        accuracy = (y == t).sum() / float(x.shape[0])
        return accuracy

     def slopeing_grad_net(self, x, t):# 初期パラメータ勾配式
        loss_c = lambda W: self.loss(x, t) # 関数化させないといけなかった 引数Wはダミー

        grads = {}
        grads['W1'] = slopeing_grad(loss_c,self.params['W1'])
        grads['b1'] = slopeing_grad(loss_c,self.params['b1'])
        grads['W2'] = slopeing_grad(loss_c,self.params['W2'])
        grads['b2'] = slopeing_grad(loss_c,self.params['b2'])

        return grads

ニューラルネットワークから損失関数までの処理を変数loss_cに入れて、それと初期化されたパラメータを各層・各種類ごとに前回作成したslopeing_gradで勾配を求めて、それをgradsと言う辞書にまとめてreturnで返します。

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

量子情報理論の基本の基として、おさえておくべき事項の中で、２つの量子状態の違いについてのお話が残っていました。よく使われる指標として「フィデリティ(fidelity)」と「トレース距離(trace distance)」というものがあります。今回は、フィデリティについて勉強してみます。どのように定義されるのか、また、どんな性質があるのかということを説明した後で、量子計算シミュレータqlazyを使って、その重要な性質について、実際に計算して実感してみたいと思います。

参考にさせていただいたのは、以下の文献です。

1. ニールセン、チャン「量子コンピュータと量子通信（３）」オーム社（2005年）
2. 石坂、小川、河内、木村、林「量子情報科学入門」共立出版（2012年）
3. 富田「量子情報工学」森北出版（2017年）

フィデリティの定義

内積の絶対値

フィデリティ(fidelity)は、日本語では「忠実度」と訳されますが、要は２つの状態の類似度を表す指標です。純粋状態の類似度として、直感的にわかりやすいのは内積です。２つの純粋状態を各々$\ket{\phi},\ket{\psi}$とすると、内積の絶対値は以下のように表現でき、これをフィデリティ$F(\rho,\sigma)$の定義とします。

F(\ket{\phi},\ket{\psi}) = |\braket{\phi}{\psi}|  \tag{1}

通常、純粋状態はノルム１に正規化されているので、純粋状態のフィデリティはヒルベルト空間上の２つのベクトルのなす角のコサインの絶対値であると解釈できます。したがって、両者直交している場合は０（最小値）となり、平行の場合は１（最大値）となります。

混合状態のフィデリティは、これの拡張版として定義されます。前回の記事で説明した「純粋化」をここで登場させます。任意の混合状態は適当な参照系をアンシラとして追加することで、形式的に純粋状態にすることができるというお話でした。いま、２つの混合状態を$\rho,\sigma$とし、各々を純粋化した状態を$\ket{\phi},\ket{\psi}$とします。その内積の絶対値$|\braket{\phi}{\psi}|$をフィデリティと定義します、と単純に言えれば良いのですが、そう簡単にはいきません。純粋化で得られる状態には参照系に関するユニタリの自由度があり、一意には決まらないのでした。そこで、あり得る純粋状態のすべてを考えて、その最大値をフィデリティの定義とします。つまり、フィデリティ$F(\rho,\sigma)$を

F(\rho,\sigma) = \max_{\ket{\phi},\ket{\psi}} |\braket{\phi}{\psi}|  \tag{2}

のように定義してみます。しかし、$max$が含まれている表式では使い勝手がよろしくないので、何とかしたいと思うわけです。そのため、ちょっと回り道をして、以下の２つのことを確認しておきます。

• ユニタリの自由度を含んだ純粋化の表現
• 演算子のトレースノルム

ユニタリの自由度を含んだ純粋化の表現

まず、純粋化の結果として得られる純粋状態を、ユニタリの自由度を含んだ形に表現し直します。
前回の記事で見たように、注目系Aにおける混合状態、

\rho_{A} = \sum_{k=1}^{d} \lambda_{k} \ket{k}_{A} \bra{k}_{A}  \tag{3}

の純粋化は、参照系Rにおける適当な正規直交基底$\{ \ket{k}_{R} \}$を使って、

\ket{\phi} = \sum_{k=1}^{d} \sqrt{\lambda_{k}} \ket{k}_{A} \ket{k}_{R}  \tag{4}

のように表すことができます。これは、以下のように書くこともできます¹。

\ket{\phi} = (\sqrt{\rho} \otimes I_{R}) \sum_{k=1}^{d} \ket{k}_{A} \ket{k}_{R}  \tag{5}

ところが、ユニタリの自由度があるので、参照系R上で定義されるユニタリ演算子$U_{R}$を用いて、

\ket{\tilde{\phi}} = (\sqrt{\rho} \otimes I_{R}) \sum_{k=1}^{d} \ket{k}_{A} U_{R} \ket{k}_{R}  \tag{6}

と書くこともできます²。

この式は、行列$U_{R}^{T}$を行列$U_{A}$とおき、$U_{A}$を注目系Aでのユニタリ演算子と同一視することで、以下のようにも表現できます。

\ket{\tilde{\phi}} = (\sqrt{\rho} \otimes I_{R}) \sum_{k=1}^{d} U_{A} \ket{k}_{A} \ket{k}_{R}  \tag{7}

これはそれほど自明ではないので、以下で証明してみます。

【証明】

\begin{align}
& \sum_{k=1}^{d} \ket{k}_{A} U_{R} \ket{k}_{R} \\
&= \sum_{k=1}^{d} \sum_{i=1}^{d} \ket{k}_{A} \ket{i}_{R} \bra{i}_{R} U_{R} \ket{k}_{R} \\
&= \sum_{k=1}^{d} \sum_{i=1}^{d} (U_{R})_{ik} \ket{k}_{A} \ket{i}_{R} \\
&= \sum_{k=1}^{d} \sum_{i=1}^{d} (U_{R}^{T})_{ki} \ket{k}_{A} \ket{i}_{R} \\
&= \sum_{i=1}^{d} \sum_{k=1}^{d} (U_{R}^{T})_{ik} \ket{i}_{A} \ket{k}_{R} \\
&= \sum_{i=1}^{d} \sum_{k=1}^{d} \bra{i}_{A} U_{A} \ket{k}_{A} \ket{i}_{A} \ket{k}_{R} \\
&= \sum_{k=1}^{d} U_{A} \ket{k}_{A} \ket{k}_{R}  \tag{8}
\end{align}

したがって、

\begin{align}
\ket{\tilde{\phi}} &= (\sqrt{\rho} \otimes I_{R}) \sum_{k=1}^{d} \ket{k}_{A} U_{R} \ket{k}_{R}  \\
&= (\sqrt{\rho} \otimes I_{R}) \sum_{k=1}^{d} U_{A} \ket{k}_{A} \ket{k}_{R}  \tag{9}
\end{align}

となります。（証明終）

演算子のトレースノルム

次に、回り道の２つ目です。「演算子のトレースノルム」を定義したいのですが、その前に「演算子のノルム（または絶対値）」を定義します。線形演算子$A$のノルム$|A|$は、

|A| \equiv \sqrt{A^{\dagger} A}  \tag{10}

と定義されます。$|A|$と書きますが、行列式ではありません。演算子のノルムは数値ではなく演算子ですので、ご注意ください。$A^{\dagger} A$のスペクトル分解を

A^{\dagger} A = \sum_{i} |\lambda_{i}|^{2} \ket{i} \bra{i}  \tag{11}

とすると、式(10)は、

|A| = \sqrt{A^{\dagger} A} = \sum_{i} |\lambda_{i}| \ket{i} \bra{i}  \tag{12}

と書くことができます。演算子のトレースノルム$||A||$は、これのトレースと定義されます。すなわち、

||A|| \equiv Tr |A| = Tr \sqrt{A^{\dagger} A} = \sum_{i} |\lambda_{i}|  \tag{13}

です。

$A$のトレースノルムは、ユニタリ演算子$V$を用いて、実は、

||A|| = \max_{V} |Tr(AV)|  \tag{14}

のように表すこともできます。以下で証明します。

【証明】

まず、$A$をユニタリ演算子$U$を使って、

A = U|A|  \tag{15}

のように極分解します³。

\begin{align}
Tr(AV) &= Tr(U|A|V) = Tr(|A|VU) \\
&= \sum_{i} \bra{i} \sum_{j} |\lambda_{j}| \ket{j} \bra{j} VU \ket{i} \\
&= \sum_{i} \sum_{j} |\lambda_{j}| \delta_{ij} \bra{j} VU \ket{i} \\
&= \sum_{i} |\lambda_{i}| \bra{i} VU \ket{i} \tag{16}
\end{align}

において、$VU$はユニタリ（つまり、ベクトルの長さを変えない演算）なので、

|\bra{i} VU \ket{i}| \leq 1  \tag{17}

が成り立ちます。したがって、式(16)は、

Tr(AV) \leq \sum_{i} |\lambda_{i}| = ||A||  \tag{18}

となり、

||A|| = \max_{V} |Tr(AV)|  \tag{14}

が証明されました。ちなみに、式(17)からわかるように、$VU = I$のときに最大値になります。（証明終）

Uhlmannの定理

さて、それでは、式(2)の$max$を使わない表式を導出する準備が整いました。注目系Aにおける２つの混合状態を$\rho,\sigma$とし、参照系Rを使った純粋化を各々、以下のように書くことにします。ここで、式(7)の表現を使います。

\begin{align}
\ket{\tilde{\phi}} &= (\sqrt{\rho} \otimes I_{R}) \sum_{k} U_{A} \ket{k}_{A} \ket{k}_{R} \\
\ket{\tilde{\psi}} &= (\sqrt{\sigma} \otimes I_{R}) \sum_{k} U_{A} \ket{k}_{A} \ket{k}_{R} \tag{19}
\end{align}

これを使って、式(2)のフィデリティを計算します。

\begin{align}
F(\rho,\sigma) &= \max_{U_{A},V_{A}} |\braket{\tilde{\phi}}{\tilde{\psi}}| \\
&= \max_{U_{A},V_{A}} |\sum_{k,l} \bra{k}_{A} \bra{k}_{R} U_{A}^{\dagger} \sqrt{\rho} \sqrt{\sigma} V_{A} \ket{l}_{A} \ket{l}_{R}| \\
&= \max_{U_{A},V_{A}} |\sum_{k} \bra{k}_{A} U_{A}^{\dagger} \sqrt{\rho} \sqrt{\sigma} V_{A} \ket{k}_{A}| \\
&= \max_{U_{A},V_{A}} |Tr(U_{A}^{\dagger} \sqrt{\rho} \sqrt{\sigma} V_{A})| \\
&= \max_{U_{A},V_{A}} |\sqrt{\rho} \sqrt{\sigma} W)| \tag{20}
\end{align}

ここで、$W=V_{A} U_{A}^{\dagger}$とおきました。式(14)を使うと、

F(\rho,\sigma) = ||\sqrt{\rho} \sqrt{\sigma}||  \tag{21}

となります。これで、maxを使わないフィデリティの表式が得られました。これを「Uhlmannの定理」と言います⁴。

フィデリティの性質

フィデリティの定義がわかったので、次にいくつかの重要な性質について説明します。

(1) 交換について対称

F(\rho,\sigma) = F(\sigma,\rho)  \tag{22}

が成り立ちます。定義より明らかです。

(2) 同じ状態に対するフィデリティは１

\rho = \sigma \Leftrightarrow F(\rho,\sigma) = 1  \tag{23}

が成り立ちます。密度演算子がエルミートでトレースが１という性質を使えば、簡単にわかります。

【証明】

F(\rho,\rho) = ||\sqrt{\rho} \sqrt{\rho}|| = ||\rho|| = Tr|\rho| = Tr(\rho) = 1 \tag{24}

（証明終）

(3) 直交している状態のフィデリティは０

\rho \sigma = 0 \Leftrightarrow F(\rho,\sigma) = 0  \tag{25}

が成り立ちます。定義より明らかです。

(4) フィデリティの取りうる値は０から１

式(2)の右辺の内積の絶対値は、

$0 \leq |\braket{\phi}{\psi}| \leq 1$なので、フィデリティの取りうる値は０から１です。

(5) 直積状態のフィデリティは積

F(\rho_{1} \otimes \rho_{2},\sigma_{1} \otimes \sigma_{2}) = F(\rho_{1},\sigma_{1}) F(\rho_{2},\sigma_{2})  \tag{26}

が成り立ちます。

【証明】

\begin{align}
F(\rho_{1} \otimes \rho_{2},\sigma_{1} \otimes \sigma_{2})
&= ||(\sqrt{\rho_{1}} \otimes \sqrt{\rho_{2}}) (\sqrt{\sigma_{1}} \otimes \sqrt{\sigma_{2}})|| \\
&= ||\sqrt{\rho_{1}} \sqrt{\sigma_{1}} \otimes \sqrt{\rho_{2}} \sqrt{\sigma_{2}}|| \\
&= Tr |\sqrt{\rho_{1}} \sqrt{\sigma_{1}} \otimes \sqrt{\rho_{2}} \sqrt{\sigma_{2}}| \\
&= Tr (\sqrt{\rho_{1}} \sqrt{\sigma_{1}} \otimes \sqrt{\rho_{2}} \sqrt{\sigma_{2}}) \\
&= Tr (\sqrt{\rho_{1}} \sqrt{\sigma_{1}}) Tr(\sqrt{\rho_{2}} \sqrt{\sigma_{2}}) \\
&= Tr |\sqrt{\rho_{1}} \sqrt{\sigma_{1}}| Tr|\sqrt{\rho_{2}} \sqrt{\sigma_{2}}| \\
&= F(\rho_{1},\sigma_{1}) F(\rho_{2},\sigma_{2})  \tag{27}
\end{align}

（証明終）

(6) あらゆる物理過程でフィデリティは減少しない（単調性または非減少性）

ある物理過程を表す量子チャネルを$\Gamma$とすると、

F(\rho,\sigma) \leq F(\Gamma(\rho),\Gamma(\sigma))  \tag{28}

が成り立ちます⁵。

【証明】

$\rho,\sigma$の純粋化を各々$\ket{\phi},\ket{\psi}$とします。環境系Eの適当な状態として$\ket{0}_{E}$をもってきて追加し（テンソル積をとり）、$U$でユニタリ変換して、環境系をトレースアウトしたものが量子チャネルです。全体系の純粋状態に対する変化は、

\begin{align}
\ket{\phi} \ket{0}_{E} & \rightarrow U \ket{\phi} \ket{0}_{E} \\
\ket{\psi} \ket{0}_{E} & \rightarrow U \ket{\psi} \ket{0}_{E} \tag{29}
\end{align}

であり、フィデリティは純粋化した状態の内積の絶対値の最大値だったので、

\begin{align}
F(\Gamma(\rho), \Gamma(\sigma)) &\geq |\bra{\phi} \bra{0}_{E} U^{\dagger} U \ket{\psi} \ket{0}_{E} | \\
&= |\braket{\phi}{\psi}| = F(\rho,\sigma) \tag{30}
\end{align}

です。（証明終）

(7) フィデリティの強凹性

F(\sum_{i} p_{i} \rho_{i}, \sum_{i} q_{i} \sigma_{i}) \geq \sum_{i} \sqrt{p_{i} q_{i}} F(\rho_{i},\sigma_{i})  \tag{31}

が成り立ちます⁶。

【証明】

$\ket{\phi_i},\ket{\psi_i}$は、$F(\rho_{i},\sigma_{i})=\braket{\phi_{i}}{\psi_{i}}$を満たすように選んだ$\rho_{i},\sigma_{i}$の純粋化とします。適当な正規直交系$\{ \ket{i} \}$をもってきて、以下の状態を定義します。

\begin{align}
\ket{\phi} &= \sum_{i} \sqrt{p_{i}} \ket{\phi_{i}} \ket{i} \\
\ket{\psi} &= \sum_{i} \sqrt{q_{i}} \ket{\psi_{i}} \ket{i} \tag{32}
\end{align}

ここで、$\ket{\phi},\ket{\psi}$は各々$\sum_{i} p_{i} \rho_{i}$および$\sum_{i} q_{i} \sigma_{i}$の純粋化なので、

\begin{align}
F(\sum_{i} p_{i} \rho_{i}, \sum_{i} q_{i} \sigma_{i}) &\geq |\braket{\phi}{\psi}| \\
&= \sum_{i} \sqrt{p_{i} q_{i}} \braket{\phi}{\psi} \\
&= \sum_{i} \sqrt{p_{i} q_{i}} F(\rho_{i},\sigma_{i})  \tag{33}
\end{align}

です。（証明終）

シミュレータで確認

それでは、上で示したフィデリティの性質のうち、６番目の「あらゆる物理過程で減少しない」という性質に注目し、それが本当なのかどうかをシミュレータで確認してみます。具体的には、２つの密度演算子をランダムに作成して、ランダムに作った量子チャネル（参照系＋環境系を追加して純粋化した状態に対してランダムなユニタリ変換を実施し最後にトレースアウトする、というやり方で定義しました）を通した結果、フィデリティが確かに減少しない（＝増加する）ということを見てみます。

全体のPythonコードは以下です。

import random
import numpy as np
from scipy.stats import unitary_group
from qlazypy import QState, DensOp

def random_densop(qnum_tar,qnum_ref,qnum_env):

    dim_pur = 2**(qnum_tar+qnum_ref)
    vec_pur = np.array([0.0]*dim_pur)
    vec_pur[0] = 1.0
    mat_pur = unitary_group.rvs(dim_pur)
    vec_pur = np.dot(mat_pur, vec_pur)

    dim_env = 2**qnum_env
    vec_env = np.array([0.0]*dim_env)
    vec_env[0] = 1.0

    vec_whole = np.kron(vec_pur,vec_env)

    qs = QState(vector=vec_whole)
    de = DensOp(qstate=[qs],prob=[1.0])

    qs.free()
    return de

def random_unitary(qnum):

    dim = 2**qnum
    mat = unitary_group.rvs(dim)

    return mat

if __name__ == '__main__':

    # settings
    qnum_tar = 2  # system A : target system
    qnum_ref = 2  # system R : reference system
    qnum_env = 2  # system E : environment system
    qnum_whole = qnum_tar + qnum_ref + qnum_env

    # two random states in system A+R+E (A+R:set randomly, E:set |0> initialy)
    de1_whole = random_densop(qnum_tar,qnum_ref,qnum_env)
    de2_whole = random_densop(qnum_tar,qnum_ref,qnum_env)

    # two states in system A (trace out R+E)
    de1_ini = de1_whole.partial(id=list(range(qnum_tar)))
    de2_ini = de2_whole.partial(id=list(range(qnum_tar)))

    # fidelity for initial states
    fid_ini = de1_ini.fidelity(de2_ini)

    # unitary transformation for whole system
    U = random_unitary(qnum_whole)
    de1_whole.apply(U)
    de2_whole.apply(U)

    # two states in system A (trace out R+E)
    de1_fin = de1_whole.partial(id=list(range(qnum_tar)))
    de2_fin = de2_whole.partial(id=list(range(qnum_tar)))

    # fidelity for final states
    fid_fin = de1_fin.fidelity(de2_fin)

    # result
    print("* fidelity(ini) =", fid_ini)
    print("* fidelity(fin) =", fid_fin)

    if fid_ini < fid_fin:
        print("OK!")
    else:
        print("NG!")

    # free memory
    de1_whole.free()
    de2_whole.free()
    de1_ini.free()
    de2_ini.free()
    de1_fin.free()
    de2_fin.free()

何をやっているか、順に説明します。

# settings
qnum_tar = 2  # system A : target system
qnum_ref = 2  # system R : reference system
qnum_env = 2  # system E : environment system
qnum_whole = qnum_tar + qnum_ref + qnum_env

まず、量子ビット数を適当に設定します。密度演算子が定義されている注目系Aの量子ビット数を２、それを純粋化するための参照系Rの量子ビットを２、注目系Aの周囲を取り巻く環境系Eの量子ビットを２とします。

# two random states in system A+R+E (A+R:set randomly, E:set |0> initialy)
de1_whole = random_densop(qnum_tar,qnum_ref,qnum_env)
de2_whole = random_densop(qnum_tar,qnum_ref,qnum_env)

全体系に対する純粋状態をランダムに作ります。random_densop関数で実行しています。その中身を見てみます。

dim_pur = 2**(qnum_tar+qnum_ref)
vec_pur = np.array([0.0]*dim_pur)
vec_pur[0] = 1.0
mat_pur = unitary_group.rvs(dim_pur)
vec_pur = np.dot(mat_pur, vec_pur)

dim_env = 2**qnum_env
vec_env = np.array([0.0]*dim_env)
vec_env[0] = 1.0

vec_whole = np.kron(vec_pur,vec_env)

qs = QState(vector=vec_whole)
de = DensOp(qstate=[qs],prob=[1.0])

最初の４行で、注目系Aと参照系Rを合わせた系で定義される $\ket{0}^A\ket{0}^R$ に相当するベクトルを作成し、scipyの関数で作ったランダムなユニタリ演算を実施して、ランダムな純粋状態を作り、変数vec_purに格納します。注目系Aの密度演算子をまず作って純粋化すべきところですが、純粋化したものをランダムに用意することにしました（すみませせん、処理を端折りました。同じことなので）。次に、環境系Eに関して、$\ket{0}^{E}$に相当するベクトルを作成し、変数vec_envに格納します。２つのベクトルのテンソル積（クロネッカー積）をnumpyの関数を使って実行し、変数vec_wholeに格納します。これで、注目系でのランダムな密度演算子に対応した、全体系の純粋状態ができたことになります。これに基づき、密度演算子のインスタンスdeを生成し、リターンします。

main部に戻ります。

# two states in system A (trace out R+E)
de1_ini = de1_whole.partial(id=list(range(qnum_tar)))
de2_ini = de2_whole.partial(id=list(range(qnum_tar)))

全体系の密度演算子から参照系と環境系をトレースアウトして、注目系に注目します。密度演算子に対するpatialメソッドに引数として注目系の量子番号リストを与えることで実行できます。今の場合、引数には[0,1]が入ることになります。これで、２つのランダムな初期密度関数が得られたので、、、

# fidelity for initial states
fid_ini = de1_ini.fidelity(de2_ini)

で、本日のメインイベントであるフィデリティを計算します。qlazypyの密度演算子クラスDensOpにfidelityメソッドを追加しました(v0.0.27)。この１行で終わりです。中ではnumpy.linlagを使っていて、固有値問題を解いて得られた固有値のルートを計算するとか、トレースを計算するとか、やっています。

# unitary transformation for whole system
U = random_unitary(qnum_whole)
de1_whole.apply(U)
de2_whole.apply(U)

random_unitary関数でランダムなユニタリ行列を取得します。関数定義を見ていただければわかる通り、scipyの関数を召喚しています。このユニタリを全体系の密度演算子に施すため、密度演算子クラスDensOpのapplyメソッドを使っています。これで、初期状態がランダムに定義された量子チャネルを通ったことになります。

# two states in system A (trace out R+E)
de1_fin = de1_whole.partial(id=list(range(qnum_tar)))
de2_fin = de2_whole.partial(id=list(range(qnum_tar)))

量子チャネルを通った後の状態に対して、部分トレースをとり、注目系での密度演算子を取得します。

# fidelity for final states
fid_fin = de1_fin.fidelity(de2_fin)

先ほどと同様、フィデリティを計算します。

# result
print("* fidelity(ini) =", fid_ini)
print("* fidelity(fin) =", fid_fin)

if fid_ini <= fid_fin:
    print("OK!")
else:
    print("NG!")

結果を表示します。フィデリティが増加していれば「OK!」、そうでなければ「NG!」を表示します。

さて、実行結果です。

* fidelity(ini) = 0.7344321087768828
* fidelity(fin) = 0.9503724941444357
OK!

というわけで、フィデリティが増加していることがわかります。何度も実行しましたが必ず増加しました。

おわりに

「はじめに」で述べましたが、量子状態の違いを定量化する指標として、今回説明した「フィデリティ」の他に、「トレース距離」というものもあります。これは名前の通り、状態がどれだけ離れているかを表す指標です。なので「フィデリティ」とはちょうど逆の関係にあります。この指標もよく使われるということなので、次回は「トレース距離」について勉強したいと思います。

以上

---

1. $\rho=\sum_{k} \lambda_{k} \ket{k}^{A} \bra{k}^{A}$から、容易にわかると思います。 ↩

2. これは大丈夫ですよね。ユニタリの自由度があるということを素直に表しているだけです。 ↩

3. 任意の線形演算子$A$に対して、あるユニタリ演算子（正確には等距離演算子）$U$があって、$A = U|A|$と表すことができるという定理があります。詳細は量子情報科学入門の付録をご参照ください。直感的には、$A^{\dagger} A$を計算してみれば、確かにそうだということがわかると思います。任意の複素数$z$を$z=|z|e^{i\theta}$と表すことができますが、それに少し似ています。 ↩

4. ニールセン、チャンと量子情報科学入門では、式(21)をフィデリティの定義として式(2)が成り立つことをもって「Uhlmannの定理」という説明がされていました。今回の記事では、量子情報工学の議論進行に従い、逆の流れで証明しました。その方が直感的にわかりやすいと思ったからです。 ↩

5. 最初にあった識別可能な特徴をもった２つの状態が、環境系と相互作用しながら次第に平衡状態に達して識別不能になっていくという風に考えれば良いと思います。一瞬、大小関係が逆じゃね？と思った方は速攻誤解を正しましょう（自分のことです、汗）。 ↩

6. なぜこれを「凹性」というのか？$\space 0 \leq p \leq 1$となる$p$に対し、関数$f(x)$が$f(px_{1} + (1-p)x_{2}) \leq pf(x_{1}) + (1-p)f(x_{2})$という性質を満たすとき、その関数$f$を凸関数と言い（関数のカーブを思い描くと上に凸ですよね）、不等号が逆の場合、凹関数と言います（こっちは逆に下に凸ですね、すなわち上に凹）。このアナロジーで理解しておけば良いと思います。ここでは「強凹性」のみを取り上げましたが、「強」がつかない「凹性」もあります。ニールセン、チャンにいくつかのバリエーションが紹介されています。 ↩

shiracamus さん

pythonの勉強をしていく過程で、Qiitaもはじめました。なにか投稿しようと思い、python勉強のために作成したものを投稿します。

はじめに

最初に数字当てゲームを知ったのは
｢退屈なことはPythonにやらせよう ―ノンプログラマーにもできる自動化処理プログラミング｣
という本でした。

この本には、1から20までの数字を6回までに当てるゲームのコードが、載っていました。

random.randint()で選ばれた数字を、int ( input () )で入力してもらう。というものです。

これをもとに、数字の範囲を広げました。(+9981!!)
が、しかし、当たるわけがありません。

そこで、自分なりにヒントのコードを加えて、当てやすくしたり、時間計測を追加したりして、ゲームとしてなんとか成り立つようにしました。

初心者なので細かくコメントを入れています。(本来は適切ではない……)

1．ヒントの表示

数字入力時に外した場合

その数字より、大きいか、小さいかを返します。
これは本にも載っていましたが、数字が近かった場合(200以内)、表示を変えるようにしました。

# 解答者が選んだ数字が大きいか、小さいか返す関数
def high_low():
    # 数字が200以内なら文言が変わる
    if (guess < secret_number) and (guess + 200 >= secret_number):
        return ("\nもう少し、大きい数字です。(200以内)")
    if (guess > secret_number) and (guess <= secret_number + 200):
        return ("\nもう少し、小さい数字です。(200以内)")
    # 数字が大きいか、小さいかの文言
    if guess < secret_number:
        return ("\nもっと、大きい数字です。")
    if guess > secret_number:
        return ("\nもっと、小さい数字です。")

数字を外して、｢ヒントを使いますか?｣で｢1:はい｣を選んだ場合

外した回数によって、ヒントの内容が変わります。

･ 1回外した場合 → 偶数か奇数か
･ 2回外した場合 → 3で割れるか、5で割れるか(Fizz Buzz)
･ 3回外した場合 → 各桁の数字を足す
･ 4回外した場合 → 下2桁目を表示
･ 5回外した場合 → ±10以内の数字をランダムで表示

｢2:いいえ｣を選んだ場合は、ヒントは表示されません。

# 各ヒントをまとめた関数
def hint(select):
    # 1回目で外した場合、偶数か奇数かを表示
    if guesses_taken == 1:
        print ( "\n偶数か奇数かを表示する……" )
        # 2で割り、余りがなかったら、True
        if not (secret_number % 2):
            return ("…… 偶数!!")
        # 2で割り、余りがあったら、True
        if (secret_number % 2):
            return ("…… 奇数!")

    # 2回目で外した場合、3で割れるか、5で割れるか、を表示
    if guesses_taken == 2:
        print ( "\n3で割れるか、5で割れるか、表示する……" )
        # 3で割っても5で割っても余りがない
        if not (secret_number % 3) and not (secret_number % 5):
            return ("…… 両方で割れる")
        # 3で割ったら余りがない
        if not (secret_number % 3):
            return ("…… 3で割り切れる(5で割れない)")
        # 5で割ったら余りがない
        if not (secret_number % 5):
            return ("…… 5で割り切れる(3で割れない)")
        return ("…… どちらでも割れない")

    # 3回目で外した場合、各桁の数字を表足す
    if guesses_taken == 3:
        # 答えの数字を文字列に変えて、dに格納
        # dにある文字列を、1文字ずつ数字に変えリスト化
        digits = [int ( d ) for d in str ( secret_number )]
        return ("\n各桁の数字を足す……"
                "\n…… " + str ( sum ( digits ) ))  # 各数字を足して文字列型にする

    # 4回目で外した場合、下2桁目を表示
    if guesses_taken == 4:
        return ("\n下2桁目を表示する……"
                "\n…… " + str ( secret_number )[-2])  # 後ろから2文字目を取得

    # 5回目で外した場合、±10以内のランダム数字を表示
    if guesses_taken == 5:
        return ("\n±10以内の数字をランダムで表示する……"
                "\n…… " + str ( random.randint ( secret_number - 10, secret_number + 10 ) ))

2．計測時間の表示と、使用したヒントの表示

調べたら、time.perf_counter ()が精密らしいので実装しました。

start_time = time.perf_counter ()
"""コード"""
end_time = time.perf_counter ()

# 計測した時間の計算(秒数)
tim = end_time - start_time
# 秒数から分数のみ(秒数なし)を計算
mint = tim // 60
# 秒数から秒数のみ(分数なし)を計算
secd = (tim % 60) - mint
# 分数は小数点以下を非表示(計算の時点で必ず0のため)、秒数は小数点第一位まで表示
print ( "所要時間は {:.0f} 分 {:.1f} 秒".format ( mint, secd ) )

｢ヒントを使いますか?｣で｢1:はい｣を選んだ場合に、空の使用リストに格納して表示

# 使用したヒントを表示するための記録用、空リスト
record = []

 # ヒントを、順番に列挙したタプルを用意
    used_ability = ("偶数か奇数か", "3で割れるか、5で割れるか", "各桁の数字を足す", "下2桁目を表示", "±10以内の数字をランダムで表示")
    print ( "\nヒントを使いますか?\n 1:はい / 2:いいえ" )
    yesno = int ( input () )
    # yesnoが1だった場合
    if yesno == 1:
        # 解答数に応じた能力を使う
        ability = hint ( guesses_taken )
        print ( ability )
        # 能力のタプルから、解答回数に応じた能力名を、onebyoneに格納
        onebyone = used_ability[guesses_taken - 1]
        # 順番にrecord = []に格納にしてリスト化
        record.append ( onebyone )
    # yesnoが1以外だった場合pass
    if not yesno == 1:
        pass

# 使用したヒントを表示
print ( "\n使用したヒント" )
print ( record )

コード

guessTheNumber.py

# randomをインポート
import random
# timeをインポート
import time


# 解答者が選んだ数字が大きいか、小さいか返す関数
def high_low():
    # 数字が200以内なら文言が変わる
    if (guess < secret_number) and (guess + 200 >= secret_number):
        return ("\nもう少し、大きい数字です。(200以内)")
    if (guess > secret_number) and (guess <= secret_number + 200):
        return ("\nもう少し、小さい数字です。(200以内)")
    # 数字が大きいか、小さいかの文言
    if guess < secret_number:
        return ("\nもっと、大きい数字です。")
    if guess > secret_number:
        return ("\nもっと、小さい数字です。")


# 各ヒントをまとめた関数
def hint(select):
    # 1回目で外した場合、偶数か奇数かを表示
    if guesses_taken == 1:
        print ( "\n偶数か奇数かを表示する……" )
        # 2で割り、余りがなかったら、True
        if not (secret_number % 2):
            return ("…… 偶数!!")
        # 2で割り、余りがあったら、True
        if (secret_number % 2):
            return ("…… 奇数!")

    # 2回目で外した場合、3で割れるか、5で割れるか、を表示
    if guesses_taken == 2:
        print ( "\n3で割れるか、5で割れるか、表示する……" )
        # 3で割っても5で割っても余りがない
        if not (secret_number % 3) and not (secret_number % 5):
            return ("…… 両方で割れる")
        # 3で割ったら余りがない
        if not (secret_number % 3):
            return ("…… 3で割り切れる(5で割れない)")
        # 5で割ったら余りがない
        if not (secret_number % 5):
            return ("…… 5で割り切れる(3で割れない)")
        return ("…… どちらでも割れない")

    # 3回目で外した場合、各桁の数字を表足す
    if guesses_taken == 3:
        # 答えの数字を文字列に変えて、dに格納
        # dにある文字列を、1文字ずつ数字に変えリスト化
        digits = [int ( d ) for d in str ( secret_number )]
        return ("\n各桁の数字を足す……"
                "\n…… " + str ( sum ( digits ) ))  # 各数字を足して文字列型にする

    # 4回目で外した場合、下2桁目を表示
    if guesses_taken == 4:
        return ("\n下2桁目を表示する……"
                "\n…… " + str ( secret_number )[-2])  # 後ろから2文字目を取得

    # 5回目で外した場合、±10以内のランダム数字を表示
    if guesses_taken == 5:
        return ("\n±10以内の数字をランダムで表示する……"
                "\n…… " + str ( random.randint ( secret_number - 10, secret_number + 10 ) ))


# 使用したヒントを表示するための記録用、空リスト
record = []

# 0 ~ 10,000までの数字をランダムに決め、secret_numberに格納
secret_number = random.randint ( 0, 10000 )
print ( "0から10,000までの数字を当ててください。" )

# 時間の計測開始
start_time = time.perf_counter ()

# 最大6回繰り返す
for guesses_taken in range ( 1, 7 ):
    print ( "\n数字を入力してください。" )
    # 入力してもらった数字をguessに格納
    # input()は文字列を返すので、int()で整数値に変換
    guess = int ( input () )

    if guesses_taken == 6:
        break

    # 数字が当たればbreakでforループから抜け出す
    if guess == secret_number:
        break

    print ( high_low () )

    # ヒントを、順番に列挙したタプルを用意
    used_ability = ("偶数か奇数か", "3で割れるか、5で割れるか", "各桁の数字を足す", "下2桁目を表示", "±10以内の数字をランダムで表示")
    print ( "\nヒントを使いますか?\n 1:はい / 2:いいえ" )
    yesno = int ( input () )
    # yesnoが1だった場合
    if yesno == 1:
        # 解答数に応じた能力を使う
        ability = hint ( guesses_taken )
        print ( ability )
        # 能力のタプルから、解答回数に応じた能力名を、onebyoneに格納
        onebyone = used_ability[guesses_taken - 1]
        # 順番にrecord = []に格納にしてリスト化
        record.append ( onebyone )
    # yesnoが1以外だった場合pass
    if not yesno == 1:
        pass

# 時間の計測終了
end_time = time.perf_counter ()

# 当たった場合とハズレた場合で文言を変える
if guess == secret_number:
    print ( "\nグッド! " + str ( '{:,}'.format ( guesses_taken ) ) + "回で当たり!" )
if guess != secret_number:
    print ( "\n結果だけを求めていると、人は近道をしたがるものだ……………"
            "\n近道した時、真実を見失うかもしれない"
            "\nやる気もしだいに失せていく"
            "\n大切なのは『真実に向かおうとする意志』だと思っている"
            "\n……正解は " + str ( '{:,}'.format ( secret_number ) ) )

# 使用したヒントを表示
print ( "\n使用したヒント" )
print ( record )

# 計測した時間の計算(秒数)
tim = end_time - start_time
# 秒数から分数のみ(秒数なし)を計算
mint = tim // 60
# 秒数から秒数のみ(分数なし)を計算
secd = (tim % 60) - mint
# 分数は小数点以下を非表示(計算の時点で必ず0のため)、秒数は小数点第一位まで表示
print ( "所要時間は {:.0f} 分 {:.1f} 秒".format ( mint, secd ) )

おわりに

ちょっとした気持ちで、範囲を0から10,000にしたら、あれもつけたい、これもつけたいで、想像以上に時間を取られました。

本にも書いてありましたが、必勝法は二分探索(binary search)
二分探索で迫れば、最後のヒントを見なくても、結構当てられます。

pythonもQiitaも、まずは1歩

前提

このシリーズは、以下の環境で動作確認などを行っています。環境によっては細かい部分が異なることがありますが、わからないことはコメントへどうぞ。

環境

OS: Windows 10 Pro (1809, 17763.737)
CPU: Intel(R) Core(TM) i7-8500Y@1.50GHz x4
RAM: 8GB

準備するもの

• テキストエディタ
  • Visual Studio Code(https://code.visualstudio.com/)
• Python
  • Python3.8.0(https://www.python.org/downloads/release/python-380/)
• ライブラリ
  • OpenCV(pip install opencv-python)
  • 他に必要なのがあったら追記します

インストールとか準備

Visual Studio Code

自分の環境にあったものを落としましょう。今回は最新版の1.39.2 を使います。
ダウンロードしたら、"VSCodeUserSetup-x64-1.39.2.exe"みたいなファイルがあると思います。（環境によってちょっとずつ違います）これを実行しましょう。そして、ダイアログに従ってインストールします。

プラグイン

拡張機能です。便利なので入れときましょう。（というか必要）
赤枠のところをポチッとして、検索すれば入ります。

日本語パック - https://marketplace.visualstudio.com/items?itemName=MS-CEINTL.vscode-language-pack-ja
Python - https://marketplace.visualstudio.com/items?itemName=ms-python.python
Indent-Rainbow - https://marketplace.visualstudio.com/items?itemName=oderwat.indent-rainbow
Bracket Pair Colorizer 2 (Beta) - https://marketplace.visualstudio.com/items?itemName=CoenraadS.bracket-pair-colorizer-2

Python3.8.0

これも同じです。自分の環境のものを落として実行します。僕の場合は"python-3.8.0-amd64.exe"でした。
そして、コマンドプロンプトで動作確認。

実行結果

C:\Users\xxx> py
 Python 3.8.0 (tags/v3.8.0:fa919fd, Oct 14 2019, 19:37:50) [MSC v.1916 64 bit (AMD64)] on win32
 Type "help", "copyright", "credits" or "license" for more information.
>>>

OSが違っても、表示される内容はほぼ同じです。バージョンがあっていて起動できればOK。
(Linux, MacOSXはpython3だったりpythonだったりします。)

OpenCV

py -m pip install opencv-python
以上。かーんたん♪

動作確認

C:\Users\xxx> py
 Python 3.8.0 (tags/v3.8.0:fa919fd, Oct 14 2019, 19:37:50) [MSC v.1916 64 bit (AMD64)] on win32
 Type "help", "copyright", "credits" or "license" for more information.
>>> import cv2
>>>

何も表示されなければOK。

次回

次回は、顔認証の学習に必要な画像を作るためのプログラムです。お楽しみに～

GitHub

自分で全部コードを打つのはちょっと・・・
というめんどくさがり屋さんのためにGitHubにソースなど一式を上げときます。
記事の更新とほぼ同時にあげます。
Hiro527/OpenCV-Py3-Face - https://github.com/Hiro527/OpenCV-Py3-Face