こんにちは、代表のみでっとです。
今回は、個人的に初めて参加し、そして最終的に入賞を果たすことができたデータコンペ「第4回Brain(s)コンテスト」に参加してみての感想や、解法についてまとめていきます！（コンテストの詳細はハイパーリンクからどうぞ。）
結構張り切ってまとめたので、是非最後まで読んでいっていただけると嬉しいです！！
普段は大学の友人と共に、初心者向けにPythonの入門記事や発展記事をほのぼのと書いていますので、よろしければチェックしていっていただけると嬉しいです。
また、もし周りにPython始めたての方がいらっしゃれば、紹介していただけるとすごく喜ぶのでよろしくお願いします！！
（↓この団体を作った経緯をまとめてあります、こちらもよろしければ…!）
また、今後も定期的に記事の配信や、データコンペへの参加記を書いていく予定なので、よかったらtwitterのフォローもお願いします！（宣伝ばっかりですみません💦）
目次
富士フィルム主催Brain(s)コンテスト
そもそもBrain(s)コンテストってどんな大会なの？
Brain(s)コンテストというのは、富士フイルム株式会社が若手のエンジニアを目指している方々に対して、人工知能（AI）、機械学習、データサイエンスにより興味を持ってもらい、将来の優秀な技術者を輩出するための取り組みの一環として行われているコンテストです。
コンテストの内容はすごくざっくり言うと、富士フィルム株式会社が完全オリジナルで作成した問題に対して、どれだけ精度良く予測できるモデルを作れるか！というもので、今回のテーマは「富士フイルム×マテリアルズ・インフォマティクス」でした。
マテリアルズ・インフォマティクスって？
私個人の話になるのですが、大学院で材料研究を行っていてマテリアルズ・インフォマティクスに興味があり今回のコンペに参加させていただいた経緯があるので、ここでは詳しく説明させていただければと思います！
マテリアルズ・インフォマティクスというのは最近AI界隈でホットな研究分野の一つで、新規の物性をもつ材料（ex. 従来よりも水をよく弾く材料（超撥水材料））をこれまで依存していた熟練研究者の勘や経験に頼ることなく、機械学習の技術を用いて従来よりもぐっと短い期間で開発できると期待されていて、様々な企業様で取り組まれています。
また、マテリアルズ・インフォマティクスは従来のように、
- 材料の特性が発現するメカニズムを解明 → 2. 物性を予測 → 3. 新規材料開発
というプロセスではなく、
- 材料データベースを数理モデルに組み込み物性を予測 → 2. メカニズムの解明 → 3. 新規材料開発
というプロセスを経るため、従来の材料開発とは逆の順序であることが特徴の一つです。
「数理モデルは作ったから、あとは任せたぞ～」と数多くのパラメータから特定の物性値をコンピューターに計算させて予測させてるの、なんだかすごく先進的なことをやっている感じがして面白いですよね！！（押し強め）
今回のコンテストでは、このプロセスのうち核となる、
「材料データベースを数理モデルに組み込み物性を予測」
という部分を、材料データベースのみ与えられたところから実際に体験でき、最終的に精度のよいモデルを作成するためにどういうアプローチをしていけばいいのか（たまたま上手くいっただけかもしれませんが、、）を学べたので、大変有意義な経験ができたのかなと思っています。
改めまして、本コンテストを開いてくださいました富士フィルム株式会社様に感謝申し上げます、ありがとうございました！
具体的なコンテスト内容
今回のコンペでは、Q1~Q3までの3つの問題があり、それぞれ以下のような内容でした。
Q1. SMILES記法で表された分子に含まれるCの数を同定するプログラムを作成せよ。
Q2. SMILES記法で表された分子の水への溶解度（Water Solubility）を予測せよ。
Q3. SMILES記法で表された分子のエームズ試験の結果を予測せよ。
最終的にはこれら3つの問題の総合順位で（実質Q2とQ3）最終順位が決まり、参加者の中で優秀な成績を収めて上位に入賞した暁には、富士フイルム株式会社が誇る最高峰のカメラ「X Series FUJIFILM X100F」や「X Series FUJIFILM XF10」などの豪華景品をGET!できるということで、ここもやはり少なからずモチベーションの一部になっていました。
最終結果
結果としては、個人成績としてQ2で4位、Q3で8位という望外な成績で初参加のコンテストを終えられ、実力以上の結果が出せたのかなと思っています。
最高峰のカメラである「X Series FUJIFILM XF10」を賞品としていただけるということで、ご連絡をいただいたときは本当に嬉しかったです。
最近は自粛生活の影響であまり外で写真を撮る機会がないのですが、外に安心して出かけられるようになったらたくさん活躍してもらおうかなと内心ワクワクしています。
しかしながら、最終的な順位をよく見てみるとQ2の1位~3位の方々と比較すると0.4%弱（ほんの僅かな差に思えますが、絶壁の様に大きな壁です）も精度が低く、まだまだアプローチのしがいがあったのではないかなと思っています（特に特徴量エンジニアリング）。
同時に、Q3も8位でまだまだ汎用的なモデルを作成する修行が足りていないなと感じました…。
そのため、上位の方々の解法や、富士フィルム株式会社の現役AIエンジニアの方の解法を表彰式にてお聞きするのを楽しみにしていたのですが、運悪く参加できなかったためとても残念です….。（懇親会楽しみにしています！！）
Slack上でのやり取りで一部の方々が解法をブログにあげていただいているのを知り、少し救われました！ありがとうございます！
（以下のリンクからアクセスできるので、よけろしければこちらも合わせてどうぞ！）
富士フィルムコンペ総合1位解法・アプローチ（化合物の溶解度予測2位、変異性予測2位）
Chem Infomatics系の基礎事項まとめ(富士フィルム Brain(s)コンテスト参加記)
それでは、それぞれの課題に対して、個人的にどのようなアプローチをして解いていったのか、次の章から詳しく解説していこうと思います！
なお、公式サイトを確認したところ、今回のコンペで作成したソースコードは公開しても良いようだったので（データセットはNG）、ソースコードのみ公開させていただきます。
私自身、Python歴が1年6ヶ月程度で機械学習の勉強を本格的にやり始めたのは半年前なので、コードが汚い！といった点が多々あるかもしれません。ごめんなさい。（というかあります、最終的には動けばいいやと思って関数化せずに脳筋プレーでごり押ししました(笑)）
コード書くときはこうしたほうがいいよ！などアドバイス等あれば、いつでも大歓迎なのでよろしくお願いします！
Q1. SMILES記法で表された分子に含まれるCの数を同定するプログラムを作成せよ。
SMILES記法って何？
私自身、学部生の時代に4年間理系として様々なことを学んできましたが、聞いたことがなかったので多分、化学を専門でやられている方でないと知らないのかな（？）と思います。
最初、問1を見たときに頭の中に？？？が浮かんで、一瞬挫折しそうになったのを覚えています(笑)
WikipediaでSMILES記法について検索すると、以下のような記述がありました。
「SMILES記法（スマイルスきほう、英語: simplified molecular input line entry system）とは、分子の化学構造をASCII符号の英数字で文字列化した、構造の曖昧性の無い表記方法である。」引用ソース：Wikipedia
つまり、簡単に言うと分子（H2OとかCH3COOH）を表記する決まりの一つだよ！という風にとらえていただけるといいかなと思っています。
Q1.の解法
つまり今回のタスクの場合、例えばプロパンのSMILES表記 CCC に対しては、3、フランのSMILES表記o1cccc1に対しては、4、と出力できるプログラムを書ければよいということですね。
数行のコードで記述できるので、皆さんもよければ挑戦してみて下さい！
正直、この問は後から考えると、Q2でSMILES記法のCの数が有効な特徴量だからとりあえずやってみなさい。という神からのご達しなのかなと感じました。(笑)
こちらはソースコードはありませんが、Q2でnumber_of_carbons関数として実装しているので、参考にされたい方はどうぞ。
Q2. SMILES記法で表された分子の水への溶解度（Water Solubility）を予測せよ。
さて。ついにここまで来ました。
正直このコンペの中で一番時間を使って、最後までなんとか精度を上げようと頑張りました。感覚的にはQ1 1%、Q2 95%、Q3 4%みたいな比率だったので、最終的に1位になれなかったのは本当に悔しかったです。
こちらの問題は、ある時点から毎日パブリックリーダーボード（通称Public Leaderboard、誰のスコアが高いか上位の順位のみ公開される掲示板）が更新されることになっていたので、自分のスコアが全体から見てどの位置にあるのか把握ができたのでその分熱が入りました。
今回の場合、コンテスト終了の5日前にリーダーボードの更新が止まってしまうということだったので、個人的に
コンテスト終了5日前まで：
単一モデルでの予測、特徴量エンジニアリング、パラメータチューニング
コンテスト終了まで：
アンサンブルによる精度向上
をしようという風に決めて取り組んでいました。最終公開時のリーダーボードでは7位と比較的よい順位につけられていて、最後の追い込みとしてアンサンブルを行って1%弱は上がることを想定していたのですが、現実問題はそんなに甘くなかったです…（反省）
おススメの参考書籍
唐突ですが今回のコンペではQ2とQ3を進めるにあたって、個人的に以下の書籍を大変参考にしました。
データ分析界隈では有名な参考書で（通称Kaggle本と呼ばれています）、もしこれからデータコンペに参加してみたい！データ分析ガチ勢になりたい！という方がいらっしゃれば、是非一度手に取ってみてください。
この本は個人的に知り合いの大手企業の現役AIエンジニアの方におススメされた書籍で、参考書として本棚においていらっしゃるとのことだったので、手元に一冊置いておいて間違いない書籍かな思います！おススメです✨
|
|
![[商品価格に関しましては、リンクが作成された時点と現時点で情報が変更されている場合がございます。]](https://hbb.afl.rakuten.co.jp/hgb/1cdc626a.1e86cc5a.1cdc626b.729e8501/?me_id=1278256&item_id=18610772&pc=https%3A%2F%2Fthumbnail.image.rakuten.co.jp%2F%400_mall%2Frakutenkobo-ebooks%2Fcabinet%2F0127%2F2000007790127.jpg%3F_ex%3D240x240&s=240x240&t=picttext "[商品価格に関しましては、リンクが作成された時点と現時点で情報が変更されている場合がございます。]")
与えられたデータセットとタスク
Q2のデータセットは、
・SMILES記法で記された分子
・マテリアルズ・インフォマティクスや、ケモ・インフォマティクス用の有名なオープンソース（RDKit）を用いて計算された物性値
が分子ごとにまとめられたもので、それぞれの分子に対して計200個以上の物性値が既に計算されていました。
このように多くの物性値がまとめられたデータセットが与えられた中、今回Q2で課せられたタスクは、
「200個以上ある物性値の中の一つである、水に対する溶解度(WaterSolubility)を他の物性値から予測して、実際の値と近いを出力する数理モデルを作成せよ！」
というものでした。（評価指数は決定係数R^2）
実際にマテリアルズ・インフォマティクスを進めていくとき、未知の分子に対してより正確な物性値を予測できると新規物質の予測精度が上がるので、数理モデルの精度が高いほうがよいのは言うまでもないですよね。
つまりQ2では、「マテリアルズ・インフォマティクスという分野において最重要タスクである精度を向上させる。」ことを実際にやってみて、社会的に価値のある数理モデルを作成してみよう！という意図で作成されたのかな。と勝手に想像しながら取り組んでいました。
個人的には、実務でもこのようなプロセスが数理モデル作成が進められているのか、というイメージがつかめたのはこの問で大変参考になった部分の一つですね。
Q2.の解法
特徴量エンジニアリング
まず初めに、特徴量エンジニアリングを行いました。
特徴量エンジニアリングとは既存の特徴量から新しい特徴量を作成することを指し、テーブルデータのコンペでは非常に重要な要素であり、良い特徴量を作成できたかによって順位が決まることがほとんどとなっています。
そのため、今回は以下の特徴量の作成を行いました。（特徴量が大事なのを知っていつつも、モデルを作るのに楽しくなってしまっておろそかにしてしまった感が否めないです。結果的にはもっと検討しておいたほうが良かった。）
- SMILES記法で示された分子のCarbonの数（Q1.で作成したプログラム）
- SMILES記法で示された分子の化学結合の種類、原子数
- SMILES記法で示された有機化合物の化学結合の種類が全体に占める割合
以下、ソースコードになります。
#SMILES記法で示された分子のCarbonの数の追加（関数の定義）
def number_of_carbons(df):
carbon = []
for i in range(len(df)):
construction_list = df.at[i, "SMILES"]
n = construction_list.count("C") + construction_list.count("c") - construction_list.count("l")
carbon.append(n)
df["Number_of_Carbons"] = carbon
#SMILES記法で示された分子の化学結合の種類、原子数を出力する関数の定義
def types_of_chemical_bonds(df):
number_of_atoms = []
number_of_SINGLE_bonds = []
number_of_DOUBLE_bonds = []
number_of_TRIPLE_bonds = []
number_of_AROMATIC_bonds = []
for i in range(len(df)):
#原子数および、化学結合の種類を格納するリストを定義
number_of_atom = 0
chemical_bonds = []
smiles = Chem.MolFromSmiles(df.at[i, "SMILES"])
for atom in smiles.GetAtoms():
number_of_atom += 1
for j in range(number_of_atom):
chemical_bonds.append(str(smiles.GetBonds()[j].GetBondType()))
number_of_SINGLE_bonds.append(chemical_bonds.count("SINGLE"))
number_of_DOUBLE_bonds.append(chemical_bonds.count("DOUBLE"))
number_of_TRIPLE_bonds.append(chemical_bonds.count("TRIPLE"))
number_of_AROMATIC_bonds.append(chemical_bonds.count("AROMATIC"))
number_of_atoms.append(number_of_atom)
#pandasのDataFrame型に変換
df["number_of_atoms"] = number_of_atoms
df["number_of_SINGLE_bonds"] = number_of_SINGLE_bonds
df["number_of_DOUBLE_bonds"] = number_of_DOUBLE_bonds
df["number_of_TRIPLE_bonds"] = number_of_TRIPLE_bonds
df["number_of_AROMATIC_bonds"] = number_of_AROMATIC_bonds
#SMILES記法で示された分子の化学結合の種類が全体に占める割合を出力する関数の定義
def percentage_of_chemical_bonds(df):
percentage_of_SINGLE_bonds = []
percentage_of_DOUBLE_bonds = []
percentage_of_TRIPLE_bonds = []
percentage_of_AROMATIC_bonds = []
for i in range(len(df)):
percentage_of_number_of_SINGLE_bonds = df.at[i,"number_of_SINGLE_bonds"] / df.at[i,"number_of_atoms"]
percentage_of_number_of_DOUBLE_bonds = df.at[i, "number_of_DOUBLE_bonds"] / df.at[i, "number_of_atoms"]
percentage_of_number_of_TRIPLE_bonds = df.at[i, "number_of_TRIPLE_bonds"] / df.at[i, "number_of_atoms"]
percentage_of_number_of_AROMATIC_bonds = df.at[i, "number_of_AROMATIC_bonds"] / df.at[i, "number_of_atoms"]
percentage_of_SINGLE_bonds.append(percentage_of_number_of_SINGLE_bonds)
percentage_of_DOUBLE_bonds.append(percentage_of_number_of_DOUBLE_bonds)
percentage_of_TRIPLE_bonds.append(percentage_of_number_of_TRIPLE_bonds)
percentage_of_AROMATIC_bonds.append(percentage_of_number_of_AROMATIC_bonds)
df["percentage_of_SINGLE_bonds"] = percentage_of_SINGLE_bonds
df["percentage_of_DOUBLE_bonds"] = percentage_of_DOUBLE_bonds
df["percentage_of_TRIPLE_bonds"] = percentage_of_TRIPLE_bonds
df["percentage_of_AROMATIC_bonds"] = percentage_of_AROMATIC_bonds/補足：
この特徴量についてですが、化学の専門書で調べたところ、溶かしたい有機化合物と溶液の比誘電率εが近いと溶けやすいとの記述があったため（水の比誘電率は78）、双極子モーメント、DFT計算、Onsager’s equationとかで何とか導出できないかなと思い丸ごと2日費やしたのですが結局できなかったのが心残りです。
ただし、物質の極性基の周りがかさ高過ぎたりすると溶媒和できなくなり水に溶けなくなるので、分子の立体構造を考える必要があり、より精度よく予測するのであれば原子位置の情報を追加するなどできたのかなー、コンペが終わってからぼんやりと考えていました。
数理モデルの作成（コンテスト終了5日前まで）
特徴量の作成後、早速モデルの作成に取り掛かりました。
まず初めに手を付けたのは、とりあえず初手。とKaggle本で述べられているくらいデータコンペではお馴染みの勾配ブースティング木（GBDT：Gradient Boosting Decision Tree）と呼ばれる手法でした。
特徴としては、使いやすく、ある程度の精度がパラメータチューニング無しでも出る点で、このモデル単体で上位に入賞される方もよくいるそうで、王道っていうイメージがぴったりなモデルです。
（ここでは勾配ブースティングの詳しい説明は省きますが、詳細に知りたい方はKaggle本で特徴の説明から、パラメータチューニングの方法、ソースコードの例が分かりやすくまとめられているのでそちらをどうぞ！）
個人的に数理モデルを一から作り上げるのは初めての経験だったので、本に載っているソースコードを1行ずつ追いながら内容の理解をし、同時にクロスバリデーションの実装を並行して進めつつ最初のモデルを作成し終え、遂に最初のモデルを提出して、結果をワクワクして待っていました!!
最初のスコアは0.78弱（確かこのくらい）で、特にパラメータをいじっていないのですが既にある程度精度の高いモデルが得られ、さすが勾配ブースティング木やなぁ～って感じました(笑)
その後、ベイズ最適化や、Kaggle本の方針に従いながらあれこれパラメータチューニングを行い精度の向上を目指しました。徐々に精度が向上していき0.81程度のスコアが安定して取れるようになってきたのですが、ここで詰まってしまい何をすればよいかにとん挫してしまったため、なんとなくクロスバリデーションを行う際の乱数値（randomstate）をいじったところ、クロスバリデーションによるスコアがなんと0.03弱上がる乱数値を発見し、まさかと思って作成したモデルを提出してみたところ、0.82の壁を遂に超えることが出来ました。
その後、1~10000までの乱数値それぞれでのクロスバリデーションスコアを求め、スコアが高かったものを片っ端から提出して精度を検証していく、という脳筋モードに突入しました。(笑)
結果、randomstateが687のときに提出したモデルが一番精度が高く、0.8366という高スコアをたたき出すことが出来、8/2時点で暫定4位になり遂にリーダーボードに名を残すことが出来ました。
そしてコンペ終了5日前まで、他の乱数値を試してみたりパラメータをさらに細かくチューニングしてみたのですがなかなか精度が上がらず、ディープラーニング等の様々な手法を用いて予測精度の向上を試みましたが、暫定7位、最終手段として残していたアンサンブルをここで解放することにしました。
以下、試行錯誤し完成した数理モデル作成用のコードになります。（意外と短い）
import xgboost as xgb
from sklearn.metrics import r2_score as R2
from sklearn.model_selection import KFold
import pickle
import pandas as pd
#先程定義した関数（ここでは省略）
#データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=687)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
# 特徴量と目的変数をxgboostのデータ構造に変換
dtrain = xgb.DMatrix(tr_x, label=tr_y)
dvalid = xgb.DMatrix(va_x, label=va_y)
# xgbのハイパーパラメータの設定
#各値はベイズ最適化でパラメータ探索を行った
params = {'objective': 'reg:squarederror',
'eta': 0.1,
'gamma': 0.45,
'alpha': 0,
'lambda': 1.0,
'min_child_weight': 1.1,
'max_depth': 5,
'subsample': 0.8,
'colsample_bytree': 0.8,
'random_state': 77}
num_round = 900
# モデルの学習
# バリデーションデータもモデルに渡し、学習の進行とともにスコアがどう変わるかモニタリングを行う
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]
model = xgb.train(params, dtrain, num_round, evals=watchlist)
# バリデーションデータでのスコアの確認
va_pred = model.predict(dvalid)
score = R2(va_y, va_pred)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {score:.4f}')
##モデルの保存
pickle.dump(model, open('mymodel.pkl', 'wb'))/数理モデルの作成（コンテスト終了まで）
先に述べたように、コンテスト終了まではアンサンブルを行い精度の向上を目指しました。
アンサンブルって？
アンサンブルとは、複数のモデルを組み合わせてモデルを作ることや、予測を行うことを指し、KaggleやSignateなどのコンペで上位1%に入るモデルでは、ほとんど行われているようです。
中には数100以上の独立したモデルを組み合わせて入賞をされた方がいらっしゃったりと、実際の実務では予算、時間の制約で出来ないことが多いですが精度を突き詰めるという上では有効な手法となります。
これは極端な例になりますが、例えば、ある分子の水への溶解度が0.5だったとします。このとき、モデル1では0.4、モデル2では0.6と予測していた場合、単独のモデルでは0.1の誤差が生じてしまいますが、2つのモデルの平均値をとると0.5となり正確に値を予測できるようになりました！
そのためアンサンブルでは、予測値の表現力を高めるために多様性に富むモデルを組み合わせると有効というのが通説で、
- 勾配ブースティング木
- ニューラルネットワーク
- 線形モデル
- ランダムフォレスト
など、予測値の境界が異なるモデルを使うとより正確に予測できる可能性が高まります。今回のコンペでは提出できるファイルサイズに10MBという制限があったため組み合わせられるモデルの数に限界はあったものの、出来る範囲で試行錯誤して進めていきました。
初手では、最もシンプルなアンサンブル手法で、複数モデルの予測値の平均を出力としてみました。精度の向上を期待したのですが、逆に精度が落ちる結果となり、自分の作成したモデルの多様性が足りないことが判明したので線形回帰や、ランダムフォレストなど他手法の勉強を進めモデルの作成を行いました。
様々なモデルの作成後、再度平均値をとってみましたが精度向上しなかったため、それぞれのモデルの予測値を特徴量とする線形モデルを第2層目のモデルとして構築してみました。
このモデルを提出したところ、0.836 → 0.838と0.02も精度が向上しました。
その後、勉強を進めていったところ、線形モデルにおいて非線形性、相互作用を表現するためには、こちらで特徴量を作成し組み込む必要があることが分かったので、以下のパラメータを特徴量として作成しました。
それぞれのモデルの予測値の2乗、平均、加重平均、分散、標準偏差等
最終的に、これらの特徴量に加えて、アンサンブルするモデルの作成、選定を繰り返し行い0.84259というスコアが得られ、最終4位にランクインすることが出来ました！
以下、作成したコードになります。（ライブラリのインポート部分、特徴量追加部分は割愛します。）
#モデル1
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=687)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
# 特徴量と目的変数をxgboostのデータ構造に変換
dtrain = xgb.DMatrix(tr_x, label=tr_y)
dvalid = xgb.DMatrix(va_x, label=va_y)
# xgbのハイパーパラメータの設定
#各値はベイズ最適化でパラメータ探索を行った
params = {'objective': 'reg:squarederror',
'eta': 0.1,
'gamma': 0.45,
'alpha': 0,
'lambda': 1.0,
'min_child_weight': 1.1,
'max_depth': 5,
'subsample': 0.8,
'colsample_bytree': 0.8,
'random_state': 77}
num_round = 900
# モデルの学習
# バリデーションデータもモデルに渡し、学習の進行とともにスコアがどう変わるかモニタリングを行う
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]
model = xgb.train(params, dtrain, num_round, evals=watchlist)
# バリデーションデータでのスコアの確認
va_pred = model.predict(dvalid)
score = R2(va_y, va_pred)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {score:.4f}')
##モデルの保存
pickle.dump(model, open("mymodel.pkl", "wb"))#モデル2
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=2195)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
# 特徴量と目的変数をxgboostのデータ構造に変換
dtrain = xgb.DMatrix(tr_x, label=tr_y)
dvalid = xgb.DMatrix(va_x, label=va_y)
# xgbのハイパーパラメータの設定
#各値はベイズ最適化でパラメータ探索を行った
params = {'objective': 'reg:squarederror',
'eta': 0.1,
'gamma': 0.45,
'alpha': 0,
'lambda': 1.0,
'min_child_weight': 1.1,
'max_depth': 5,
'subsample': 0.8,
'colsample_bytree': 0.8,
'random_state': 77}
num_round = 900
# モデルの学習
# バリデーションデータもモデルに渡し、学習の進行とともにスコアがどう変わるかモニタリングを行う
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]
model = xgb.train(params, dtrain, num_round, evals=watchlist)
# バリデーションデータでのスコアの確認
va_pred = model.predict(dvalid)
score = R2(va_y, va_pred)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {score:.4f}')
##モデルの保存
pickle.dump(model, open("mymodel_xgb_2.pkl", "wb"))#モデル3
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=687)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
# 特徴量と目的変数をxgboostのデータ構造に変換
dtrain = xgb.DMatrix(tr_x, label=tr_y)
dvalid = xgb.DMatrix(va_x, label=va_y)
# xgbのハイパーパラメータの設定
#各値はベイズ最適化でパラメータ探索を行った
params = {'objective': 'reg:squarederror',
'booster': 'dart',
'eta': 0.1,
'gamma': 0.45,
'alpha': 0,
'lambda': 1.0,
'min_child_weight': 1.1,
'max_depth': 5,
'subsample': 0.8,
'colsample_bytree': 0.8,
'random_state': 77}
num_round = 1000
# モデルの学習
# バリデーションデータもモデルに渡し、学習の進行とともにスコアがどう変わるかモニタリングを行う
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]
model = xgb.train(params, dtrain, num_round, evals=watchlist)
# バリデーションデータでのスコアの確認
va_pred = model.predict(dvalid)
score = R2(va_y, va_pred)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {score:.4f}')
##モデルの保存
pickle.dump(model, open("mymodel_xgb_3.pkl", "wb"))#モデル4
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=687)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
# 特徴量と目的変数をxgboostのデータ構造に変換
dtrain = xgb.DMatrix(tr_x, label=tr_y)
dvalid = xgb.DMatrix(va_x, label=va_y)
# xgbのハイパーパラメータの設定
#各値はベイズ最適化でパラメータ探索を行った
params = {'objective': 'reg:squarederror',
#'learning_rate': 0.1,
'eta': 0.1,
'gamma': 0.45,
'alpha': 0,
'lambda': 1.0,
'min_child_weight': 1.1,
'max_depth': 5,
'subsample': 0.8,
'colsample_bytree': 0.8,
#'base_score': 1.627179,
'scale_pos_weight': 0,
#'random_state': 77,
'seed': 0}
num_round = 920
# モデルの学習
# バリデーションデータもモデルに渡し、学習の進行とともにスコアがどう変わるかモニタリングを行う
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]
model = xgb.train(params, dtrain, num_round, evals=watchlist)
# バリデーションデータでのスコアの確認
va_pred = model.predict(dvalid)
score = R2(va_y, va_pred)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {score:.4f}')
##モデルの保存
pickle.dump(model, open("mymodel_xgb_4.pkl", "wb"))#モデル5
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=27)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
# 特徴量と目的変数をxgboostのデータ構造に変換
dtrain = xgb.DMatrix(tr_x, label=tr_y)
dvalid = xgb.DMatrix(va_x, label=va_y)
# xgbのハイパーパラメータの設定
#各値はベイズ最適化でパラメータ探索を行った
params = {'objective': 'reg:squarederror',
'eta': 0.1,
'gamma': 0.4,
'alpha': 0,
'lambda': 1.0,
'min_child_weight': 4,
'max_depth': 5,
'subsample': 0.9,
'colsample_bytree': 1.0,
'colsample_bylevel':0.1,
'random_state': 77}
num_round = 1000
# モデルの学習
# バリデーションデータもモデルに渡し、学習の進行とともにスコアがどう変わるかモニタリングを行う
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]
model = xgb.train(params, dtrain, num_round, evals=watchlist, early_stopping_rounds=100)
# バリデーションデータでのスコアの確認
va_pred = model.predict(dvalid, ntree_limit=model.best_ntree_limit)
score = R2(va_y, va_pred)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {score:.4f}')
##モデルの保存
pickle.dump(model, open("mymodel.pkl", "wb"))#モデル6
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=687)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
#値の標準化
scaler = StandardScaler()
tr_x = scaler.fit_transform(tr_x)
va_x = scaler.transform(va_x)
#線形回帰で予測
model_1 = LinearRegression()
model_1.fit(tr_x, tr_y)
model_1_pred = model_1.predict(va_x)
model_1_score = R2(va_y, model_1_pred)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {model_1_score:.4f}')
#モデルの保存
pickle.dump(model_1, open("mymodel_LR.pkl", "wb"))#モデル7
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=687)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
#ランダムフォレストで予測
model_2 = RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=20,
max_features='auto', max_leaf_nodes=None,
min_impurity_decrease=0.0, min_impurity_split=None,
min_samples_leaf=1, min_samples_split=3,
min_weight_fraction_leaf=0.0, n_estimators=300, n_jobs=1,
oob_score=False, random_state=2525, verbose=0,
warm_start=False)
model_2.fit(tr_x, tr_y)
model_2_pred = model_2.predict(va_x)
model_2_score = R2(va_y, model_2_pred)
#グリッドサーチで適切なパラメータを算出
"""
from sklearn.model_selection import GridSearchCV
search_params = {
'n_estimators': [5, 10, 20, 30, 50, 100, 300],
'random_state': [2525],
'n_jobs': [1],
'min_samples_split': [3, 5, 10, 15, 20, 25, 30, 40, 50, 100],
'max_depth': [3, 5, 10, 15, 20, 25, 30, 40, 50, 100]
}
gsr = GridSearchCV(
RandomForestRegressor(),
search_params,
cv=3,
n_jobs=-1,
verbose=True
)
gsr.fit(tr_x, tr_y)
print(gsr.best_estimator_)
"""
# 決定係数R2でモデルの性能評価
print(f'R2_score: {model_2_score:.4f}')
#モデルの保存
pickle.dump(model_2, open("mymodel_RDF.pkl", "wb"))#アンサンブル
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('dataset.csv')
number_of_carbons(df)
types_of_chemical_bonds(df)
percentage_of_chemical_bonds(df)
# データセットから説明変数と回帰対象の変数を取り出し(今回はWater Solubilityを回帰)
train_x = df.drop(['SMILES','Water Solubility'], axis=1)
train_y = df['Water Solubility']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=7)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
kf_meta = KFold(n_splits=5, shuffle=True, random_state=687)
tr_idx_meta, va_idx_meta = list(kf_meta.split(train_x))[0]
tr_x_meta, va_x_meta = train_x.iloc[tr_idx_meta], train_x.iloc[va_idx_meta]
tr_y_meta, va_y_meta = train_y.iloc[tr_idx_meta], train_y.iloc[va_idx_meta]
# 特徴量と目的変数をxgboostのデータ構造に変換
dvalid = xgb.DMatrix(va_x, label=va_y)
dvalid_meta = xgb.DMatrix(va_x_meta, label=va_y_meta)
model_xgb = pickle.load(open("mymodel_xgb.pkl", "rb"))
model_xgb_2 = pickle.load(open("mymodel_xgb_2.pkl", "rb"))
model_xgb_3 = pickle.load(open("mymodel_xgb_3.pkl", "rb"))
model_xgb_4 = pickle.load(open("mymodel_xgb_4.pkl", "rb"))
model_xgb_5 = pickle.load(open("mymodel_xgb_5.pkl", "rb"))
model_RDF = pickle.load(open("mymodel_RDF.pkl", "rb"))
model_LR = pickle.load(open("mymodel_LR.pkl", "rb"))
xgb_pred = model_xgb.predict(dvalid)
xgb_pred_2 = model_xgb_2.predict(dvalid)
xgb_pred_3 = model_xgb_3.predict(dvalid)
xgb_pred_4 = model_xgb_4.predict(dvalid)
xgb_pred_5 = model_xgb_5.predict(dvalid)
RDF_pred = model_RDF.predict(va_x)
LR_pred = model_LR.predict(va_x)
xgb_pred_meta = model_xgb.predict(dvalid_meta)
xgb_pred_2_meta = model_xgb_2.predict(dvalid_meta)
xgb_pred_3_meta = model_xgb_3.predict(dvalid_meta)
xgb_pred_4_meta = model_xgb_4.predict(dvalid_meta)
xgb_pred_5_meta = model_xgb_5.predict(dvalid_meta)
RDF_pred_meta = model_RDF.predict(va_x_meta)
LR_pred_meta = model_LR.predict(va_x_meta)
average_meta = (xgb_pred_meta + xgb_pred_2_meta + xgb_pred_3_meta + xgb_pred_4_meta + xgb_pred_5_meta + RDF_pred_meta + LR_pred_meta)/7
weighted_average_meta = (xgb_pred_meta * 2 + xgb_pred_2_meta + xgb_pred_3_meta + xgb_pred_4_meta + xgb_pred_5_meta + RDF_pred_meta + LR_pred_meta)/7
variance_meta = (average_meta-xgb_pred_meta)**2 + (average_meta-xgb_pred_2_meta)**2 + (average_meta-xgb_pred_3_meta)**2 + (average_meta-xgb_pred_4_meta)**2 + (average_meta-xgb_pred_5_meta)**2 + (average_meta-RDF_pred_meta)**2 + (average_meta-LR_pred_meta)**2
stdev_meta = variance_meta**0.5
stack_pred = np.column_stack((xgb_pred_meta,
xgb_pred_meta**2,
xgb_pred_2_meta,
xgb_pred_2_meta**2,
xgb_pred_3_meta,
xgb_pred_3_meta**2,
xgb_pred_4_meta,
xgb_pred_4_meta**2,
xgb_pred_5_meta,
xgb_pred_5_meta**2,
RDF_pred_meta,
RDF_pred_meta**2,
LR_pred_meta,
average_meta,
weighted_average_meta,
stdev_meta
))
scaler = MinMaxScaler()
stack_pred_scaled = scaler.fit_transform(stack_pred)
meta_model = LinearRegression(normalize=True)
meta_model.fit(stack_pred_scaled, va_y_meta)
average = (xgb_pred + xgb_pred_2 + xgb_pred_3 + xgb_pred_4 + xgb_pred_5 + RDF_pred + LR_pred)/7
weighted_average = (xgb_pred * 2 + xgb_pred_2 + xgb_pred_3 + xgb_pred_4 + xgb_pred_5 + RDF_pred + LR_pred)/7
variance = (average-xgb_pred)**2 + (average-xgb_pred_2)**2 + (average-xgb_pred_3)**2 + (average-xgb_pred_4)**2 + (average-xgb_pred_5)**2 + (average-RDF_pred)**2 + (average-LR_pred)**2
stdev = variance**0.5
stack_test_pred = np.column_stack((xgb_pred,
xgb_pred**2,
xgb_pred_2,
xgb_pred_2**2,
xgb_pred_3,
xgb_pred_3**2,
xgb_pred_4,
xgb_pred_4**2,
xgb_pred_5,
xgb_pred_5**2,
RDF_pred,
RDF_pred**2,
LR_pred,
average,
weighted_average,
stdev
))
stack_test_pred_scaled = scaler.transform(stack_test_pred)
meta_pred = meta_model.predict(stack_test_pred_scaled)
score = R2(meta_pred, va_y)
# 決定係数R2でモデルの性能評価
print(f'R2_score: {score:.4f}')
##モデルの保存
pickle.dump(meta_model, open("mymodel_meta.pkl", "wb"))
pickle.dump(scaler, open("scaler.pkl", "wb"))Q3. SMILES記法で表された分子のエームズ試験の結果を予測せよ。
こちらの問題はデータセットが与えられず、自分でデータを集めてきて、与えられたSMILE記法に対してのエームズ試験の結果（1:陽性、0:陰性）の結果を予想するモデルを作成するタスクでした。
取り組むに当たり、エームズ試験のデータは以下のリンクから持ってきました。
Benchmark Data Set for In Silico Prediction of Ames Mutagenicity
Q3.の解法
SMILES記法とエームズ試験の結果だけでは特徴量がなく有効的なモデルの作成が出来ないので、ネット情報を参考にしながら個々のSMILES記法に対して200個以上の物性値を出力できるプログラムを作成しました。
今回の場合も、Q2と同じく勾配ブースティング木でアプローチを行い、その後線形回帰、ランダムフォレストなどでモデルの作成を行いました。
こちらの問題では最終順位が秘密裏のデータセットを用いて行われるということで、汎用的なモデルを作成できるのかどうかに焦点が当てられていて、個人的にはQ2よりも難易度は高かったように思えます。
以下、ソースコード例になります。（どのソースコードで最終的なスコアをとれたのかが分からないため、勾配ブースティング木を用いたモデルを掲載します。PL:0.871）
#データセットの作成
import pandas as pd
from rdkit import Chem
from rdkit.Chem import Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors
#各種データを格納するための空のリストを作製
mols_list =[]
smiles_list = []
activity_list = []
#Activity.csvの読み込み
df = pd.read_csv('Activity.csv')
#SMILES記法から作製したMOLオブジェクトがNoneではなかった場合、MOLオブジェクトとエームズ試験の結果をリストに追加
for i in range(len(df)):
mols = Chem.MolFromSmiles(df.at[i, "SMILES"])
if mols is not None:
mols_list.append(mols)
smiles_list.append(df.at[i, "SMILES"])
activity_list.append(df.at[i, "Activity"])
#Rdkitを用いて各種物性値の計算
names = [x[0] for x in Descriptors._descList]
names = names[:200]
calc = MoleculeDescriptors.MolecularDescriptorCalculator(names)
descs = [calc.CalcDescriptors(mol) for mol in mols_list]
#新規のPandas.DataFrameオブジェクトを作製し、データセット（datasetsQ3.csv）を作製し出力
df = pd.DataFrame(descs, columns=names, index=smiles_list)
df.head()
df["Activity"] = activity_list
df.to_csv('datasetsQ3.csv', sep=',')#モデル
import xgboost as xgb
from sklearn.metrics import log_loss
from sklearn.model_selection import KFold
import sys
import pickle
import pandas as pd
from rdkit import Chem
# データセットの読み込み
#各有機化合物のCの数をデータとして追加
df = pd.read_csv('datasetsQ3.csv')
# データセットから説明変数と回帰対象の変数を取り出し(今回はActivityを回帰)
train_x = df.drop(['SMILES','Activity'], axis=1)
train_y = df['Activity']
# 学習用と評価用にデータを分割(クロスバリデーション)
kf = KFold(n_splits=5, shuffle=True, random_state=71)
tr_idx, va_idx = list(kf.split(train_x))[0]
tr_x, va_x = train_x.iloc[tr_idx], train_x.iloc[va_idx]
tr_y, va_y = train_y.iloc[tr_idx], train_y.iloc[va_idx]
# 特徴量と目的変数をxgboostのデータ構造に変換
dtrain = xgb.DMatrix(tr_x, label=tr_y)
dvalid = xgb.DMatrix(va_x, label=va_y)
#xgbのハイパーパラメータの設定
#各値はベイズ最適化でパラメータ探索を行った
params = {'objective': 'binary:logistic',
'eta': 0.01,
'gamma': 0.0,
'alpha': 0.0,
'lambda': 1.0,
'min_child_weight': 1.0,
'max_depth': 10,
'subsample': 0.8,
'colsample_bytree': 0.8,
'random_state': 77}
num_round = 600
# モデルの学習
# バリデーションデータもモデルに渡し、学習の進行とともにスコアがどう変わるかモニタリングを行う
watchlist = [(dtrain, 'train'), (dvalid, 'eval')]
model = xgb.train(params, dtrain, num_round, evals=watchlist)
# バリデーションデータでのスコアの確認
va_pred = model.predict(dvalid)
score = log_loss(va_y, va_pred)
# 対数損失でモデルの性能評価
print(f'log_loss: {score:.4f}')
##モデルの保存
pickle.dump(model, open("mymodel.pkl", "wb"))初めてデータコンペに参加してみた感想
今回初めてデータコンペに参加してみて感じたのは、「社会に出たときに必要な実践的な技術、アプローチが学べる」という点です。
もちろん、実際の予算や期間を考えるとこのように1か月単位で精度を上げるのは現実的ではないと思うのですが（1ヵ月取り組んだが、0.78→0.84と向上したのは6%程度）、どの切り口でモデルを選択すればよいのか、モデルの作り方、どういう特徴量が精度向上に有効なのか、主要ライブラリの使い方（pandas、numpy、scikit-learn、keras、xgboost…）などなど、もし将来データ分析に携わる職を志望しているのであれば、参加して得られるものは計り知れないなと感じました。
私自身、将来はマテリアルズ・インフォマティクスによるデータドリブンの材料開発を主導できる研究者になり、環境問題、エネルギー効率の問題に取り組んでいくことを一つの選択肢として考えているので、データ分析の技術＋α、実際にマテリアルズ・インフォマティクスを学べたので、本コンペは貴重な経験になりました。
全くデータ分析をしたことがない方でも（私自身もそうでした）、紹介させてもらったKaggle本のような書籍を読みながら進めていけばデータ分析の本質に触れながら本格的なモデルを作れるようになるので、まだ参加したことがないけど興味はある！という方は是非一度参加してみてください！（習うより慣れよ、っていう言葉がぴったりだと思います）
今後はこの経験を糧に様々なデータコンペに参加してみてデータ分析技術を学びつつ、精度の高いモデル作成を試みていきたいと思っています！！
P.S. 出来るだけ書籍等を参考にしながら間違いのないように努めましたが、データ分析初心者でまだまだ未熟なこともあり記事内容に不備があるかもしれません。理解を深めるために、もし内容に不備があった場合コメント欄にてその点ご指摘いただけると幸いです。また、記事の感想などありましたら合わせてコメント欄に書いていただけると嬉しいです！！
とても長文になってしまいましたが、ここまで読んでくださった方本当にありがとうございました！
ENHANCER代表 みでっと
|
|