化学⼯学特論２ 第１４回 2019年1月21日 (月) 0 理⼯学部 応用化学科 データ化学⼯学研究室 専任講師 ⾦⼦ 弘昌

これまでを少し思い出しましょう 第1,2,3,4回 Pythonの基礎・アルゴリズム • フィボナッチ数列、ソート、アレニウスの式など 第5回 化学⼯学計算 • 連続槽型反応器 一次反応・逐次反応、律速 第6回 データ分布、検定 • 少ないデータから、真の平均の範囲を推定する 第7回 相関係数 • 少ないデータから、真の相関係数の範囲を推定する 第8回 実験計画法 • 多くの実験候補の中から、どうやって少数の候補を選ぶか︖ 第9回 化学構造の扱い • 化学構造の数値表現、分⼦量の計算とか 1

これまでを少し思い出しましょう 第10回 データの可視化 • 多くの変数で表現されたデータを、平⾯上に⾒えるようにする ⁃ 主成分分析、tSNE 第11回 クラスタリング、クラス分類 • クラスタリング: サンプルを塊ごとに自動的に分ける ⁃ 階層的クラスタリング • クラス分類: サンプルをクラスの情報(正解)に基づいて分類する ⁃ 決定木、ランダムフォレスト 第12回 回帰分析 • 説明変数 X と目的変数 y との間の定量的関係を求める ⁃ PLS, LASSO 2

データ解析がどのように応用されるのか︖ 分⼦設計 材料設計 プロセス管理・スペクトル解析 • Chemoinformatics (ケモインフォマティクス) Chemistry + Informatics • Materialsinformatics (マテリアルズインフォマティクス) Materials + Informatics • Processinformatics (プロセスインフォマティクス) Process + Informatics 3

回帰モデル・クラス分類モデルの応用例 4 X: 説明変数 y: 目的変数 回帰モデル クラス分類モデル y = f( X ) モデリング データベース yの推定値 新しいデータ xnew 構造記述⼦*など 予測 物性・活性など x1 x2 1 2 2 1 3 3 データ1 データ2 データ3 例) X: 2変数 データ数: 3 線形モデル *化学構造の情報を数値化したもの 例) 分⼦量、炭素原⼦の数、 ベンゼン環の数 y 活性なし 活性なし 活性あり 1 2 3 1 2 3 クラス分類 モデル

回帰モデル・クラス分類モデルの応用例 5 X: 説明変数 y: 目的変数 回帰モデル クラス分類モデル y = f( X ) モデリング データベース yの推定値 新しいデータ xnew 構造記述⼦*など 予測 物性・活性など x1 x2 1 2 2 1 3 3 データ1 y 5.1 3.9 9.2 データ2 データ3 例) X: 2変数 データ数: 3 線形モデル y = x1 + 2x2 + 誤差 回帰モデル *化学構造の情報を数値化したもの 例) 分⼦量、炭素原⼦の数、 ベンゼン環の数

モデルの構築と予測 構築すべきモデル • フィッティング精度 (モデル構築用データに対する精度) が高い ⁃ データにはノイズが内在するため、完璧に適合するモデルは 意味が無い • 予測精度 (新しいデータに対する予測性能) が高い • 分かりやすい 予測する際に望まれること • 精度良く (誤差が小さいように) 予測したい • モデルの信頼性を知りたい (どれくらいの予測誤差になりそうか︖) ⁃ 推定したいデータによって、モデルの信頼性は異なる → モデルの適用範囲 • 迅速に予測したい 6

誤差 7 精度が高い 誤差が小さい 例) 世界で一番精度の高い時計・・・誤差は160億年に±1秒 モデル誤差 測定誤差  測定ごとにばらつく誤差  一定の傾向を持った誤差 例) 温度計への汚れ付着  実測値と回帰モデルによる計算値および推定値との 間の誤差 モデル構築においては、測定誤差とモデル誤差を考慮に入れる必要がある 仕方のない誤差もある (誤差 0 にはできない)

回帰モデルを構築する 8 x y ある x から y を推定する回帰モデルを構築する y = f( x ) の f を決定する データ

どのモデルが良いか︖ 9 x y x y x y A: 線形モデル y = a1 x + a2 B: 非線形モデル1 C: 非線形モデル2 フィッティング誤差の小さい C が良いか︖ モデル誤差 (フィッティング誤差)

別のデータを用いてモデルを検証する 10 x y : モデル構築に用いたデータ (トレーニングデータ) : モデル構築に用いていないデータ (バリデーションデータ) データ

検証した結果 11 x y x y x y モデル誤差 (予測誤差) A: 線形モデル ○ B: 非線形モデル1 ◎ C: 非線形モデル2 × 線形モデルで十分か 予測精度が 高いモデル モデルがトレーニング データに過度に適合 (オーバーフィッティング) 適切な回帰モデルを構築するためには、 フィッティング誤差だけでなく、予測誤差を考慮に入れることが重要

誤差をまとめてどう考えるか︖ 複数の誤差を一つの指標で表して、モデルの良し悪しを検討する 回帰分析 • 決定係数 r2 • 根平均二乗誤差 (Root Mean Squared Error, RMSE) • 平均絶対誤差 (Mean Absolute Error, MAE) など クラス分類 • 混同⾏列 (confusion matrix) を計算したのちの、 正解率、精度、検出率、誤検出率、Kappa係数など 12 回帰モデル・クラス分類モデルを評価・比較するためのモデルの検証 (Model validation) https://datachemeng.com/modelvalidation/

回帰分析 決定係数 r2 目的変数のばらつきの中で、回帰モデルによって説明できた割合 1に近いほど回帰モデルの”性能”が高い • どんな“性能”かは、r2 を計算したデータセット・推定値による 相関係数 r を二乗したものとは異なる 異なるデータセットの間で r2 を比較してはいけない 13 ( ) ( ) 2 ( ) ( ) EST 2 1 2 ( ) A 1 1 n i i i n i i y y r y y = = − = − −   y(i)︓i 番目のサンプルにおける 目的変数の値 yEST (i)︓i 番目のサンプルにおける 目的変数の推定値 yA ︓目的変数の平均値 n︓サンプル数

回帰分析 RMSE 平均的な誤差の大きさ 0 に近いほど回帰モデルの”性能”が高い • どんな“性能”かは、RMSE を計算したデータセット・推定値による 異なるデータセットの間で RMSE を比較してはいけない データセットが同じであれば、r2 が大きいほど RMSE は小さい 外れ値 (異常に誤差が大きいサンプル) があると、その値の影響を 受けやすく、RMSE が大きくなりやすい 14 ( )2 ( ) ( ) EST 1 n i i i y y RMSE n = − = 

回帰分析 MAE 平均的な誤差の大きさ 0 に近いほど回帰モデルの”性能”が高い • どんな“性能”かは、MAE を計算したデータセット・推定値による 異なるデータセットの間で MAE を比較しないほうがよい 外れ値 (異常に誤差が大きいサンプル) の影響を受けにくい 15 ( ) ( ) EST 1 n i i i y y MAE n = − = 

クラス分類 混同⾏列・正解率・精度・検出率 混同⾏列 (confusion matrix) 16 予測されたクラス 1 (Positive, 陽性) -1 (Negative, 陰性) 実際の クラス 1 (Positive, 陽性) True Positive (TP) False Negative (FN) -1 (Negative, 陰性) False Positive (FP) True Negative (TN) 正解率 = TP + TN TP + FN + FP + TN 検出率 = TP TP + FN 精度 = TP TP + FP 誤検出率 = FP FP + TN など

もう一つのモデルの検証方法 トレーニングデータとバリデーションデータに分けるほどサンプルが 多くないとき トレーニングデータとバリデーションデータに分けることで、データに 偏りがでてしまうのが嫌なとき • クロスバリデーション 17

クロスバリデーション 例) 3-fold クロスバリデーション 18 X 比較指標の 計算 変数 サンプル y X1 X3 y1 y3 X2 y2 X1 y1 X2 y2 X3 モデル1 y3p y1 y3 y2 y1p y3p y2p ① X2 y2 X3 y3 X1 モデル2 y1p ② X3 y3 X1 y1 X2 モデル3 y2p ③ ① ③ ②

クロスバリデーションの補足 Leave-one-out クロスバリデーション • サンプルを１つ除いて、残りのサンプルでモデルを構築し、 除いたサンプルを推定する、ということをサンプル数だけ繰り返す • 特にサンプル数が多いときに、すべてのサンプルでモデルを構築し、 すべてのサンプルを推定することと似てしまうため、望ましくない 2-fold, 5-fold, 10-foldが一般的 データ数が多すぎて、計算時間がかかりすぎてしまうときは、 トレーニングデータとバリデーションデータとを分ける方法のがよい • ディープニューラルネットワーク 19

モデル選択が必要な⼿法の例 線形 • 部分最小二乗法 (Partial Least Squares, PLS) • リッジ回帰 (Ridge Regression, RR) • Least Absolute Shrinkage and Selection Operator (LASSO) • Elastic Net 非線形 • Support Vector Machine (SVM) [線形も] • Support Vector Regression (SVR) [線形も] • 決定木 (Decision tree) • ランダムフォレスト (Random Forest, RF) • ディープニューラルネットワーク 20 https://datachemeng.com/forstudentsresearchers/

モデル選択した後は︖ トレーニングデータを用いてモデル構築 バリデーションデータを用いた検証、もしくはクロスバリデーションにより モデル選択 テストデータを用いて選択されたモデルの最終チェック • PLS と SVR とでどちらがよいか︖ • モデルがどのくらいの予測性能を持っているか︖ モデルの適用範囲を考慮して、モデルを用いた推定 21

どんなXの値でもモデルに入⼒してよいのか︖ 22 物性(活性)推定モデル y=f(X) X y 入⼒ 出⼒ 入⼒してはいけない X の値がある︕ モデルの適用範囲・適用領域 (Applicability Domain, AD) [1-3] [1] I. V. Tetko, et al., J. Chem. Inf. Model. 2008, 48, 1733. [2] D. Horvath, G. Marcou, A. Varnek, J. Chem. Inf. Model. 2009, 49, 1762. [3] H. Kaneko, M. Arakawa, K. Funatsu, AIChE J. 2011, 57, 1506.

モデルの適用範囲・適用領域のイメージ 夏は暑い • 夏であっても涼しい日はある • クーラーの効いた部屋にいれば暑くない • 北極や南極などでは夏でも寒い 炭化水素のデータで構築した水溶解度を推定するモデル • アルコールの水溶解度を正しく推定できるか︖ 23

モデルの適用範囲・適用領域のイメージ 24 x y : トレーニングデータ : 推定したいデータ : モデル x と y の真の関係 誤差大 誤差大 誤差大 適用範囲内 適用範囲内

モデルの適用範囲・適用領域 25 QSAR [1-3] QSPR [4-6] 『モデルが十分な性能を発揮できるデータ領域』 を定めよう︕ モデルの適用範囲・適用領域 (Applicability Domain, AD) モデルの適用範囲・適用領域 (Applicability Domain, AD) 予測したいデータによって、モデルの信頼性は異なる [1] R.P. Sheridan, et al., J. Chem. Inf. Comput. Sci., 44, 1912-1928 (2004) [2] I.V. Tetko, et al., Drug Discov. Today, 11, 700-707 (2006) [3] D. Horvath, et al., J. Chem. Inf. Model, 49, 1762-1776 (2009) [4] P. Bruneau, N.R. McElroy, J. Chem. Inf. Model, 46, 1379-1387(2006) [5] A. Schwaighofer, et al., J. Chem. Inf. Model, 47, 407-424 (2007) [6] Alexandre V., Igor B., J. Chem. Inf. Model, 52, 1413-1437 (2012) 適切にモデルの適用範囲を設定し、 推定するときは適用範囲内かどうか判断する必要がある

AD の設定 トレーニングデータの範囲 [1] トレーニングデータの中⼼からの距離 [2,3] データ密度 [4,5] アンサンブル学習 [6] 26 [1] H. Kaneko, et al., Comput. Chem. Eng. 35 (2011) 1135–1142. [2] S. Dimitrov, et al., J. Chem. Inf. Model. 45 (2005) 839–849. [3] I. Sushko, et al., J. Chem. Inf. Model. 50 (2010) 2094–2111. [4] I.I. Baskin, et al., Mol. Inf. 29 (2010) 581–587. [5] H. Kaneko, et al., Chemometr. Intell. Lab. Syst. 58 (2001) 109–130. [6] H. Kaneko, et al., J. Chem. Inf. Model. 54 (2014) 2469-2482.

トレーニングデータの範囲 モデル構築用データにおける各説明変数 X の範囲 変数間の相関が大きい場合は、凸包や主成分分析後のスコアを利用 27 x1 x2 : トレーニングデータ : 予測データ x1 x2

トレーニングデータの中⼼からの距離 モデル構築用データの平均までの距離 変数間の相関が大きい場合は、マハラノビス距離を用いる 28 x1 x2 : トレーニングデータの平均 : トレーニングデータ : 予測データ

データ密度 モデル構築用データが密に存在する領域が適用範囲内 29 x1 x2 : トレーニングデータ : 予測データ

データ密度 k最近傍法 (k-Nearest Neighbor method, k-NN) One-Class Support Vector Machine (OCSVM) データ密度が高い モデルの適用範囲内 30 https://datachemeng.com/knn/ https://datachemeng.com/ocsvm/

k-NN法 31 例) k = 3 x1 x2 データ密度 : 高い データ密度 : 低い 距離 : 小さい 距離 : 大きい

アンサンブル学習 複数個のモデルを構築し、それらによる予測値の分散(ばらつき)を用いて 適用範囲を評価 例えば、 • トレーニングデータのサンプル • 記述⼦ をランダムに選択して複数モデルを構築 32 データセット サブデータセット1 サブデータセット2 サブデータセット3 例) モデルを3つ構築した場合 モデル3 予測データ x モデル2 モデル1 ypred 1 ypred 2 ypred 3 ばらつき 大 なら 適用範囲外

データ密度で AD を設定してみましょう︕ 33

スペクトルデータの特徴 波⻑ (波数) が近いと、吸光度 (強度) の値も似ている ノイズが含まれる 吸光度 (強度) の極大値 (ピーク) 以外のデータも重要 34

スペクトルデータの前処理 平滑化 (スムージング) • スペクトル・時系列データを “均す (ならす)” ことでノイズを低減する • やりすぎて極大値・極小値の情報が消えないように注意する 微分 • スペクトル・時系列データの傾きを計算することで、 ⁃ ベースラインを補正する ⁃ 新しいスペクトル情報を抽出する ⁃ 時間変化を得る • 一次微分、二次微分、三次微分、・・・ • 微分するとノイズが大きくなるので注意する 35

単純移動平均 (スペクトルデータ) ある波⻑ (波数) の前後 n 点での強度 (吸光度) の平均値を、 平滑化後の値にする • 波⻑ごとに計算する • (2n+1) を 窓枠の数 と呼ぶ • 端っこの波⻑については、(2n+1) 点とれないこともある 36 波⻑ (波数) 強度 (吸光度) 平均値 (2n+1) 点︓窓枠

線形加重移動平均 (スペクトルデータ) ある波⻑ (波数) の前後 n 点での強度 (吸光度) について、 対象の波⻑から離れるにつれて、線形に重みが小さくなる加重平均の 値を、平滑化後の値にする • (2n+1) を 窓枠の数 と呼ぶ 37 ( ) ( ) ( ) ( ) 1 1 1 1 S, 2 1 1 2 1 2 1 1 2 1 i n i n i i i i n i n i x x n x nx n x x x x n n n − − + − + + − + + + + − + + − + + + = + + + − + + − + + + ⋯ ⋯ ⋯ ⋯ ある波⻑ i における強度を xi とし、平滑化後の値を xS,i とすると、

指数加重移動平均 (スペクトルデータ) ある波⻑ (波数) の前後 n 点での強度 (吸光度) について、 対象の波⻑から離れるにつれて、指数関数的に重みが小さくなる 加重平均の値を、平滑化後の値にする • 波⻑からある程度離れると、重みはほぼ 0 になるため、 窓枠をある程度大きくしておけば、細かい数字は気にしなくてよい 38 2 2 2 1 1 2 S, 2 2 1 i i i i i i x x x x x x α α α α α α α α − − + + + + + + + + = + + + + + + ⋯ ⋯ ⋯ ⋯ ある波⻑ i における強度を xi とし、平滑化後の値を xS,i とすると、 α を 平滑化係数 とよぶ

微分 隣の波⻑・時刻における値との差分をとることで、一次微分 一次微分の値について、隣の波⻑・時刻における値との差分を とることで、二次微分 ・・・ 39

Savitzky-Golay (SG) 法 [1,2] データの平滑化と微分とを同時に⾏う方法 • 窓枠のデータを多項式で近似して、多項式の計算値を 平滑化後の値とする • 多項式の微分係数を微分後の値とする ⁃ 波⻑や時刻ごとに計算 スペクトル解析の分野における前処理の方法として一般的 時系列データの前処理法としても有効 [3,4] 40 [1] A. Savitzky, M.J.E. Golay, Anal. Chem. 36, 1627-1639, 1964. [2] 吉村 季織, 高柳 正夫, Journal of Computer Chemistry, Japan, 11, 149-158, 2012 [3] H. Kaneko, K. Funatsu, Ind. Eng. Chem. Res., 54, 12630-12638, 2015. [4] H. Kaneko, K. Funatsu, J. Chem. Eng. Jpn., 50, 422-429, 2017

SG法の例 41 1100 1150 1200 1250 1300 0 0.5 1 1.5 2 2.5 3 強度 波長 [nm] 1100 1150 1200 1250 1300 0 0.5 1 1.5 2 2.5 3 吸光度 波長 [nm] 1100 1150 1200 1250 1300 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 吸光度 波長 [nm] 1100 1150 1200 1250 1300 -4 -2 0 2 4 x 10-3 吸光度 波長 [nm] 元のスペクトル SG法後 SG (1次微分) SG (2次微分)

SG法 (スペクトルデータ) 42 t : 波⻑ x: 強度 x = t の多項式 = a2 t2 + a1 t + a0 (例) 窓枠の数  多項式の次数  窓枠の数 を事前に決めなければならない

SG法でスムージングをしてみよう︕ 43

⼿法・ハイパーパラメータ・微分次数はどうする︖ ４つの⼿法とハイパーパラメータの値の候補 • 単純移動平均︓窓枠の数 (5, 11, 21, 31, …, 201) • 線形加重移動平均︓窓枠の数 (5, 11, 21, 31, …, 201) • 指数加重移動平均︓平滑化係数 (0.01, 0.02, …, 1) • SG法︓ 多項式の次数 (1, 2, 3, 4) 窓枠の数 (5, 11, 21, 31, …, 201) 微分次数 (場合によってはその組み合わせ) をどのように決めるか︖ 44 ① モデルの検証により選択する ② ノイズの正規分布性により選択する

① モデルの検証による選択 各⼿法・各ハイパーパラメータの値・各微分係数の値で、 回帰分析・クラス分類のモデルの検証を⾏い、 最も検証結果のよい組み合わせを選択する • たとえば、 ⁃ クロスバリデーション推定値の r2 が最も大きい組み合わせ ⁃ バリデーションデータの r2 が最も大きい組み合わせ • モデルの検証︓http://datachemeng.com/modelvalidation/ 45

① モデルの検証による選択 特徴 メリット • モデルの検証の仕方によっては、推定性能の高いモデルを構築できる ⼿法・ハイパーパラメータの値・微分係数 を選択可能 デメリット • 教師ありデータが必要 • モデリングを何回も⾏わなくてはならない (時間がかかる) 46

② ノイズの正規分布性による選択 平滑化前後の値を引くことで、平滑化によって “均(なら)された” ノイズの値を計算できる ノイズは正規分布であると仮定すると、平滑化によって減少したノイズの 分布も正規分布に従う必要がある コルモゴロフ–スミルノフ検定などの正規分布性の検定により、ノイズが 正規分布に従う⼿法・ハイパーパラメータの組み合わせを選択 選択された⼿法・ハイパーパラメータの組の中で、 標準偏差が最も大きい ( = ノイズが最も減少した) 組を選択 詳しくは下の論⽂を参照のこと 47 H. Kaneko, K. Funatsu, J. Chem. Eng. Jpn., 50, 422-429, 2017

② ノイズの正規分布性による選択 特徴 メリット • 教師データ不要 • モデリング不要 (時間がかからない) デメリット • 微分次数は選択できない • 選択の際、モデルの推定性能は考慮されていない 48