乙Py先生のプログラミング教室
初学者のためのプログラミング学習サイト

機械学習

今回は仮想通貨の自動売買ロジックの
組み立て方のチュートリアルコードを
動かしてみました。


解説動画はこちら



はじめに


今回は書籍
「日給300万円のSS級トレーダーが明かすbotterのリアル」
のチュートリアルコードを動かす方法などについてです。


書籍画像
No description has been provided for this image


チュートリアルのコードはこちら
githubのコード

こちらのコードは
そのままでは動かない可能性もあるので
それ含めて動かせる様にしてみました。


チュートリアル環境の構築


1.Dockerを用いる方
Readme通りにやるだけです。
# コードのクローン
git clone https://github.com/richmanbtc/mlbot_tutorial.git

# Dockerの起動
cd mlbot_tutorial
docker-compose up -d

# Jupyterの起動
# http://localhost:8888 をブラウザで開く

2.自分のJupyter環境 or Colabでやる方
必須ライブラリをインストールしておく必要があります。
# Jupyterの人は全部
ccxt : 仮想通貨取引用のライブラリ
TA-lib : 指標の計算ライブラリ

# Colabの人は不要
numba : Pythonの高速化ライブラリ
scikit-learn : 機械学習用ライブラリ
lightGBM : 機械学習のモデル構築ライブラリ

なおインストール方法は割愛します。
インストールが自力で進める方だけ
進んでください。


チュートリアルの進め方


ここからはチュートリアルのコードを
そのまま進める形になります。

チュートリアルのコードを
そのまま実行してみて下さい。

チュートリアルはこのような
内容になっていました。

1.ライブラリのインポート
2.データ取得
3.データ読み込み
4.前処理
5.目的変数の計算
6.モデルの学習
7.バックテスト


このうち
2.データ取得
3.データ読み込み
に関しては環境によっては
うまく行えない可能性がありました。

ライブラリの差異などにより
うまくデータが生成されませんでした。

なので代わりのコードを
置いておきます。

データ取得のコード
import urllib.request
import os
import time
from datetime import datetime, timedelta
# SSLの問題が有ったら追加
import ssl
ssl._create_default_https_context = ssl._create_unverified_context

def get_date_range(start_date, end_date):
    start = datetime.strptime(start_date, '%Y-%m-%d')
    end = datetime.strptime(end_date, '%Y-%m-%d')
    date_range, current_date= [],start
    while current_date <= end:
        date_range.append((current_date.year, current_date.month, current_date.day))
        current_date += timedelta(days=1)
    return date_range

def get_data(market,start_date,end_date,data_dir = "data"):
    data_dir = f"{data_dir}/{market}"
    url_base = 'https://api.coin.z.com/data/trades/{0}/{1}/{2:02}/{1}{2:02}{3:02}_{0}.csv.gz'
    dates = get_date_range(start_date, end_date)
    for d in dates:
        year,month,day = d
        url = url_base.format(
            market,
            year,
            month,
            day
        )
        file_name = os.path.basename(url)
        if not os.path.exists(f"{data_dir}/{year}"):
            os.makedirs(f"{data_dir}/{year}")
        save_path = os.path.join(f"{data_dir}/{year}", file_name)
        urllib.request.urlretrieve(url, save_path)
        time.sleep(1.37)
# マーケットと開始日、終了日を指定
market = "BTC_JPY"
start_date = "2024-01-01"
end_date = "2024-02-22"
get_data(market,start_date,end_date)

実行するとこのコードが置いてある
ディレクトリ配下にdataが配置されます。

データの読み込み
import os
import glob

def make_df(file_path,interval_sec):
    df = pd.read_csv(file_path)
    df = df.rename(columns={'symbol': 'market',})
    df['price'] = df['price'].astype('float64')
    df['size'] = df['size'].astype('float64')
    df['timestamp'] = pd.to_datetime(df['timestamp'], utc=True)
    df['side'] = np.where(df['side'] == 'BUY', 1, -1).astype('int8')
    df['timestamp'] = df['timestamp'].dt.floor('{}S'.format(interval_sec))
    df.set_index("timestamp", inplace=True)
    volume = df.groupby('timestamp')['size'].sum().rename('volume')
    ohlc = df['price'].resample('1T').ohlc()
    df_merged = ohlc.merge(volume, left_index=True, right_index=True,how="left")
    df_merged = df_merged.fillna(method="ffill")
    df_merged = df_merged.rename(columns={"open":"op","high":"hi","low":"lo","close":"cl"})
    return df_merged
# データの保存場所を指定
data_dir = "data/BTC_JPY/2024/"
interval_sec = 60
df = pd.DataFrame()

for file_path in sorted(glob.glob(data_dir + "*.gz")):
    tmp_df = make_df(file_path,60)
    df = pd.concat([df,tmp_df],axis=0)

こちらのコードを実行すると
変数 df が生成されます。

チュートリアルで用いられている
データを読み込んだ変数になっているので
そのまま使用する事ができると思います。

データが読み込めたかどうかは
次のコードで確認してみて下さい。
df[["cl"]].plot(figsize=(12,4))
plt.show()
download


前処理部分のコードの追加

コードの前処理部分の最初で
手数料を記入するコードがあります。

チュートリアルが作成されてからだいぶ経っており
手数料などが更新されていないため
適宜追加してみて下さい。

maker_fee_history = [
    {
        # https://coin.z.com/jp/news/2020/08/6482/
        # 変更時刻が記載されていないが、定期メンテナンス後と仮定
        'changed_at': '2020/08/05 06:00:00Z',
        'maker_fee': -0.00035
    },
    {
        # https://coin.z.com/jp/news/2020/08/6541/
        'changed_at': '2020/09/09 06:00:00Z',
        'maker_fee': -0.00025
    },
    {
        # https://coin.z.com/jp/news/2020/10/6686/
        'changed_at': '2020/11/04 06:00:00Z',
        'maker_fee': 0.0
    },

    ### 追加
    {
        # 現在値
        'changed_at': '2023/08/05 06:00:00Z',
        'maker_fee': -0.0003
    },
]

あとはチュートリアルを
上から実行していけば
動いていくだろうと思います。

累積リターンのところで
次のような画像が出ればほぼほぼ成功してると思います。

download-1



特徴量の重要度の算出


一般的な機械学習では
予測にどの変数が寄与したのか
重要度を算出することが多いです。

こちらがチュートリアルの
モデルの重要度を算出するコードになります。
feature_importance = model.feature_importances_

# 特徴量名と重要度を紐づける
feature_importance_df = pd.DataFrame({'Feature': features, 'Importance': feature_importance})

# 重要度の降順でソート
feature_importance_df = feature_importance_df.sort_values(by='Importance', ascending=False).reset_index(drop=True)

# 結果の表示(トップ10)
feature_importance_df.head(10)
Feature Importance
0 HT_DCPERIOD 145
1 ADXR 137
2 BBANDS_upperband 125
3 KAMA 124
4 BBANDS_lowerband 116
5 LINEARREG 114
6 MFI 109
7 MACD_macdsignal 107
8 HT_PHASOR_quadrature 106
9 HT_DCPHASE 102


この意味合いとしては
chatGPTによると次の様な意味になっています。



HT_DCPERIOD: Hilbert Transform - Dominant Cycle Period
(ヒルベルト変換 - 優勢サイクル期間)と呼ばれる指標
時系列データの優勢周期を表します。この特徴量は
トレンドのサイクル性を捉えるために使用されます。


ADXR: Average Directional Movement Index Rating
(平均方向運動指数レーティング)は
トレンドの強さを測定する指標

特に、ADX(Average Directional Index)に基づいて計算され
トレンドの強さを示します。


KAMA: Kaufman's Adaptive Moving Average
(カウフマンの適応的移動平均)は移動平均線の一種
価格変動の速さに応じて重みを調整し
より滑らかな移動平均線を提供します。

BBANDS_upperbandおよびBBANDS_lowerband: Bollinger Bands
(ボリンジャーバンド)は、移動平均線を中心に上下に
標準偏差を加えたバンドをプロットする指標

上側バンドと下側バンドは、価格の上下の範囲を示し
価格の過熱やサポート/レジスタンスのレベルを示すのに使用されます。


LINEARREG: Linear Regression(線形回帰)は
価格や指標の値の線形トレンドを推定するための手法です。

この特徴量は、過去のデータからの線形回帰に基づいて
将来のトレンドを予測するために使用されます。


HT_PHASOR_quadratureおよびHT_PHASOR_inphase:
Hilbert Transform - Phasor Components
(ヒルベルト変換 - フェーザーコンポーネント)は
サイン波の位相成分と直交成分を表す指標です。

これらの特徴量はサイクル性や位相関係を捉えるために使用されます。

MFI: Money Flow Index(マネーフローインデックス)は
取引量と価格の関係を用いて買い圧力と売り圧力を計算する指標
主に過買いや過売りの状態を示すのに使用されます。


BETA: ベータは、株式のリスクを示す指標で
市場全体(代表的な指数など)との相関関係を指します。

ベータが1より大きい場合、株価の変動が市場の変動よりも
大きいことを示し、逆にベータが1より小さい場合はその逆を意味します。


ULTOSC: Ultimate Oscillator(アルティメットオシレータ)は
短期、中期、長期の3つの期間を用いてトレンドの強さを計算する指標
主に過買いや過売りの状態を示すのに使用されます。


HT_DCPHASE: Hilbert Transform - Dominant Cycle Phase
(ヒルベルト変換 - 優勢サイクル位相)は
優勢サイクルの位相を表す指標です。この特徴量は
トレンドの位相情報を捉えるために使用されます。


STOCHF_fastk: Stochastic Fast %K
(ストキャスティクスファスト%K)は
価格の変動範囲内での位置を示す指標
主に過買いや過売りの状態を示すのに使用されます。


おまけ

自分もチュートリアルを参考に
別のモデルを構築してみました。

説明変数は同一
目的変数を5分後の価格
に設定し

学習期間は
予測行の前期間step分のみ(step=30)
というモデルにし

1分ずつスライドして
学習しては予測を行うようにしました。

これが予測結果です。

実測と予測値の差分のプロット
download-3

実測と予測値の差分のヒストグラム
download-2

実測と予測のプロット
download-4


学習と予測を近いところで行っているので
ズレすぎることは少ないですね。

最後に予測が現在価格を上まっている条件で
5分後価格と現在価格の差分の累計をプロット
download-5

見事に右肩上がり
予測はかなり機能している様に見えました。


最後に

あくまでシミュレーション上ですが
予測が機能している様に見えるので
売買のロジックに組み込むことも出来るかもしれません。

簡易なロジックとしては
5分後価格予測が現在価格を上まっていたら買い
5分後に売る

5分後価格予測が現在価格を下まっていたら売り
5分後に買う

というような感じですね。

仮の予測モデルは出来ているので
あとは売買ロジック周りを構築すれば
仮想通貨の自動売買botが出来上がります。

その辺りも、もう少し進めていけたらと
思っています。

今回は仮想通貨の自動売買bot本の
チュートリアルを動かしてみたでした。

それでは。

今回はPython言語で機械学習用のUIが
サクサク作れるGradioについてです


解説動画はこちら






Gradioについて

機械学習用のモデルを一般の人が使うには
いちいちHTMLやJavascriptとかを書いて
適切なUIを実装しないと使えません


Gradio を使用すると
すべてPythonでデモを構築して
共有するところまでできちゃいます



導入方法

コマンドラインで
インストールコマンドを実行するだけです
pip install gradio



主なGradio操作方法

1.どこか作業用のディレクトリを用意する
2.ディレクトリ内にGradioのコードを書いた
 Pythonファイルを作る(app.py)
3.ターミナルなどからコマンドで実行する
 gradio ファイル名



サンプル起動用のファイル

デモサイトのソースです
sample.py

import gradio as gr

def greet(name):
    return "Hello " + name + "!"

demo = gr.Interface(fn=greet, inputs="text", outputs="text")

ファイルを作成したらターミナルを起動して
ファイルの有る場所に移動し
次のコマンドを実行します

gradio sample.py



アプリケーションが起動したら
以下のURLにアクセスすると見る事ができます

http://127.0.0.1:7860 or 7861
http://localhost:7860 or 7861

ポート番号は起動時にターミナルに表示されるので
それを入力しましょう

スクリーンショット 2023-07-08 16.57.27

こんな感じの画面が出たら成功です



Gradioの仕組み

デモソース
import gradio as gr

def greet(name):
    return "Hello " + name + "!"

demo = gr.Interface(fn=greet, inputs="text", outputs="text")
gradioをimportして
gr.Interface でインターフェイスを作成します


画面に表示されるのがインターフェイスです
スクリーンショット 2023-07-08 16.57.27


インターフェイスは3つの必須パラメータで
初期化を行います

1 . fn: UIをラップする関数

2 . inputs: 入力に使用するコンポーネント
(例"text": 、"image"または"audio")

3 . outputs: 出力に使用するコンポーネント
(例"text"、"image"、 または"label")

関数の引数や戻り値は
対応する入力ソースの数や
出力するソースの数に合わせる必要があります

入力と出力、そして何をするかを関数で定義して
gr.Interfaceの引数に指定して実行するだけ
という非常にシンプルな作りになっています



コンポーネントについて

GradioではUI変更用にさまざまな部品が用意されていて
これを使用してさまざまなUIを作成することができます

部品のサンプルソースです
import gradio as gr

def greet(name, is_morning, temperature):
    salutation = "Good morning" if is_morning else "Good evening"
    greeting = f"{salutation} {name}. It is {temperature} degrees today"
    celsius = (temperature - 32) * 5 / 9
    return greeting, round(celsius, 2)

demo = gr.Interface(
    fn=greet,
    inputs=["text", "checkbox", gr.Slider(0, 100)],
    outputs=["text", "number"],
)
demo.launch()
スクリーンショット 2023-07-08 17.01.22

入力として
・テキスト
・チェックボタン
・スライダー
が配置されています

部品を配置してUIをサクサク作る事ができます

これらを簡単に作成出来るのも
Gradioの強みですね




サンプルアプリ

huggingfaceの画像判別モデルを使用して
画像判別を行うUIサンプルアプリです

今回は
huggingface/google/vit-base-patch16-224
というモデルを使用する事にします

サンプルアプリのソース
import gradio as gr

demo = gr.Interface.load(
    "huggingface/google/vit-base-patch16-224",
    examples=[
    "./images/ファイル名",
    "./images/ファイル名",
    ])
    
demo.launch()

事前に画像ファイルを用意しておき
imagesというディレクトリ の中に配置しておきます

examples の部分のファイル名を
正しく指定しておきましょう

こんな感じの画面が出ます


スクリーンショット 2023-07-08 17.05.54

画像ファイルをクリックして「送信」をクリックすると
画像の判別結果が右側に出る仕組みです

このように判別モデルなんかも
数行のコードで使う事ができてしまうUIが
簡単に作成できてしまうので
非常に便利です

作ったコードを
Google colabなどで配布すれば
簡単に共有も出来てしまいます


huggingfaceにモデルを置いておけば
どこでも使う事ができるので
オープンソースで公開を検討している場合は
非常にお手軽に使う事ができます



今回は機械学習用のUIが
サクサク作れてしまう
Gradio をご紹介しました

それでは

今回はたった2行で
機械学習の精度検証が出来ちゃう
Lazypredictを試してみました

解説動画はこちら


このライブラリは
機械学習用のモデルを自動作成して
複数のモデルを比較してくれます

これがgithubです
Lazypredict

ドキュメントとしては
それほど情報量が多くないようです

コードなどもシンプルで
かなり使いやすい印象ですね


GoogleColabで動かせるので
使いたい方はコードを参考にしてみてください。

まず最初はライブラリのインストールです
# ライブラリのインストール
pip install lazypredict
30秒くらいでインストールも終わります。

GoogleColabなら追加のインストールが
無くても動くと思います。

このライブラリは
数値を予測する回帰モデルと
種別を予測する判別モデル
2種類の機械学習モデルの比較を行えます



最初は判別モデルをみてみましょう
データはscikit-learnのデータセットを使います

load_wine
ワインの品種に関するデータセットで
14列178個のデータがあり
3クラス分類で予測します


説明変数
alcohol アルコール濃度
malic_acid リンゴ酸
ash 灰
alcalinity_of_ash 灰のアルカリ成分
magnesium マグネシウム
total_phenols 総フェノール類量
flavanoids フラボノイド(ポリフェノールらしい)
nonflavanoid_phenols 非フラボノイドフェノール類
proanthocyanins プロアントシアニジン(ポリフェノールの一種らしい)
color_intensity 色の強さ
hue 色合い
od280/od315_of_diluted_wines ワインの希釈度合い
proline プロリン(アミノ酸の一種らしい)

目的変数
ワインの品種

こんな感じのデータです
読み込みのコードはこちら
from sklearn import datasets
import pandas as pd
df = datasets.load_wine(as_frame=True).frame

print('行列数 : ' , df.shape)
print('ターゲットの種別数 : ' , df['target'].nunique())
df.head()

これでデータが用意できたので
次は学習データの分割です
トレーニング7 , テストサイズ3で分割します
# 学習データの分割
from sklearn.model_selection import train_test_split
X = df.iloc[0: , 0:-1]
Y = df['target']

# トレーニング7 , テストサイズ3で分割
x_train, x_test, y_train, y_test = train_test_split(X , Y , test_size=0.3 , random_state=0)
x_train.shape, x_test.shape, y_train.shape, y_test.shape
((124, 13), (54, 13), (124,), (54,))


最後にlazypredictで
予測モデルの作成です

分類モデルを複数自動で生成し
その精度検証の結果と
予測結果を出してくれます

モデルの作成、検証部分は2行でいけます
分類モデルの検証はLazyClassifierを使います
# lazypredictのインポート
import lazypredict
# 回帰予測
#from lazypredict.Supervised import LazyRegressor
# 分類問題
from lazypredict.Supervised import LazyClassifier

# 予測モデルの作成(たったの2行だけ)
clf = LazyClassifier(verbose=1,ignore_warnings=True,predictions=True)
models , predictions = clf.fit(x_train , x_test , y_train , y_test)

最後の行の変数
models , predictions
に精度検証の結果と
予測結果が格納される仕組みです


作成した予測モデルの比較を見てみましょう

# モデルの比較
models
スクリーンショット 2022-10-01 16.01.17

こんな感じで予測モデルを一気に検証してくれます


テストデータの予測結果も見てみましょう
# テストデータ で予測
predictions
スクリーンショット 2022-10-01 16.02.36

列方向は各モデルで
行がテストデータです


回帰モデルも同様に試してみましょう

load_boston
ボストン住宅価格のデータです
from sklearn import datasets
import pandas as pd

data = datasets.load_boston()
boston_df = pd.DataFrame(data=data["data"] , columns=data["feature_names"])
boston_df['MEDV'] = data["target"]

print('行列数 : ' , boston_df.shape)
print('ターゲットの平均値 : ' , boston_df['MEDV'].mean())
boston_df.head()

学習データの分割も同様です
# 学習データの分割
from sklearn.model_selection import train_test_split
X = boston_df.iloc[0: , 0:-1]
Y = boston_df['MEDV']

# トレーニング7 , テストサイズ3で分割
x_train, x_test, y_train, y_test = train_test_split(X , Y , test_size=0.3 , random_state=0)
x_train.shape, x_test.shape, y_train.shape, y_test.shape
回帰モデルの検証も2行でいけますが
回帰の場合はLazyRegressorを使用します

# lazypredictのインポート
import lazypredict
# 回帰予測
from lazypredict.Supervised import LazyRegressor
# 分類問題
#from lazypredict.Supervised import LazyClassifier

# 予測モデルの作成(たったの2行だけ)
reg = LazyRegressor(verbose=1,ignore_warnings=True,predictions=True)
models , predictions = reg.fit(x_train , x_test , y_train , y_test)
# モデルの比較
models
スクリーンショット 2022-10-01 15.42.01

こんな感じの結果になりました



まとめ


データを突っ込めばどの手法の精度が高くなるのかを
結構な速さで出してくれるので
初めて予測モデルを作るデータには
最適かもしれません


AutoMLで似たようなものに
PyCaretというのが有りますが
Lazypredictの方がバグが出なくて
良いかもしれないなーと思いました

検証をパパッとやりたい方には
かなり良いライブラリかと思いますので
試したい方はコードを参考にしていただければと思います

それでは

先日、食べログ関係の裁判結果が
ニュースになっていたので
気になって調べちゃいました。

解説動画はこちら 


さて、最近こんな判決がニュースになっていました。

東京地裁は16日
グルメサイト「食べログ」の運営会社カカクコムに対し
3840万円の損害賠償金の支払いを命じる判決を出した。

気になってしまったので
二年前に作った
食べログの星を予測するモデルを
引っ張り出してきました。

当時の記事はこちら

さて、データも昔のデータが有ったので
そのまま流用してモデル作成です。

RandomForest
MAE :  0.08949504442045006
RMSE :  0.1153333553457383

download

精度としては
3.2と予測したら
実際が3.1 - 3.3くらいに落ち着くような
結果を出すモデルになります。

たまに大きく外れますが
まあ、誤差でしょう・・・。

これを踏まえて
今回の裁判を起こした側のお店のデータを取って
このモデルにかけてみました。

すると
スクリーンショット 2022-06-25 17.45.13

平均で0.18
最大0.4ほどずれる結果となりました。

予測よりも実測が全体的に低くなる
結果となりました。

download-1

誤差を見てみると
高い点を予測した方が
ズレが大きいですねー

download-2

予測モデルの検証では誤差が少ないものが
最頻値でしたが、この店舗のデータだと
結構外れていますね。


ここからは個人的な推測になりますが
この問題は一方的な意見のみを
見るべきではないと思っています。

採点アルゴリズムというのは
いわば企業秘密であり
秘伝のタレもようなものです
コーラの元みたいな。

他者に教える事は出来ないものであるので
いわば生命線となります。

これにいちゃもんつけて来た訳になるので
この際は白黒はっきりしてもらいたい訳です。

この手のアルゴリズムには
不正防止の観点から
有る程度の秘匿性を持って
運営を行う必要があります。

攻略できてしまったら
不正に点を釣り上げようとすることも
出来てしまう訳なので。

というわけで一方的に
食べログが悪いという判断はできません。

考えられる事象としては
・何らかの理由で故意に下げた
・不正の排除の結果下がった

というのが可能性としてあるのかなと思います。



おいしいお店はおいしく採点されるべきだし
おいしくないお店は
低くなっていないとおかしいです。

お店や食べログの利益は二の次
最優先に守るべきは
使うユーザー側の利益です。

ですので、裁判では
この点を白黒はっきりしてもらいたいですね。

今後もおもしろくなっていきそうです。

今回はこれまでです
それでは。



今回はプログラミングにおける
距離の測り方をまとめてみました。

解説動画はこちら




さて今回はプログラミングなどで用いられる
距離の測り方についてです。

wikipediaでは
距離(きょり、distance)
ある2点間に対して測定した長さの量をいう
ということになっていますが
距離と言っても結構たくさんの距離が存在します。


今回扱う距離は次のようなものです。

1.数学的な距離
 ユークリッド距離
 マンハッタン距離
 チェビシェフ距離
 ミンコフスキー距離
 マハラノビス距離

2.文字列の類似度
 ハミング距離
 レーベンシュタイン距離
 ジャロ・ウィンクラー距離

3.数学/統計学/機械学習における類似度
 コサイン類似度
 ピアソンの相関係数

4.集合同士の類似度
 ジャッカード係数
 ダイス係数
 シンプソン係数



それでは一つ一つの距離の測り方を
見ていくこととしましょう。

Python言語では
scipy.spatial ライブラリの distanceを用いると
簡単に距離を求める事ができるようになっています。

始めにライブラリを読み込みしておきましょう。
from scipy.spatial import distance
import numpy as np

import matplotlib.pyplot as plt
%matplotlib inline




1.数学的な距離

・ユークリッド距離

まずはじめは一般的に距離と言うと
思い浮かばれるのかこれだと思います。

2次元では、2点間の直線距離のことで
各データの差の二乗の合計値のルート
が距離になります。

スクリーンショット 2022-04-02 17.27.59
直角三角形の斜辺を求めるやり方と一緒ですね。
図表を作図して2点を表示してみましょう。

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(0,0,marker='o',c='blue')
plt.plot(4,3,marker='o',c='red')
plt.plot([0,4],[0,3],linestyle='dashed',c='green')
plt.title('Euclid Distance')
plt.axis('square')
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.show()
download-4

これがユークリッド距離になります。

これを求めるには次のコードで
求めることができます。
from scipy.spatial import distance

a = [0,0]
b = [4,3]

# ユークリッド距離を求める
euclid_distance = distance.euclidean(a, b)
print(euclid_distance)
5.0



3次元の場合は数値を増やすことで対応できます。
from scipy.spatial import distance

a = [0,0,0]
b = [1,1,1]

# ユークリッド距離を求める
euclid_distance = distance.euclidean(a, b)
print(euclid_distance)
1.7320508075688772



・マンハッタン距離

マンハッタンや京都のような
碁盤の目のような街を移動する時の距離です。

それぞれの次元の差の合計値ですね
スクリーンショット 2022-04-02 17.28.07

作図はこんな感じです。
import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(0,0,marker='o',c='blue')
plt.plot(4,3,marker='o',c='red')
plt.plot([0,4],[0,0],linestyle='dashed',c='green')
plt.plot([4,4],[0,3],linestyle='dashed',c='green')
plt.title('Manhattan Distance')
plt.axis('square')
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.show()
download-3

次のコードで距離を求めることができます。
import numpy as np

a = np.array([0, 0])
b = np.array([4, 3])

# マンハッタン距離を求める
manhattan_distance = np.linalg.norm(a - b , ord=1)
print(manhattan_distance)
7.0


次元が増えても対応できます。
import numpy as np

a = np.array([0, 0, 0])
b = np.array([4, 4, 4])

# マンハッタン距離を求める
manhattan_distance = np.linalg.norm(a - b , ord=1)
print(manhattan_distance)
12.0




・チェビシェフ距離

次元の差の絶対値の最大値のことで
座標次元ごとで求めた差の中で
一番大きいものになります。

スクリーンショット 2022-04-02 17.28.13

作図はこんな感じですが
マンハッタン距離と同じです。
import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(0,0,marker='o',c='blue')
plt.plot(4,3,marker='o',c='red')
plt.plot([0,4],[0,0],linestyle='dashed',c='green')
plt.plot([4,4],[0,3],linestyle='dashed',c='green')
plt.title('Chebyshev Distance')
plt.axis('square')
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.show()
download-2

距離は次のコードで求めることができます。
from scipy.spatial import distance

a = [0 , 0]
b = [4 , 3]

# チェビシェフ距離を求める
chebyshev_distance = distance.chebyshev(a , b)
print(chebyshev_distance)
4

次元が増えても対応できます。
from scipy.spatial import distance

a = np.array([0, 0, 0])
b = np.array([1, 2, 4])

# チェビシェフ距離を求める
chebyshev_distance = distance.chebyshev(a , b)
print(chebyshev_distance)
4


・ミンコフスキー距離

ユークリッド、マンハッタン、チェビシェフ距離を
一般化したようなもので、pによって取りうる値が変わります。

p=1にするとマンハッタン距離
p=2にするとユークリッド距
p=∞にするとチェビシェフ距離
になります。


コードはこんな感じです。
from scipy.spatial import distance

a = [0 , 0]
b = [4 , 3]

# ミンコフスキー距離を求める p=1
minkowski_distance = distance.minkowski(a,b,p=1)
print(minkowski_distance)

# ミンコフスキー距離を求める p=2
minkowski_distance = distance.minkowski(a,b,p=2)
print(minkowski_distance)

# ミンコフスキー距離を求める p=99
minkowski_distance = distance.minkowski(a,b,p=99)
print(minkowski_distance)

pの値を変えることで距離が変化します。



・マハラノビス距離

一般的にあまりなじみのない距離の名前かと思いますが
多変数間の相関関係を取り入れて
既知のサンプルとの関係を明らかにする距離です

データ群の平均や共分散行列を用いる方法で
異常検知などに使われるケースが多いようです。

サンプルデータを作図してみましょう。

# サンプルデータの可視化
import numpy as np
import matplotlib.pyplot as plt

x1 = np.array([(x, x + np.random.normal()) for x in np.arange(-5, 5, 0.1)])
x2 = np.array([(x, 12 + x + np.random.normal()) for x in np.arange(-10, 0, 0.1)])

plt.plot(x1[:,0], x1[:,1] , '.' , c='blue')
plt.plot(x2[:,0], x2[:,1] , '.' , c='orange')
plt.plot(0 , 0 , marker='o' , c='red',markersize=10)
plt.title('Mahalanobis Distance')
plt.axis('square')
plt.ylim(-20, 20)
plt.xlim(-20, 20)
plt.show()
download

青いデータが正常なもの
オレンジが異常データだと仮定します。

これにユークリッド距離で
等高線を描いてみましょう。
# ユークリッド距離で見てみる
import numpy as np
from scipy.spatial import distance
import matplotlib.pyplot as plt

# サンプルデータ作成
x1 = np.array([(x, x + np.random.normal()) for x in np.arange(-5, 5, 0.1)])
x2 = np.array([(x, 12 + x + np.random.normal()) for x in np.arange(-10, 0, 0.1)])

# ユークリッド距離を計算する
m = x1[:,0].mean(), x1[:,1].mean()
y = np.linspace(-20, 20, 100)
z = np.array([[(i, j, distance.euclidean([i, j], m)) for i in y] for j in y])

# 等高線を表示する
z1,z2,z3 = z.transpose()[0],z.transpose()[1], z.transpose()[2]
cont = plt.contour(z1,z2,z3, levels=[5,10,15], colors=['k'])
cont.clabel(fmt=' d=%1s', fontsize=12)

plt.plot(x1[:,0], x1[:,1] , '.' , c='blue')
plt.plot(x2[:,0], x2[:,1] , '.' , c='orange')
plt.plot(0 , 0 , marker='o' , c='red',markersize=10)
plt.title('Mahalanobis Distance')
plt.axis('square')
plt.ylim(-20, 20)
plt.xlim(-20, 20)
plt.show()
download-1

正常な青いデータを囲もうとすると
ユークリッド距離だとなかなかうまくいきません。

今度は青いデータを用いて
マハラビノス距離を考えてみましょう。

次のコードは簡単なデータでの
マハラビノス距離の計算です。
import numpy as np
from scipy.spatial import distance

# サンプルを用意
x = np.array([[0, 0], [1, 1], [2, 2]])

# サンプルの分散共分散行列の逆行列を計算する
vcm = np.linalg.pinv(np.cov(x.T))

#  [0 , 0]と[3 , 3]のマハラノビス距離を計算する
mahalanobis_distance = distance.mahalanobis([0, 0] , [3, 3], vcm)
print(mahalanobis_distance)
2.9999999999999996

マハラビノス距離は2点のデータと
分散共分散行列のデータが必要になります。


先ほどのサンプルデータに
マハラビノス距離で等高線を描いてみましょう。

# マハラノビス距離で見てみる
import numpy as np
from scipy.spatial import distance
import matplotlib.pyplot as plt

# サンプルデータ作成
x1 = np.array([(x, x + np.random.normal()) for x in np.arange(-5, 5, 0.1)])
x2 = np.array([(x, 12 + x + np.random.normal()) for x in np.arange(-10, 0, 0.1)])

# マハラノビス距離を計算する
m = x1[:,0].mean(), x1[:,1].mean()
vcm = np.linalg.pinv(np.cov(x1.T))
y = np.linspace(-20, 20, 100)
z = np.array([[(i, j, distance.mahalanobis([i, j], m,vcm)) for i in y] for j in y])

# 等高線を表示する
z1,z2,z3 = z.transpose()[0],z.transpose()[1], z.transpose()[2]
cont = plt.contour(z1,z2,z3, levels=[5,10,15], colors=['k'])
cont.clabel(fmt=' d=%1s', fontsize=12)

plt.plot(x1[:,0], x1[:,1] , '.' , c='blue')
plt.plot(x2[:,0], x2[:,1] , '.' , c='orange')
plt.plot(0 , 0 , marker='o' , c='red',markersize=10)
plt.title('Mahalanobis Distance')
plt.axis('square')
plt.ylim(-20, 20)
plt.xlim(-20, 20)
plt.show()
download-1

今度は青いデータとオレンジのデータを
マハラビノス距離でうまく分離することができます。

こんな感じで距離を計算すれば
異常検知などに使えますね。



2.文字列の類似度

・ハミング距離

同じ文字数の二つの文字列を比較したとき
同じ位置にある異なる文字の個数になります。

二つの文字列を同じくするために
必要な文字の置換の回数(0回 - MAXは文字数分)

ABCDE と ABDFC なら3回(割合は0.6)
割合のMAXは1.0になります。
from scipy.spatial import distance

# 文字列をリストとして格納
s1 = list('abcde')
s2 = list('abdfc')

# ハミング距離を計算する(同じ長さであること)
hamming_distance = distance.hamming(s1,s2)
print(hamming_distance)

hamming_num = hamming_distance * len(s1)
print(hamming_num)
0.6
3.0



・レーベンシュタイン距離

二つの文字列がどの程度異なっているかを示す
距離の一種です。

1文字の挿入・削除・置換によって
一方の文字列をもう一方の文字列に変形するのに
必要な手順の最小回数になります。

編集距離とも呼ばれ、文字列の類似度が分かります。


例:kitten と sitting
sitten(「k」を「s」に置換)
sittin(「e」を「i」に置換)
sitting(「g」を挿入して終了)

3回の手順が必要なので
レーベンシュタイン距離は3

Pythonでは専用のライブラリが有るので
インストールしておきましょう。
!pip install python-Levenshtein

距離の計算は次のようになります。
import Levenshtein

s1 = 'abcde'
s2 = 'abdfcee'

# レーベンシュタイン距離を算出する
levenshtein_distance = Levenshtein.distance(s1, s2)
print(levenshtein_distance)
3



・ジャロ・ウィンクラー距離

ある文字列と別の文字列で一致する文字数と
置換の要不要から計算する距離になります。

距離の取りうる値が0~1で
距離の数値が大きいほど
文字列間の類似度が高くなります。


import Levenshtein

s1 = 'abcde'
s2 = 'abdfc'

# ジャロ・ウィンクラー距離を算出する
jaro_winkler_distance = Levenshtein.jaro_winkler(s1, s2)
print(jaro_winkler_distance)
0.8266666666666667

これらの距離を使い分けることで
単語の類似度などに
応用することができます。



3.数学/統計学/機械学習における類似度

・コサイン類似度

2つのベクトルが「どのくらい似ているか」という
類似性を表す尺度のことで
2つのベクトルがなす角のコサイン値になります。

1なら完全に似ている
0なら無関係
-1なら完全に似ていない
ということになります。

作図するとこんな感じです。
import matplotlib.pyplot as plt
from matplotlib import patches

# 2つのベクトルを指定
a = np.array([5.0, 5.0])
b = np.array([3.0, 9.0])

# 作図
fig, ax = plt.subplots()
plt.quiver(0, 0, a[0], a[1], angles='xy', scale_units='xy', scale=1, color='c', label='a')
plt.quiver(0, 0, b[0], b[1], angles='xy', scale_units='xy', scale=1, color='orange', label='b')
c = patches.Arc(xy=(0.9, 1.5), width=0.5, height=0.3, angle=-60, 
                theta1=0, theta2=180, linewidth=3, color="r")
ax.add_patch(c)
plt.title('Cosine Similarity')
plt.axis('square')
plt.xlim(-2,10)
plt.ylim(-2,10)
plt.legend()
plt.grid()
plt.show()
download

計算は次のコードでできます。
from scipy.spatial import distance

a = [2, 48, 9, 11]
b = [3, 51, 13, 19]

cosine_dos = 1 - distance.cosine(a , b)
print(cosine_dos)
0.9902460327205821

真逆の方向を向くベクトルで計算してみましょう。
from scipy.spatial import distance

a = [1 , 1]
b = [-1 , -1]

cosine_dos = 1 - distance.cosine(a , b)
print(cosine_dos)
-1.0




・ピアソンの相関係数

2つの変数間の関係の強さと方向性を表す
-1 ~ 1の範囲の数値になります。

1(強い正の相関)では、2つの変数が強く同方向に連動
-1(強い負の相関)では強く逆方向に連動

データ群を用意して計算してみましょう。
from scipy.spatial import distance

a = [2,3,4,5]
b = [2,3,5,6]

correlation_coefficient =  1 - distance.correlation(a , b)
print(correlation_coefficient)
0.9899494936611665

from scipy.spatial import distance

a = [2,3,4,5]
b = [5,4,3,1]

correlation_coefficient =  1 - distance.correlation(a , b)
print(correlation_coefficient)
-0.9827076298239907

相関係数は様々なライブラリで
求めることができるので
好きな方法を使うのが良いと思います。



4.集合同士の類似度

・ジャッカード係数

集合の類似度を表す指標でテキストマイニングでは
文章と文章の類似度=距離を表す指標になります。

2つの集合に含まれている要素のうち
共通要素が占める割合です。

0から1の間の値となり値が大きいほど
2つの集合の類似度は高い(似ている)ことになります。
スクリーンショット 2022-04-02 18.03.38

集合の計算方法で係数を求めることができます。

import numpy as np

# ジャッカード係数を算出する
def jaccard_coefficient(x,y):
    # 積集合を求める
    intersection = len(set.intersection(set(x), set(y)))    
    # 和集合を求める
    union = len(set.union(set(x), set(y)))                
    return intersection / union

a = np.array([1, 2, 3, 5])
b = np.array([1, 2, 3, 4, 6])

jaccard_index = jaccard_coefficient(a,b)
print(jaccard_index)
0.5



・ダイス係数

2つの集合の平均要素数と
共通要素数での割合になります。

スクリーンショット 2022-04-02 18.03.43


import numpy as np

# ダイス係数を算出する
def dice_coefficient(x,y):
    num = len(set.intersection(set(x), set(y)))
    if num!=0:
        return (2.0 * num) / (len(x) + len(y))
    else:
        return 0.0
    
a = np.array([1, 2, 3, 5])
b = np.array([1, 2, 3, 4, 6])

dice_index = dice_coefficient(a,b)
print(dice_index)
0.6666666666666666


・シンプソン係数

ダイス係数の定義式の分母を「2集合の平均要素数」から
「少ない方の要素数」にしたものです。

ダイス係数よりも差集合の要素数による影響を下げ
相対的に共通要素数を重視した類似度計算です。


スクリーンショット 2022-04-02 18.03.48

# Simpson係数
def simpson_coefficient(x,y):
    num = len(set.intersection(set(x), set(y)))
    if num!=0:
        return num / min([len(x),len(y)])
    else:
        return 0.0

a = np.array([1, 2, 3, 5])
b = np.array([1, 2, 3, 4, 6])

simpson_index = simpson_coefficient(a,b)
print(simpson_index)
0.75

集合の一部が共通要素の場合は
Jaccard係数 < Dice係数 < Simpson係数となります。

これらの係数は
 文書同士の類似度計算
 トピック同士の関連推定
 文字列同士の共起推定
など類似性や距離を測るのに用いられます。



まとめ

プログラミングでは文字列や集合
ベクトルの類似度など
距離を測るのにたくさんの方法があるので
目的に応じて使い分けられるようになると
すごく捗るでしょう。

今回はこれまでです
それでは。

このページのトップヘ