生物物理屋が数理最適化でブラックジャックに勝つ手を探してみた話 後編

## はじめに
この記事は「[ただただアウトプットを癖付けるための Advent Calendar 2024](https://qiita.com/advent-calendar/2024/newbie2024)」に投稿した記事です。

[最初の記事](https://qiita.com/yotabouya/items/e63c53c5e18b3cf6a454#%E3%81%AF%E3%81%98%E3%82%81%E3%81%AB)にも書いた通り、私は生物物理の実験を専門にしている研究者です。
最近はデータ解析のため機械学習のコード開発も行っており、幸いにもその成果が[NeurIPSに採択されました](https://neurips.cc/virtual/2024/poster/93465)。

前回から学習のため、今回はブラックジャックでの行動を最適化するという課題に取り組んでいます。
今回の記事では、前回定式化した問題を実際にpythonで実装し、最適化を行ってみた結果を報告します。

## 関連記事

#### 前の記事 「[生物物理屋が数理最適化でブラックジャックに勝つ手を探してみた話](https://qiita.com/yotabouya/private/8808d6872452d2f3ac5e)」

#### 次の記事「[生物物理屋が数理最適化でブラックジャックに勝つ手を探してみた話 後編その2](https://qiita.com/yotabouya/private/60c3dfcac8ad04ce2372)」

## 実装
次に、定式化した数理最適化問題を実装します。
実装したものを[こちらのリポジトリ](https://github.com/MasahitoUWAMICHI/optimize_blackjack)に置いています。よければご覧ください。

### ライブラリの選定
まずはpythonの数理最適化ライブラリを選定します。
ライブラリの一覧は[こちら](https://qiita.com/np_hsgw/items/0a7b55935bf944da8fc6)を参考にしました。

今回の問題設定は非線形かつ離散的な問題です。
したがって、非線形最適化ライブラリを選定する必要があります。

今回の問題設定だと、以下のような方針で最適化をおこなうことになります。
- 状態空間を定め、そこでの勝敗を計算する関数を定義する
- バースト状態とディーラーのターンにおける方策を与える
- それ以外の状態について、方策を与える
- 状態価値関数を計算する
- 評価関数を計算する
- 評価関数を最大化するような方策を求める

非線形であり、やや複雑な問題設定になるため、pyomoやGekkoあたりが適していると考えられます。

と、本来ならこのような選定になるのですが、いきなりgekkoを書こうとして挫折しました。
そのため、今回は比較的使い慣れているscipyで、まずブラックジャックの実装をしっかりおこないました。
また、連続最適化の方が慣れているので、方策は0から1の連続値とします。

### 方針

まずブラックジャック自体の実装です。

確実に必要な情報は以下でしょう。
- 手札に各数字のカードが何枚ずつあるか
- その手札は何点か
- バーストしているか
- 勝敗
- カードを引く行動の処理
- 引くカードの確率の設定
- スタンドの処理
- ディーラーの行動の処理
- 初期状態の設定

また、最適化を考えると以下も必要です。
- 方策を踏まえた、各手札状態での勝率の計算
- 評価関数の設定
- 最適化における変数の設定

### 手札に各数字のカードが何枚ずつあるか

これは、各点数のカードの枚数を整数で格納するベクトルにより表現しました。
正規のブラックジャックでは1-10点(Aは一旦すべて1点としてカウント)のカードがあるので、プレイヤーとディーラーそれぞれに、10次元のベクトルが与えられることになります。

### その手札は何点か

これが意外に面倒でした。
基本的には上で定義した手札枚数のベクトルに対応して、点数のベクトルを定義し、内積を取れば点数が出ます。これを基本点とします。
その後にAを処理する必要がありますが、これは次のようにおこないました。
- バーストの点数と基本点の差を取る
- A1枚を11点にするときの差分(11-1=10点)で上の差を割り、小数点以下を切り捨てる
    - こうすることで、手持ちのAの枚数とは関係なく、バーストせず11点にできるAの枚数の最大値が求まる
- 手持ちのAの枚数と比較し、小さい方を取る
- もし負の値であれば0にする
- 求めた枚数に、10点をかける
- 基本点に足す

### バーストしているか

これは素直に、バースト点である21を超えているか判定すれば良いです。
バーストした場合、ただちに勝敗が返されます。バーストした側の負けですね。

### 勝敗

勝敗は0と1で表現します。プレイヤーの勝利を1としておくことで、勝敗の期待値をそのまま勝率として扱えます。
バースト以外での勝敗は、手札の点数の大きい方の勝利です。
便宜上、点数が等しいときはプレイヤーの勝利としておきます。

### カードを引く行動の処理

カードを引く場合、引かれるカードの確率分布を参照します。
そして、手札ベクトルの該当するカードの枚数に1を足した状態を参照します。
今回は引かれるカードについてすべての可能性を考えます。
参照先での勝率と、参照先への遷移確率の内積を取ることで、現状態での勝率を再帰的に計算することになります。

### 引くカードの確率の設定

現状でプレイヤーとディーラーが持っている各種カードの枚数をその総数から引いた残り枚数を規格化することで確率を計算します。

### スタンドの処理

プレイヤーのターンであるか、ディーラーのターンであるかを記述する変数を用意しておきます。
初期状態ではプレイヤーのターンにしておき、スタンドしたらディーラーのターンにします。
ディーラーのスタンドに対しては、その時点での勝敗を返すこととします。

### ディーラーの行動の処理

ディーラーは、現時点での勝敗を見て行動します。
プレイヤーが勝っている状況ではカードを引きます。
ディーラーが勝っている状況ではスタンドします。

### 初期状態の設定

初期状態については特殊な方策を採用します。
ディーラーの手札が0枚のときは、ディーラーがカードを引きます。
プレイヤーの手札が0または1枚のときは、プレイヤーがカードを引きます。
これにより、初期状態の確率分布も考えることができます。

### 方策を踏まえた、各手札状態での勝率の計算

これは上でも触れてしまいましたが、カードを引くことによる手札状態の遷移確率を踏まえて再帰的に勝率を定義します。
方策はカードを引く確率として扱うので、カードを引く場合の勝率(遷移先の勝率の加重平均)とスタンドする場合の勝率(その手札状態からディーラーの行動を考えることにより決まる)をそれぞれ計算し、これらによる期待値を計算するわけです。

### 評価関数の設定

評価関数は、両者が１枚もカードを持たない状態からの勝率として定義しました。
これにより、初期状態(プレイヤーが2枚、ディーラーが1枚のカードをオープンに持つ)の確率分布を踏まえた、ゲーム自体への勝率を最適化の対象とすることができます。

### 最適化変数

変数は、各手札状態での方策となります。
これを実装するには、まずどのような手札状態があるか確認する必要があります。
そのため、まず、両者手札ゼロの状態から、カードを引くorスタンドする場合の遷移先を再帰的に探索し、現れた状態をdictに登録していきました。
この中で、初期状態実現までのステップ、バースト状態とディーラーのターンを除いた各状態を各要素に一対一で紐付けたnumpy.ndarrayを用意し、これを変数としました。
また、最適化の計算量を減らすため、引くカードの確率分布や引いたあとの手札状態、手札の点数など、方策に関係ないものについてはこの時点で計算し、dictに保存しておきました。

## 結果

正規のブラックジャックで以上の処理を試したところ、長い時間がかかってしまいました。
そのため、ひとまずはカード枚数を減らしました。
4点までのカードのみ用意し、A,2,3は4枚ずつ、4点のカードは8枚ある状況とします。
また、Aは都合により5点を取ることができ、10点を超えたらバーストです。
このように設定すると、数十分で最適化が終わりました。
結果を確認すると、0や1という、はっきりした方策が現れました。
手札の点数が多くなってきたらカードは引かない方が良い、点数が少ないうちは引いたほうが良いという現実が、ある程度反映されたと思います。
しかし、依然として0.5くらいのどっちつかずのものもありました。
カードを引くorスタンドの選択では、差は小さくとも優勢なものはあるはずなので、最適化がそこまで及ばなかったと見るのが良さそうです。