作る前に予測したいと皆思う

有機材料のPL(フォトルミネッセンス)波長をビシッと予測してみましょう。有機材料の中でもベンゼン環などが入っている芳香族化合物は光が当たると、そのエネルギーを吸収して、材料特有の発光を出します。最近市場に出ている有機ELテレビやiPhoneやGalaxyといったスマホ、ゲーム向けの高精細ディスプレイなどのパネルには有機材料の発光が使われています(正確にはエレクトロルミネッセンス)。

材料の構造を変化させると、青、緑、赤と発光色を変化させることができ、このチューニングが分子設計・合成の醍醐味でもあります。もちろん、光るだけではダメで、長い間光るか(安定性)もかなり重要です。

写真) 有機EL発光材料の図(https://newswitch.jp/p/26151)

話すときりがないので背景の説明はここまでとし、手順を説明していきます。以下の手順で進めていきましょう。

データセットの準備
データ分布の可視化、記述子計算
予測モデルの構築、ハイパーパラメータ調整
未知の材料の予測
予測結果の妥当性判断
予測モデルの解釈(SHAP分析)

といった流れで進めていきます。図にするとこんな感じです。今回の記事では①②を進めていきます。

この流れができるようになると材料を合成する前に、事前にPL波長を予測することができ、合成する材料に優先順位をつけることで、材料開発のスピードアップ、無駄になってしまう資源を減らすことができます(有機材料合成では、試薬、溶媒、装置、人件費とかなりコストがかかります)。

もちろん、基本の流れを知っていればPL材料だけでなく他の特性予測にも使えるので、興味ある方ぜひぜひ応用してみてください。それでは始めていきましょう！！

1. データセットの準備

↓実際に使うデータセットを共有しておきます。
https://drive.google.com/file/d/1JIJH0XqFydZKm8t4x1DoHDcFqIKFUir-/view?usp=sharing

オープンアクセスの論文や特許をいろいろ調べて、構造とPL波長を引っ張ってきました。下記のような構成です。

学習データにはPLという列に発光波長(単位はnm)が記入されており、未知のデータのPL列は空欄です。未知のデータなので正解がわからないからです。参考までに青が450-460nm, 緑は520-530nm, 赤は630-640nmぐらいです。

Typeという列は、可視化の際に材料の中心骨格ごとに分けて可視化したかったので用意していますが、予測モデルを作るときには使わないので、なくても大丈夫です。

SMILES列は少し複雑なので細かく説明します。材料の構造、例えばベンゼン環や原子と原子間の結合などはそのままではPCは理解することができないので、PCが取り扱いできるように別の形で変換してあげる必要があります。その手法の一つがSMILESです。文字として、分子の構造を1行の並びで表現できます。SMILESを取得するには、まず構造を何かしらのツールで描写する必要があり、大学や会社ではChemDrawという優秀なツールがあるのですが、個人利用ではお金がかかり使いにくいです。そこで無料でできる方法として、私はMarvinSketchをお勧めします。

登録は必要なのですが、登録後はMarvin Demoというインターネット上で構造をかけるツールが使えるようになります。

↓こんな感じで自由に構造を描けます

SMILESとして取り出すときは、構造全体をドラッグして選択し、右クリック → Copy as→ SMILESとするとSMILESとして構造をコピーできるので、それをエクセルの所定のセルに張り付ければ、構造を取得することができます。

こんな感じで最初のデータセットを準備することができました。次行きましょう。

2. データセットの分布を可視化

簡単に可視化しておきます。まずデータセットを読み込みましょう。

# ライブラリの読み込み
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import IPython
from IPython.display import display

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Descriptors, Draw
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem.rdFingerprintGenerator import GetMorganGenerator
from rdkit.ML.Descriptors import MoleculeDescriptors
from rdkit.Chem.Draw import IPythonConsole

import mordred
from mordred import Calculator, descriptors

from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from scipy.spatial.distance import cdist

それぞれ使っているversionは以下の通りです。

numpy: 1.26.4
pandas: 2.3.0
IPython: 8.18.1
matplotlib: 3.9.4
seaborn: 0.13.2
rdkit: 2025.03.3
mordred: 1.2.0
sckit-learn: 1.6.1
scipy: 1.13.1

リンク先のエクセルを所定のフォルダに保存して、pandasで読み込みます。

# データ読み込み、0列目(Material)をindexに
data = pd.read_csv('material_data.csv', index_col=0)

display(data.head())

Materialという列をindexに指定し、先ほど説明した、SMILES, Type, PLの列が並んでいますね。まずはPL波長の分布をみてみます。

# Typeごとのkde plot
plt.figure(figsize=(14, 7))
sns.kdeplot(data=data,
            x="PL",
            hue="Type",
            fill=True,
            common_norm=False,
            alpha=0.5)

plt.tight_layout()
plt.show()

波長域として紫外(300-400nm)から近赤外(700-900nm)までカバーしていることが分かります。ここで今回使っている骨格を紹介しておきます。

・炭化水素系(アントラセン etc.)：名前の通りで炭素と水素のみからなる材料群です。紫外～青色を示します。量子収率(PLQY)というどれくらい光るかの指標があるのですが、これらの材料はほぼ100%を示すので、効率的に発光します。

・DABNA：京大の畠山先生らによって開発された材料で、現在注目を集めている材料です。実際は特殊な発光過程を経て蛍光を呈す材料となっており、発光スペクトルの幅がかなり小さく、色純度の高いきれいな青色を示します。

・キノキサリン、BTD、キナクリドン：ヘテロ原子(炭素原子以外の原子)が入った芳香族化合物です。置換基によっては緑から赤色まで色をチューニングすることができます。いろんな分野で研究されている非常に有名な材料群です。

・BODIPY：ホウ素原子に4つの結合がくっついている材料群です。ホウ素アニオンとなることで4つの結合をつくることができます。2つはフッ素原子との共有結合、2つは酸素原子 or 窒素原子との配位結合という形式で構成されています。これらも置換基のチューニングにより緑から近赤外領域の蛍光を呈します。

・BBTD：BTDのお兄ちゃん的な存在です。近赤外領域の蛍光を呈します。この材料も精力的に研究されているのですが、なかなか有機溶媒に溶けないのはご愛敬です((+_+))

念のため箱ひげ図でも確認しておきます。

# 箱ひげ
plt.figure(figsize=(14, 7))
sns.boxplot(x='Type',
            y='PL',
            data=data)

plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

これを見ると、中心の骨格でおおよそ発光波長が決まっていることが改めて分かります。

PLの波長は溶媒によって変化するので注意です。今回は代表値を使っており溶媒の情報は載せていませんが、溶媒の極性が変化すると±10nm程度は簡単に動きます。ばらつきがある値という認識でお願いします。もちろん溶媒の情報も説明変数に入れれば、精度は上がると思うのですが、論文によって測定に使う溶媒が違ったりするので、全部の情報を記載するのは限界があります(心の声：みんな測定する溶媒揃えろよ！！)

3. 記述子の計算

長くなりましたが、いよいよ記述子の計算をしていきましょう。データセットを読み込むだけでは、機械学習モデルを構築することができません。PL波長を予測するための説明変数がないからです。説明変数を準備していきましょう。ここでは分子記述子と呼ばれる有機材料特有の説明変数を作っていきます。分子記述子とは、先ほど用意したSMILESをベースに計算できる分子の個性を表す数値のことです。

これには、例えば以下のような様々な情報が含まれます。

サイズや構造に関する情報: 分子量、原子の数、環の数、分岐の多さなど
物理化学的な特性: 分子の極性（電子の偏り）、水への溶けやすさ（logP）、水素結合のしやすさなど
3次元的な形状に関する情報: 分子の表面積や体積など

このように、一つの分子を数百〜数千種類もの多彩な数値（特徴量）で表現することで、機械学習モデルがその分子の「個性」を深く理解できるようになります。

今回使う分子記述子の計算には、RdkitとMordredの二つを使います。

■ RDKit

RDKitは、ケモインフォマティクス（化学情報学）の分野でデファクトスタンダード（事実上の標準）となっている、非常に有名で強力なオープンソースのライブラリです。

単に記述子を計算するだけでなく、分子構造の読み書き、部分構造検索、フィンガープリントの生成など、化学データを扱う上で必要な機能が幅広く揃っています。今回はこのRDKitを使い、分子量やlogP（脂溶性の指標）といった、基本的で広く使われている記述子を計算します。

公式サイト: https://www.rdkit.org/docs/index.html

■ Mordred

Mordredは、阪大の森脇先生らによって開発された分子記述子の計算に特化したPythonライブラリです。その最大の特徴は、計算できる記述子の数が膨大かつ網羅的である点です。

RDKitが基本的な記述子をカバーするのに対し、Mordredは2D記述子・3D記述子を合わせて1800種類以上もの記述子を計算することができます。これにより、より多角的に分子の「個性」を捉え、モデルの予測精度向上に繋がる可能性のある特徴量を見つけ出すことを目指します。

公式サイト (GitHub): https://github.com/mordred-descriptor/mordred

それでは実際にこれらを計算していきましょう。先のデータセットを読み込んだ状態で、下記コードを実行しましょう。

# smiles列を抽出
smiles = data['SMILES']

# 計算する記述子名の取得
descriptor_names = []

for descriptor_information in Descriptors.descList:
    descriptor_names.append(descriptor_information[0])

print('計算する記述子の数：', len(descriptor_names))

計算する記述子の数： 217

RDKitで計算できる記述子(説明変数)の数は217個だそうです。

# 記述子の計算
# 計算機
descriptor_calculator = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)

# 計算した記述子を入れるリスト
descriptors_rdkit = []

print('分子の数：', len(smiles))

for index, smiles_i in enumerate(smiles):
    print(index+1, '/', len(smiles))
    # molに
    mol = Chem.MolFromSmiles(smiles_i)
    # 追加
    descriptors_rdkit.append(descriptor_calculator.CalcDescriptors(mol))

# DataFrameに
descriptor_df = pd.DataFrame(descriptors_rdkit, index=data.index, columns=descriptor_names)

このコードはそのままコピペでOKです。何をしているかというと、記述子の計算機(MolecularDescriptorCalculator)をインスタンスとして呼び出して、その計算機が1分子ずつ記述子を計算します。そして計算結果をdescriptors_rdkitというリストに入れて、最後にDataFrame形式に変換しています。

中身を少し見てみましょう。

display(descriptor_df.head())

見づらくて恐縮ですが、どの分子もちゃんと217個の記述子が計算されてそうです。
そして、後々使うので保存しておきます。

# 保存
descriptor_df.to_csv('descriptor_rdkit.csv')

次はMordredの計算です。まずはMolファイルを作ります。これも分子を表現する形式の一つです。さっきもちゃっかり使っています。基本、記述子の計算をするインスタンスにはMol形式を入力します。

# mol file作成
mols = []
molhs = []

for smiles_i in smiles:
    # 通常のmol file
    mol = Chem.MolFromSmiles(smiles_i)
    mols.append(mol)

    # 3次元構造の入ったmol file
    molh = Chem.AddHs(mol)
    AllChem.EmbedMolecule(molh, AllChem.ETKDG())
    molhs.append(molh)

Mordredは2次元、3次元両方の記述子を計算できるので、Molファイルと水素原子が付加したMolhファイルの両方を作ります。2次元用のMolファイルは簡単で、SMILESをMolに変換するだけです。

3次元は工夫が必要で、水素原子を付加した後、ETKDG(Experimental-Torsion-angle-preference Knowledge-based Distance Geometry)と呼ばれる手法で3次元構造を作ります。名前の通りで、原子同士の距離のルールを決める手法で、知識ベース、つまり化学的に妥当な原子間の結合距離や角度の情報をルールとして使い、実験的に観測されたねじれ角のデータを元に構造を生成させます。簡単に言うと、立体反発が少なくなるように構造を最適化するイメージです。

ETKDGは電子状態を考慮していないので注意が必要です。今回はあくまでスクリーニング用に簡便に計算できる手法として、ETKDGを用いて3次元構造を決めています。正確な3次元構造が欲しいのであれば、DFT計算などの高精度な手法を使うべきです。

それでは生成したMolを計算機にぶち込みましょう。

# 2次元記述子の計算
calc_2d = Calculator(descriptors, ignore_3D=True)
df_2d = calc_2d.pandas(mols)

# 3次元記述子の計算
calc_3d = Calculator(descriptors, ignore_3D=False)
df_3d = calc_3d.pandas(molhs)

print(f"計算直後の2D記述子 shape: {df_2d.shape}")
print(f"計算直後の3D記述子 shape: {df_3d.shape}")

計算直後の2D記述子 shape: (1258, 1613)
計算直後の3D記述子 shape: (1258, 1826)

分子の数にもよりますが、数分ほど計算に時間がかかります。計算結果はDataFrameとして取得できます。ただ、文字列(False or True)、数字なのに数字として認識されていないものがあったりするので、このままではモデルには入力できません。そこで前処理を実施します。

# 2D記述子
# まず、全ての非数値列を数値化しようと試みる (True/Falseは1/0になる)
non_numeric_cols_2d = df_2d.select_dtypes(exclude='number').columns
for col in non_numeric_cols_2d:
    df_2d[col] = pd.to_numeric(df_2d[col], errors='coerce')

# それでも数値化できなかった列（エラーオブジェクトなど）があれば、その列を削除する
cols_to_drop_2d = df_2d.select_dtypes(exclude='number').columns
if not cols_to_drop_2d.empty:
    print(f"完全に数値化できなかった2D列を削除: {list(cols_to_drop_2d)}")
    df_2d.drop(columns=cols_to_drop_2d, inplace=True)

# 3D記述子でも同様の処理
non_numeric_cols_3d = df_3d.select_dtypes(exclude='number').columns
for col in non_numeric_cols_3d:
    df_3d[col] = pd.to_numeric(df_3d[col], errors='coerce')

cols_to_drop_3d = df_3d.select_dtypes(exclude='number').columns
if not cols_to_drop_3d.empty:
    print(f"完全に数値化できなかった3D列を削除: {list(cols_to_drop_3d)}")
    df_3d.drop(columns=cols_to_drop_3d, inplace=True)

print(f"変換後の2D記述子 shape: {df_2d.shape}")
print(f"変換後の3D記述子 shape: {df_3d.shape}")

変換後の2D記述子 shape: (1258, 1611)
変換後の3D記述子 shape: (1258, 1824)

このコードもコピペで構いませんが、一応解説すると、まずintやfloatでないものが入っている列を抽出します。その抽出した列に対して数値化を試みます。これは本当はfloat型やint型のような数字として認識されなければならないものが、間違ってobject型のような文字として認識されている可能性があるからです。そして数値変換の際に、errors='coarse'とすることで文字列を変換できなかったものは無理やりNaNに置き換えてしまいます。念のため、最後に残っている文字列をdropで消去して終わりです。

なぜこういった処理をするかというと、文字列が入っていると予測モデルに入力ができなかったので、最初は数字だけ抜き出すために、

df_2d = df_2d.select_dtypes('number')
df_3d = df_3d.select_dtypes('number')

このコードで実行したのですが、これだと、2次元と3次元で形状の数が同じになってしまい、せっかく3次元構造を作ったのに計算に使えないという問題が発生しました( ;∀;)
そこでチャッピーにコードを修正してもらいました。色々苦戦したわけです。

また、長くなりましたが本題に戻ります。これで終わりではなく、欠損値が多く入っている列は消去します。

# 欠損値の割合が30%以上の列を削除
threshold = 0.3  

# 2D記述子
valid_cols_2d = df_2d.columns[df_2d.isnull().mean() < threshold]
df_2d_filtered = df_2d[valid_cols_2d]

# 3D記述子
valid_cols_3d = df_3d.columns[df_3d.isnull().mean() < threshold]
df_3d_filtered = df_3d[valid_cols_3d]

print(f"NaNが{threshold*100}%未満の2D記述子の数: {df_2d_filtered.shape[1]}")
print(f"NaNが{threshold*100}%未満の3D記述子の数: {df_3d_filtered.shape[1]}")

NaNが30.0%未満の2D記述子の数: 1427
NaNが30.0%未満の3D記述子の数: 1480

結局、あまり変わらんやんけ！！という突っ込みは置いといて、ひとまず後処理終了です。それでもRDKitよりかは説明変数の数は多い結果になりました。予測モデルに入れる前には別の後処理をするので、さらに減るのですが。。。。

こちらも保存しておきます。

# 保存
df_2d_filtered.index = data.index
df_2d_filtered.to_csv('descriptor_mordred_2d.csv')

# 保存
df_3d_filtered.index = data.index
df_3d_filtered.to_csv('descriptor_mordred_3d.csv')

ひとまずこれで記述子の計算は終わりです。一回休憩しましょう。コーヒーでも飲んできて下さい

☺️☕️ 　　　　　🐵
😦☕️🐵
😇 　　　　　🐵☕️

4. 主成分解析(PCA)、t-SNE

あれ！？もう戻ってきたのですか、やる気ありますね。ではそのやる気を削がないように次進みましょう。このセッションで取り扱う内容は少し難しいので、完全に理解できなくても大丈夫です。簡単に言うと、予測モデルを作る前に、説明変数と目的変数との関係を少しカンニングします。

先ほど保存した記述子を改めて読み込んで、最初のデータセットとくっつけましょう。全部のパターンを記載すると量が多くなってしまうので、descriptor_typeという変数に使いたい記述子の名前を記入して、一つ選んでください。今回はmordred_3dで行きます。

# 入力：読み込みたい記述子のタイプを選択
descriptor_type = 'mordred_3d'  # 'rdkit' or 'mordred_2d' or 'mordred_3d'

# ベースのデータ
dataset = pd.read_csv('material_data.csv', index_col=0)

# 条件に応じて記述子ファイルを読み込む
if descriptor_type == 'rdkit':
    des = pd.read_csv('descriptor_rdkit.csv', index_col=0)
elif descriptor_type == 'mordred_2d':
    des = pd.read_csv('descriptor_mordred_2d.csv', index_col=0)
elif descriptor_type == 'mordred_3d':
    des = pd.read_csv('descriptor_mordred_3d.csv', index_col=0)
else:
    raise ValueError(f"未知のdescriptor_type: {descriptor_type}")

# 結合
dataset_full = pd.concat([dataset.reset_index(), des.reset_index(drop=True)], axis=1)
dataset_full = dataset_full.set_index('Material')

# 確認
print(dataset.shape, des.shape, dataset_full.shape)

(1258, 3) (1258, 1480) (1258, 1483)

ちゃんと元のデータセットと記述のデータセットがくっついていることが分かります。indexが重複しているので、pd.concatの際に片方を消すことを忘れずに。

このデータセットには学習データとテストデータが混ざっていたので、まず学習データのみを抽出します。dropnaの引数でsubset='PL'とすると、PLが空欄のデータ(テストデータ)を削除できるので、結果として学習データだけになります。SMILESとTypeの列も不要なのでここで削除しておきます。

# PLに欠損値が入っている行を消して、学習データのみにする
dataset_train = dataset_full.dropna(subset='PL')
# SMILESとTypeも消しておく
dataset_train = dataset_train.drop(['SMILES', 'Type'], axis=1)

print(dataset_train.shape)

(251, 1480)

学習データは251個となりました。次に先ほど述べた追加の前処理をします。

# infをNaNに置き換え
dataset_train = dataset_train.replace(np.inf, np.nan).fillna(np.nan)
dataset_train = dataset_train.drop(dataset_train.columns[dataset_train.isnull().any()], axis=1)

# 学習データのstdが0の列を特定
zero_std_cols = dataset_train.columns[dataset_train.std() == 0]

# 学習データから同じ列を削除
dataset_train = dataset_train.drop(columns=zero_std_cols)

print(dataset_train.shape)

(251, 1220)

さらに減って説明変数の数は1220個となりました。この処理は何をしているかというと、無限大(inf)となっている値が含まれていると予測モデルに入れることができないので、NaNに変換して消しています。また、標準偏差が0の変数(全てのデータで値が同じ説明変数)は標準化ができないので、これも消しておきます。

まずは簡単に計算できる相関係数を出します。説明変数と目的変数の関係が線形の関係になっていれば、1に近づく指標です。

# 相関係数の計算
corr_coef = dataset_train.corr()

# 相関係数の絶対値（PL自身との相関を除く）
pl_corr_abs = corr_coef['PL'].drop('PL').abs()

# ヒストグラム描画
sns.histplot(pl_corr_abs, bins=50)

plt.xlabel('Abs(Correlation with PL)')
plt.ylabel('Number of Features')
plt.title('Histogram of |Correlation| with PL')
plt.axvline(x=0.7, color='red', linestyle='--', label='Threshold = 0.7')
plt.legend()
plt.tight_layout()

plt.show()

相関係数のヒストグラムも作りました。相関係数が0.7のところに点線をいれています。0.7以上で強い相関があると言われています。中身が何か確認してみます。

# PLとの相関が±0.7以上の説明変数を抽出
abs_07 = pl_corr_abs[pl_corr_abs.abs() >= 0.7].sort_values(ascending=False)

# 結果を表示
print(abs_07)

MID_h      0.767951
nHetero    0.759882

2つの記述子が0.7以上の相関係数を示しています。記述子の内容に関してはここでは深く議論しません。もちろん、一つ一つの記述子の中身を理解しようとすることは大事なのですが、これは後程、予測モデル構築後にやります。

■ PCA

また長くなりましたが、いよいよ主成分解析（PCA: Principal Component Analysis）です。簡単に説明すると、

「高次元のデータを、情報をなるべく保ちながら少ない次元に圧縮する手法」です。

たくさんの説明変数を情報(分散)を多く含む軸に変換し、次元数を減らして、可視化・ノイズ除去・学習の効率化に役立てます。それでは、コード行きます。

# 目的変数
y = dataset_train['PL']
# 説明変数（標準偏差が０のものを除いたもの）
X = dataset_train.drop('PL', axis=1)

# 標準化
autoscaled_X = (X - X.mean()) / X.std()

# PCA
pca = PCA()
pca.fit(autoscaled_X)

# loading、各説明変数がその主成分に対してどれだけ寄与しているか
loading = pd.DataFrame(pca.components_.T, index=X.columns, columns=[f'PC{i+1}' for i in range(pca.n_components_)])

# score、変換後の座標
score = pd.DataFrame(pca.transform(autoscaled_X), index=X.index, columns=[f'PC{i+1}' for i in range(pca.n_components_)])

これでOKです。loadingとscoreという新しい変数がでてきますが、これも簡単に言うと、

loading : 主成分がどの変数をどれだけ使ってできたか（荷重）
score : 各データが主成分軸上のどこにあるか（新しい座標）

ということです。このscoreを使って、新しくできた第一主成分と第二主成分をプロットし、目的変数の値で色付けもします。

# 第一主成分と第二主成分の散布図
plt.scatter(score.iloc[:, 0], score.iloc[:, 1], c=y, cmap=plt.get_cmap('jet'))
plt.colorbar()

plt.xlabel('t_1(PCA)')
plt.ylabel('t_2(PCA)')

plt.show()

散布図から、目的変数が新しい特徴空間で連続的に変化していることが分かり、主成分スコアを使うことで目的変数をある程度予測できることが期待できます。この後出てくるPLSモデルでの精度が期待できるということです。

■ t-SNE

t-SNE（t-distributed Stochastic Neighbor Embedding)も実行してきましょう。こやつはかなり難しく、実装はできるのですが理論は理解できていません。簡単に説明すると(さっきから簡単にしか説明してないような。。。)

「似たものは近く、違うものは遠くになるように、データを低次元（2次元など）に配置する手法」です。

チャッピーに要約してもらいました(^^♪
困ったら生成AIに頼りましょう笑

クラスタ(データが集まっている塊)になっているかを視覚的にとらえることができるので、よく使われる手法です。実装行きます。まずは金子先生の書籍を参考に、下記関数を定義します。

def k3n_error(x_1, x_2, k):
    """
    k-nearest neighbor normalized error (k3n-error)

    可視化する前のサンプル間のユークリッド距離と、可視化後のサンプル間のユークリッド距離に基づく方法

    サンプルの近接関係は多次元空間⇔2次元平面で保持されているか？ということで、細かく分けると、

    1. 多次元空間において近いサンプル同士は、2次元平面においても近いか？
    2. 2次元平面において近いサンプル同士は、多次元空間においても近いか？

    k3n_errorは1と2を同時に評価する、2つとも近ければよい。
    サンプル同士が近いほど、k3n_errorは小さくなる

    k3n-error = k3n-Z-error + k3n-X-error

    Parameters
    ----------
    x_1: numpy.array or pandas.DataFrame
    x_2: numpy.array or pandas.DataFrame
    k: int
        The numbers of neighbor

    Returns
    -------
    k3n_error : float
        k3n-Z-error or k3n-X-error
    """
    x_1 = np.array(x_1)
    x_2 = np.array(x_2)

    x_1_distance = cdist(x_1, x_1)
    x_1_sorted_indexes = np.argsort(x_1_distance, axis=1)
    x_2_distance = cdist(x_2, x_2)

    for i in range(x_2.shape[0]):
        _replace_zero_with_the_smallest_positive_values(x_2_distance[i, :])

    identity_matrix = np.eye(len(x_1_distance), dtype=bool)
    knn_distance_in_x_1 = np.sort(x_2_distance[:, x_1_sorted_indexes[:, 1:k + 1]][identity_matrix])
    knn_distance_in_x_2 = np.sort(x_2_distance)[:, 1:k + 1]

    sum_k3n_error = (
            (knn_distance_in_x_1 - knn_distance_in_x_2) / knn_distance_in_x_2
    ).sum()
    return sum_k3n_error / x_1.shape[0] / k

def _replace_zero_with_the_smallest_positive_values(arr):
    """
    Replace zeros in array with the smallest positive values.

    Parameters
    ----------
    arr: numpy.array
    """
    arr[arr == 0] = np.min(arr[arr != 0])

長々と書いていますが、コードを覚える必要はなくコピペでいいです。何をしているかというと、t-SNEのハイパーパラメータにperplexityと呼ばれるものがあるのですが、これを最適化するための指標として、金子先生はk3n-errorというものを使っています。「高次元でのk近傍の関係」と「低次元でのk近傍の関係」が似ているかどうかをユークリッド距離ベースで評価しています。

さっそくこの関数を使って、最適なperplexityを探してみましょう。

# k(近傍数)の数をサンプル数に応じて決定
n_samples = autoscaled_X.shape[0]
k_in_k3n_error = min(30, max(5, int(np.sqrt(n_samples))))  # 5〜30の範囲に制限

# perplexityの候補
candidates_of_perplexity = np.arange(5, 105, 5, dtype=int)

# k3n-errorを用いたperplexityの最適化
k3n_errors = []

for index, perplexity in enumerate(candidates_of_perplexity):
    # print(index + 1, '/', len(candidates_of_perplexity))

    # tsneの実行
    t = TSNE(perplexity=perplexity, n_components=2, init='pca', random_state=10).fit_transform(autoscaled_X)
    # tを標準化している
    scaled_t = (t-t.mean()) / t.std(axis=0, ddof=1)

    k3n_errors.append(k3n_error(autoscaled_X, scaled_t, k_in_k3n_error) + \
                      k3n_error(scaled_t, autoscaled_X, k_in_k3n_error))

# 最適な値を抽出
optimal_perplexity = candidates_of_perplexity[np.where(k3n_errors == np.min(k3n_errors))[0][0]]

print(f"使用するk（近傍数）: {k_in_k3n_error}")
print(f"最適なperplexity: {optimal_perplexity}")

使用するk（近傍数）: 15
最適なperplexity: 85

# 可視化
plt.scatter(candidates_of_perplexity, k3n_errors, c='blue')

# 右上に注釈を追加（位置は yのmax にちょっと足したくらい）
text_str = f'k = {k_in_k3n_error}\noptimal perplexity = {optimal_perplexity}'
plt.text(0.95, 0.95, text_str,
         transform=plt.gca().transAxes,
         fontsize=10,
         verticalalignment='top',
         horizontalalignment='right',
         bbox=dict(boxstyle='round,pad=0.3', edgecolor='black', facecolor='white'))

plt.xlabel('perplexity')
plt.ylabel('k3n-errors')
plt.show()

kの値はサンプル数に応じて変更するようにしています。小さすぎるとノイズに敏感になってしまい、大きすぎると局所性が評価できないためです。そして、k3n-errorの評価指標により最適なperplexityが得られました。この値を使って、t-SNE本番を実行しましょう！！

# 本番t-SNE
t = TSNE(perplexity=optimal_perplexity, n_components=2,
         init='pca', random_state=10).fit_transform(autoscaled_X)

t = pd.DataFrame(t, index=X.index, columns=['t_1', 't_2'])

# t1とt2の散布図
plt.rcParams['font.size'] = 18
plt.scatter(t.iloc[:, 0], t.iloc[:, 1], c=y, cmap=plt.get_cmap('jet'))
plt.colorbar()

plt.xlabel('t_1(TSNE)')
plt.ylabel('t_2(TSNE)')
plt.show()

t-SNEで得られたt-SNEスコアもPCAでの主成分スコアのように可視化で使うことができます。散布図からサンプルが集まっているクラスタのようなもの見えました。また、目的変数の値も近いサンプル同士が近接していることから、似たような特徴を持つサンプルは似たようなPL波長を示すことが分かります。

まとめ

これにて第1回終了です。いかかでしたか？内容はボリューミーですが、まずは1つずつ実行してどういった流れが掴んでいただけたらと思います。気になる点あればコメントお願いします🙇‍♂️

最後まで読んでくださり有難うございました。次回はいよいよ予測モデルが色々出てきます!(^^)!

引用文献

・化学のためのPythonによるデータ解析・機械学習入門
https://datachemeng.com/post-4014/

本記事(次回以降の記事も)の流れは、金子先生の書籍をベースに記載しています。何も勉強していない状態でこの本を読むと難しく感じて、くじけてしまいそうになるのですが、ある程度のpython知識を身に着けた後に見返すと、初歩的な内容から発展的な内容まで詳しく記載されていることに気づき、良書であることが分かります。金子先生、いつもお世話になっております。。。ケモインフォマティクスをやっている人なら絶対みんな読んでいると勝手に思っています。

・化学の新しいカタチ
https://future-chem.com

こちらのサイトもかなり参考になります。分子構造の可視化などいつも参考にさせていただいています。psi4の使い方などもレクチャーされており、大変ありがたいです。

おまけ

ETKDGで作った3次元構造はpythonで可視化できます。先ほど作ったmolhsから構造を呼び出せば見れます。

# 描写
import random

# ランダムなインデックスを取得（0 ～ len(molhs)-1）
idx = random.randint(0, len(molhs) - 1)

# その分子を3Dで描画
IPythonConsole.drawMol3D(molhs[idx])

↓こんな感じ、かっこいいね！笑