Ichijoji Sagarimatsu

統計モデリング 生態学のデータ解析・麻雀の解析など

【天鳳】天鳳位なるまでシミュレーションを回してみた【閲覧注意?】

天鳳位になるまで打ち続けた場合」にどれだけの試合数がかかるのか気になったので調べてみました。

私の知る範囲ではこれについて調べた記事は見つからなかったので、自作関数を作ってシミュレーションを回してみました。

記事の前半では天鳳位なるまでチャレンジ後半では十段なるまでチャレンジの結果を紹介します。


設定としては、7段原点(1400pt)を開始地点として、(実力としての)安定段位が9.5段・9段・8.5段・8段の4通りについて、天鳳位(十段)になるまでのシミュレーションを1000回まわし、何半荘かかったかを調べました。

なお、順位分布は各安定段位になるような等差順列(1>2>3>4型)を想定しました。

天鳳位チャレンジのシミュレーション

天鳳位チャレンジ:安定9.5段(.284-.261-.239-.216)

中央値で2890半荘、平均にすると3739半荘かかりました。
中央値で2890半荘というのは、この実力の人たちを集めて打たせたら2890半荘以内に半数が天鳳位になってるということです。
逆に言うと残り半数は2890半荘時点で天鳳位になれていないということです。
実力安定9.5段というとウルトラ鉄強なんですが、意外に時間がかかりますね。イヤな予感がしてきました。

天鳳位チャレンジ:安定9段(.278-.259-.241-.222)

中央値で5936半荘、平均で8145半荘かかりました。
安定9段も超鉄強なんですが、沼ると20000半荘以上かかるようです。
そういったケースに引っ張られて、中央値と平均値の乖離は大きくなっています。

天鳳位チャレンジ:安定8.5段(.271-.257-.243-.229)

中央値で17243半荘、平均で25722半荘かかりました。
安定8.5段も十分に鉄強です。

天鳳位チャレンジ:安定8段(.265-.255-.245-.235)

中央値で89517半荘、平均で125743半荘かかりました。
一方で100人に1人は2000半荘程度で天鳳位をとれるようです。

十段チャレンジのシミュレーション

では初めて十段になるまでには何半荘かかるでしょうか。
ここからは十段チャレンジの結果を紹介します。

十段チャレンジ:安定9.5段(.284-.261-.239-.216)

中央値で1018半荘、平均で1198半荘かかりました

十段チャレンジ:安定9段(.278-.259-.241-.222)

中央値で1418半荘、平均で1767半荘かかりました

十段チャレンジ:安定8.5段(.271-.257-.243-.229)

中央値で2351半荘、平均で3044半荘かかりました

十段チャレンジ:安定8段(.265-.255-.245-.235)

中央値で5754半荘、平均で7592半荘かかりました

おわりに:天鳳との向き合い方

改めて天鳳の闇が深すぎることがわかりました。

安定9段といえば超鉄強なんですが、それでも下振れると20000半荘天鳳位を取れない一方で、安定8段が上振れれば2000半荘で取ってしまうようです。

どれだけ麻雀の腕に自信があっても、天鳳位獲得を目標にするのは気が狂ってしまう可能性があるためお勧めできません。

私も一応、一鳳凰民として天鳳位は目指していますが、自分で結果を見て正直ゲンナリしてしまいました。

まぁ、(良い意味でも悪い意味でも)所詮はゲームということですね。

純粋に天鳳ではレベルの高い麻雀が(サクサク)打てる、という点に天鳳の魅力があることは変わらないと思います。

また一方で、天鳳位を目標にするのは現実的ではないが、ある程度の実力があれば十段は達成できることもわかりました

ちなみに、各安定段位における十段→天鳳位チャレンジの成功率は以下のようになっています。

安定8.5段あれば、(十段になれさえすれば)10回に1回は天鳳位を取れるようです。

成功率に差はあれど、やはり十段は天鳳位ガチャのチケットと考えるのがよさそうですね。

この記事を読んでいただいた天鳳民がどう思うかはわかりませんが、私としては、まだ十段にもなったことはないので、とりあえず十段は目指して頑張りたいと思います。

①十段を目指す(ちゃんと実力があればある程度の試行回数で達成できそうであるため)
②数%の可能性にかけて天鳳位ガチャを目指す
③降段したら...その時考える

って感じでしょうか。

なお、本記事は1>2>3>4型の着順分布の結果のみ示しましたが、1=2=3>4型や2=3>1=4型も試してみました。
後者2つの分布では、1>2>3>4型と比べて、天鳳位達成までに若干時間がかかるようでした(中央値で300~400半荘くらい)。

最後に、ここまで書いておいてなんですが、本記事の執筆にあたっては誰の校正も受けていないので、私のシミュレーションのコードが間違えている可能性もあります。
Rのコードになりますが、再現性のため載せておきます。

for(i in 1:nsimu){
  dani <- 7      # 初期段位
  point <- 1400  # 初期ポイント
  loop <- 0      # 試合数の初期設定
  
  while (dani <= 10 & dani >= 4){   # 4段から10段の間は以下を繰り返す
    if(dani == 10){
      while (point >= 0 & point < 4000){
        tmp <- rmultinom(1, 1, juni)
        tmp.point <- 90*tmp[1] + 45*tmp[2] - 180*tmp[4]
        point <- point + tmp.point
        loop <- loop + 1
        cat(i, loop, dani, "段", point, "\n")
      }
    } else if (dani == 9){
      while (point >= 0 & point < 3600){
        tmp <- rmultinom(1, 1, juni)
        tmp.point <- 90*tmp[1] + 45*tmp[2] - 165*tmp[4]
        point <- point + tmp.point
        loop <- loop + 1
        cat(i, loop, dani, "段", point, "\n")
      }
    } else if (dani == 8){
      while (point >= 0 & point < 3200){
        tmp <- rmultinom(1, 1, juni)
        tmp.point <- 90*tmp[1] + 45*tmp[2] - 150*tmp[4]
        point <- point + tmp.point
        loop <- loop + 1
        cat(i, loop, dani, "段", point, "\n")
      }
    } else if (dani == 7){
      while (point >= 0 & point < 2800){
        tmp <- rmultinom(1, 1, juni)
        tmp.point <- 90*tmp[1] + 45*tmp[2] - 135*tmp[4]
        point <- point + tmp.point
        loop <- loop + 1
        cat(i, loop, dani, "段", point, "\n")
      }
    } else if (dani == 6){
      while (point >= 0 & point < 2400){
        tmp <- rmultinom(1, 1, juni)
        tmp.point <- 90*tmp[1] + 45*tmp[2] - 120*tmp[4]
        point <- point + tmp.point
        loop <- loop + 1
        cat(i, loop, dani, "段", point, "\n")
      }
    } else if (dani == 5){
      while (point >= 0 & point < 2000){
        tmp <- rmultinom(1, 1, juni)
        tmp.point <- 90*tmp[1] + 45*tmp[2] - 105*tmp[4]
        point <- point + tmp.point
        loop <- loop + 1
        cat(i, loop, dani, "段", point, "\n")
      }
    } else if (dani == 4){
      while (point >= 0 & point < 1600){
        tmp <- rmultinom(1, 1, juni)
        tmp.point <- 90*tmp[1] + 45*tmp[2] - 90*tmp[4]
        point <- point + tmp.point
        loop <- loop + 1
        cat(i, loop, dani, "段", point, "\n")
      }
    }
  
    # 各段位における昇段降段シミュレーションを終えて段位変動
    if (point < 0){
      dani <- dani - 1
    } else {
      dani <- dani + 1
    }
  
    # 変動後段位に応じてポイントを原点に設定
    if (dani == 10){
      point <- 2000
    } else if (dani == 9){
      point <- 1800
    } else if (dani == 8){
      point <- 1600
    } else if (dani == 7){
      point <- 1400
    } else if (dani == 6){
      point <- 1200
    } else if (dani == 5){
      point <- 1000
    } else if (dani == 4){
      point <- 800
    } 
  }
  
  loop_all[i] <- loop
  dani_end[i] <- dani
}

【麻雀】【天鳳】麻雀の「スタイル」の解析

はじめに

本記事は、「麻雀のスタイル」なるものを、主成分分析というデータ分析手法を用いて明らかにしようという記事です。

本記事では 麻雀におけるプレイヤースタイルは存在するという前提で話を進めたいと思います。

解析の流れ

データには、天鳳鳳凰卓において、長期で好成績を残しているプレイヤー(いわゆる"鉄強")36ID(一部同一人物の別アカウントを含む)のスタッツを用いました。

スタッツには、成績に影響があると考えられる以下の12項目を用いました。

・副露率
・流局時聴牌率:流局した際に聴牌している確率
和了
和了平均得点:和了時の平均打点
・副露平均得点:副露で和了した際の平均打点
和了順目
・ツモ率
・リーチ率
・リーチ先制率
・リーチ成功率
・放銃率
・放銃平均得点:放銃時の平均放銃打点

スタッツの取得にはのどっちさん(https://nodocchi.moe/tenhoulog/)を利用させていただきました。

前提として、取り上げたプレイヤーは皆鉄強であり、どのスタイルが良い(悪い)という話ではないということを理解しておいてもらえればと思います。

解析には主成分分析という手法を用いました。

結果

見方は後で説明するとして、とりあえず、結果の図を載せます。

12変数を2次元に圧縮し、プレイヤーをプロット

今回の検証には12のスタッツを用いています。
本来12の変数を図示するには12次元が必要なのですが、それを2次元まで圧縮して図示したものが上図になります。
当然、本来12次元の情報を2次元に図示するにあたっては、いくらか情報が削ぎ落とされていますが、その損失をできるだけ抑えるように上手くやってくれるのが主成分分析になります。
今回の解析においては、上の図は元々の情報の66.6%を保持しています

図に見方についての疑問点は
・縦軸と横軸はなんなの?
・矢印は何?
の2点かと思います。

軸について

横軸と縦軸は、12の変数(スタッツ)を適当に重み付けして合成された軸になります。
スタッツの矢印は、軸方向に与えている影響の大きさを示します。

横軸について

横軸に最も大きな影響を与えているスタッツは、「和了平均得点」であり、ついで「副露率」「副露平均得点」、その次に「和了率」「和了順目」でした。

すなわち、「和了平均得点」「副露平均得点」が大きいプレイヤーは図の右側に、「副露率」が高いプレイヤーは図の左側にプロットされることを示しています。
矢印が逆方向にあるのは、それぞれのスタッツが負の相関を持っていることを示しており、副露率が高いプレイヤーは打点が小さくなることを反映していると考えられます。
どうやら横軸は「副露率と打点」を示していると考えられそうです。

図の右側:面前型・低和了率・高打点
図の左側:副露型・高和了率・低打点

縦軸について

主成分分析において、横軸と縦軸は完全に無相関であることを意味します。

すなわち、横軸が「スタイル①:副露型 / 門前型」を示していると考えれば、縦軸は、それとは全く独立に取りうるスタイル②について示唆していると捉えられます。

縦軸に大きな影響を与えているスタッツは、「リーチ率」「放銃率」「流局聴牌率」ついで「放銃平均得点」「リーチ成功率」「和了率」でした。

どうやら、図の縦軸は「立直率」「放銃率」「和了率」を示していそうです。
「リーチ率」「放銃率」「流局聴牌率」「放銃平均得点」「和了率」が高いプレイヤーは図の上側に、それらが低い代わりに「リーチ成功率」は高いプレイヤーが図の下側にプロットされると考えられます。

図の上側:リーチ型・高和了・高放銃
図の下側:ダマ型・低和了・低放銃

まとめ

再掲

以上より、麻雀のスタイルは

門前リーチ型(図の右上領域)
門前ダマ型(図の右下領域)
副露リーチ型(図の左上領域)
副露ダマ型(図の左下領域)

で分類できるということがわかりました。

繰り返しますが、面白いのは横軸と縦軸は無相関であるということです。
実際、元データの相関関係を見てみると、「副露率」と「リーチ率」の相関係数は0.047であり、限りなく無相関でした。

これはすなわち、「門前リーチ派」「門前ダマ派」「副露リーチ派」「副露ダマ派」はどれもスタイルとして取り得るということです

もしも牌譜等を見て勉強するときは、スタイルが違いすぎる人よりも、ある程度自分とスタイルが近い人を参照すればいいのではないかと思います。

自分のスタッツを入力してどこにプロットされるかを示すWebアプリを実装しようかと考えましたが、一旦見送りました...もしも需要があれば作ってみたいと考えるので、気軽にコメント等くれるとうれしいです

ちなみに私は「副露ダマ型」でした。

考察(2024/3/10追記)

縦軸が「立直率」「放銃率」「和了率」を反映しているのは間違いなそうですが、それを一纏めにリーチ型 / ダマ型 と分類することに関しては議論の余地があるかもしれません。

というのも、今回はいわゆる"鉄強"ばかりを集めているので、基本的にどのプレイヤーもリーチするべき手はリーチし、ダマるべき手はダマっているはずです。

とすると、縦軸は、単に「いまこの手をリーチするか」というリーチ判断を示しているのではなく、そこに至るまでのプロセスを反映していると考えるのが妥当かもしれません。

たとえば、先制を取られた時に攻撃的な判断をするか、あっさりやめる守備的な判断を取るか、といったことです。
前者は追いついた場合に追っかけ立直することも増えるでしょうから、結果的に「立直率」「和了率」「放銃率」はどれも上昇しそうです。

その場合、縦軸は「攻撃型」か「守備型」かを反映していると捉えることもできそうです。

だとすると結果的には麻雀のスタイルは

門前攻撃型(図の右上領域)
門前守備型(図の右下領域)
副露攻撃型(図の左上領域)
副露守備型(図の左下領域)

で表されるということになります。

この表現は麻雀界隈ではよくある表現ですが、それが統計的に正しい分類であるということがわかるのも面白いですね。

鳳凰卓東南戦の真の安定段位をベイズ推定

本記事は、以下のWebアプリに関する補足記事となります。

tenhou-antei.netlify.app

本Webアプリの開発に関しては、理論面の他、パラメータの設定など、
一方通行さんのブログ記事 【麻雀】安定段位をベイズ推定する|一方通行
を参考にさせていただきました。

ここでは、概念的な説明と、結果の解釈などについて説明します。

ベイズ推定の利点

本Webアプリで実装しているベイズ推定は、天鳳における実力の推定に非常に相性が良いと考えています。

その理由は、ベイズ推定ではデータとして用いるプレイヤーの成績の他に、「事前情報」を統計的に正しく利用できるためです

ベイズ的な手法を用いない場合、推定には手元のデータを用いるしかありません。
その場合、500戦の安定段位が13段の人の真の実力段位は、「13段がもっともありそうだ」と推定されます。

しかし、鳳凰卓の人であれば、安定13段は流石に上振れだろうと感じると思います。
長期安定段位は高くても10段前後が限界であり、11段を超える人はほぼいないと考えられます。

このような事前情報を、ベイズ統計であれば適切に反映させることができます。

ベイズ推定の仕組み

ベイズ統計では、
①事前確率(事前情報)を与える
②データを与える
③データが事前確率を更新し、事後確率(推定結果)を返す

の流れで推定が行われます。

ベイズ推定のイメージ

実際に、今回の解析の流れを図示して例示します。

①事前情報とは、鳳凰卓プレイヤーの着順分布に関する事前の信念を指します。長期でトップ率35%やラス率15%はいないだろうといったことです。実際には着順分布について仮定を敷いていますが、安定段位の指標に置き換えて考えても構いません。
以下の図が今回与えている安定段位に関する事前分布です。
事前分布の95%区間は5.58 ~ 8.72段となっています。これは、任意のプレイヤーAを選んだ時、Aの実力段位は95%の確率で5.58段 ~ 8.72段の間にあり、可能性としては7段周辺のことが多いかもしれない、と予想している状態です。実力段位が9段を超えている可能性も1%程度残されています。
この事前分布の妥当性については後で検討します。

②データとしてプレイヤーAの1000戦の成績[270-260-240-230](安定8.43段)を与えます。

③事前確率は更新され、事後確率を返します。以下が更新された事後分布です(事前分布と重ね合わせています)。事後分布の95%区間は6.70 ~ 8.99段になりました。

ベイズ統計では、サンプルサイズによる影響は「不確実性」として適切に評価されます

与えるデータが多いほど、事前確率は大きく更新され、推定結果の確実性が増します。
与えるデータが少なければ、事前情報からの更新は少なく、事後確率は事前確率と大きくは変わりません。

※ここで言う「確実性」や「不確実性」は、事後分布の幅(95%推定区間の広さなど)について言及したものであり、推定結果が「正しくない」という意味ではありません。


したがって、よく「500戦程度の安定段位はアテにならない」などと言われますが、本手法を用いることで、「どれだけアテにならない」のか定量的に評価できます。


なお、本解析では実力が一定であるという仮定のもと推定を行っています。

実際には実力が向上することも考えられるため、恣意的でさえなければ全成績を入力する必要はないと考えています
鬼打ち勢の方は直近数千戦の結果を入力することも推奨します。
その分データが少なくなるため、推定結果の不確実性が増すことは避けられませんが、「現在の実力」の推定が可能と考えます。

結果の解釈

以上を実装するに当たって問題となるのは、どのような事前分布を与えるかです。

ただ、その話題に移る前に、結果の解釈についての補足説明をしたいと思います。

例として[1400-1400-1400-1300] の成績を入力した結果

推定中央値

推定中央値に関する解釈は、ザックリ言えば、

解釈① 実力としてもっともありそうな安定段位

になります。
ただし、この解釈は語弊を招くかもしれません。
というのも、「相対的にそのあたりの確率密度が高い」だけであり、それ以外の可能性(それより低い可能性も高い可能性も)大いにあるからです。
ですので、ここでは別の表現による解釈を紹介したいと思います。

解釈② あなたの長期安定段位は50%の確率で[推定中央値]以上になる

こちらの方が、より正確というか、語弊を招かないかと思います。
言い換えれば、あなたが新しく天鳳アカウントを作り、同じ実力のままプレイしなおした時、長期(この場合は数千戦と考えて構いません)安定段位は50%の確率で [推定中央値] を超える、という意味になります。

95%推定区間

事前分布と入力データから割り出された
あなたの真の実力段位が95%の範囲で含まれる区間 を意味します。

上側2.5%点と下側2.5%点から割り出される数値になりますので、
推定区間の下限値に関しては「97.5%以上の確率であなたの真の実力段位は [下限] 以上である」と言えます。

したがって、これが7段を超えている人は、
紛れも無い鳳凰、あるいは 統計的に有意な鳳凰 と言うことができます。
これはなかなか誇れることではないかと思います。

鳳凰の猛者的なやつ

私が独自に決めた分類に従って表示しています。
95%区間の下限値、中央値や打数などを勘案したものになります。
あくまでも私がザックリ決めた分類なので、必ずしも正確な実力を示すものではありません
95%区間の下限値などを参照しているので、打数の影響も大きいでしょう
(たとえば打数が数百戦だと "豆鳳" 認定確率は高いと思われます)。
つまり遊び心なので、あまり真に受けず楽しんでもらえればと思います。

事前情報に関する設定

実装に関して具体的にどのような事前分布(事前情報)を与えたかを説明します。
これに関しても、 【麻雀】安定段位をベイズ推定する|一方通行 に明瞭な説明がありますが、ここでも説明します。

事前確率に関しては、ディリクレ分布と呼ばれる確率分布のパラメータ  \alpha を指定することで制御しています。

ディリクレ分布というのは、よく「サイコロの出目の出やすさ」を決める確率分布と説明され、ここでは着順分布を決める確率分布になります。

 \alphaには、一方通行さんに倣い \alpha=200を指定しました。
これが具体的にどのような事前情報かというと、本条件でシミュレートされるプレイヤーの真の安定段位の分布が先ほどの図になります(再掲)。

7段をピークとして、上位1%の人間で実力安定9.07段、上位0.1%の人間で実力安定9.85段となっています。

実際に鳳凰卓東南戦を2000戦以上プレイしたプレイヤー1087人(2024年3月時点)の安定段位の分布(出典:男冥利さん http://otokomyouri.com/toppage.aspx)を図示すると以下になります。

鳳凰卓を2000戦以上打てている時点である程度強い人が多いといった背景も考えられるため、平均が7以上になっていたり、下が切れているような分布になっていますが、上限は概ね一致しており、事前分布としては上々かと思われます。

(正確には、上の水色のヒストグラムは事前分布からシミュレートされた実力段位の分布であり、黄緑のヒストグラム実際に2000戦プレイした結果の分布なので、水色の実力の人たちが2000戦プレイすると結果の分散はより大きくなりそうですが、ここでは大体OKということにしたいと思います)

まとめると、事前分布としてディリクレ分布のパラメータには \alpha = 200を指定しました。
これは「95%のプレイヤーの実力が5.6段〜8.7段に収まり、上位1%の実力が9段以上」と指定していると考えてもらえればと思います。

さいごに

いかがでしたでしょうか。
面白いと思ってもらえたり、自信・やる気につながってもらえれば嬉しいです。

途中にも書きましたが、個人差はあれど、麻雀も実力は向上するものだと思います
なので、本解析で満足のいく結果が得られなくても、今後上手く打つことを意識して、ある程度の打数を重ねた時、今の自分の実力のチェックに本サービスを利用してもらえればと思います。

また、「95%推定区間が広すぎてほとんど何もわかんないじゃん!」と思った人は、自分のしているゲームはそういうゲームであることを認識してもらいつつ、打数も増やして確実性の高い結果を得てもらえればと思います。

個人的には天鳳が一番面白い麻雀ゲームだと思っています。

皆が運の有識者になってくれることを願いつつ、楽しく遊んでいけたらと思います。

NAGA記録供養

私はNAGA解析にかけた対局のNAGA度・悪手率を記録していました。
といっても、何か目的意識があったわけではなく、なんとなく習慣として続けている程度のものだったのですが... 。

データもそこそこ溜まってきたので、これを機にデータを整理し、有効活用して供養したいと思います。

調べたこと

① 過去9ヶ月間のNAGA度、悪手率の推移を可視化しました
ラスの後は打牌のクオリティが落ちるかを検証しました
③ 日ごとのコンディションの良し悪しは存在するかを検証しました
ポイント状況で打牌のクオリティが変わるかを検証しました

NAGA記録データについて

今回はNAGAの守備型(ガンマタイプ)を参照した 2023年6月 ~ 2024年2月の記録をデータとして使用します。
ちなみに、私が守備型のガンマタイプを参照している理由は、ガンマが一番懐が広い(NAGA度が高く出る)気がするためです。

(整形したデータ↓)

2023年6月以降、私は鳳凰卓東南戦を1100半荘プレイしており、そのうちNAGA解析にかけたのは451半荘でした。

検証

検証① 月ごとのNAGA度、悪手率の推移

早速結果を図示します。NAGA度については線形モデル(LM)、悪手率については一般化線形モデル(GLM)で近似直線(曲線)も同時に図示しています。

NAGA度

悪手率

全期間での平均NAGA度は91.5、平均悪手率は2.92% でした。

......まぁ、フーンって感じですね。

近似線によると、統計的に有意な実力の変化は見られませんでした。
実力が向上しているという結果になることを望んでいましたが、まぁ麻雀ですので、そんな簡単に上手くはなれないということでしょう。

しかし、心の目で見ればNAGA度は上昇傾向、悪手率は減少傾向にあります(実際近似線はわずかにその傾向を示しています)。

分析手(ぶんせきて)としてはこの程度の傾向はランダムのブレによる可能性が高いと言わざるを得ませんが、打ち手(うちて)としてはそこに意味を見出したくなってしまいます。

検証② ラスのあとは打牌のクオリティが落ちるか

ここでは、前の対局の着順が次の対局の内容に影響を及ぼしているかを調べました。

たとえば、ラスを引いた次の半荘はNAGA度が低くなるのではないか、ということです。

一方で、「寝たら気分はリセットされる」と考え、その日初めてプレイする対局は、「前の着順なし」として別に分けて考えました。
こちらも、NAGA度と悪手率についてそれぞれ検証しました。

NAGA度

悪手率

データ数が多くはないこともあり、統計的に有意な結果は得られませんでした。
しかし意外にも、ラスを引いた次の対局の内容が悪くなる、という結果は得られないどころか、NAGA度・悪手率ともにむしろ内容が良くなっている傾向が見られました

検証③ 日ごとのコンディションの良し悪しは存在するか

ここでは日ごとのコンディション、というものが存在するかを調べました。

「なんか今日は上手く打ててる気がする」「今日はがよく見える...」的な話はあるとも聞きますが、NAGA度・悪手率に日ごとの影響は検出されるでしょうか。

この節では図らしい図はないので、代わりに技術的な話を書いておこうと思います(興味ない人は飛ばしてください)。


以下のモデルⅠとモデルⅡのどちらがより優れているかを調べました。

モデルⅠ : 日ごとのコンディションによる影響がないモデル
このモデルでは簡単に、すべての対局が独立試行であると仮定します。
(NAGA度を指標とした場合には)対局結果のNAGA度( y)が真の実力としてのNAGA度( \mu)から標準偏差 \sigma正規分布誤差をもって生成されると考えます(簡単に言えば、単なる正規分布モデルです)。

 y \sim Normal (\mu, \sigma)


モデルⅡ : 日ごとのコンディションによる影響があるモデル

モデルⅠではすべての試行を独立としましたが、ここにその日の実力  \mu_{\ day} の階層を挿入します。
すなわち、真の実力としてのNAGA度( \mu)からその日の実力  \mu_{\ day} が生成され、さらにそこから対局結果のNAGA度( y)が生成される、というモデルです(簡単に言えば、正規分布×2です)。

 y \sim Normal (\mu_{\ day}\ , \sigma)

 \mu_{\ day} \sim Normal (\mu, \sigma_{\ day}\ )

どちらが優れているかの判断には、WAICと呼ばれる情報量基準を用いました。
また、悪手率を指標とした場合には、悪手率は「打牌選択回数における悪手回数」の割合データなので、正規分布ではなく二項分布をもちいましたが、それ以外は同じです。

結果↓


WAICは値が小さいほど優れたモデルであることを示します。

したがって、今回の検証では、NAGA度、悪手率どちらの指標を用いた場合でも、
「日ごとのコンディションの良し悪し」は検出されませんでした

とはいえ、NAGA度については、両モデルのWAICの差は0.2程度であり、大きな差はありませんでした。
WAIC(AIC)の概念的な観点から解釈すれば「その日のコンディションを考慮したモデルではほんの少しだけ予測精度が向上したが、モデルを複雑化するほどではなかった」といったところでしょうか。

検証④ ポイント状況で打牌のクオリティが変わるか

ここではポイント状況と打牌選択の質の関係を検討しました。

打ち手としては、体感的に、ポイントが3桁になってくると萎え萎えに、折り返し超えるとイケイケドンドンな気持ちになる気がするので、そこを基準に分類しました。

(結果の図↓)

NAGA度

悪手率

データ数の問題もあるので、統計的に有意な差はありませんでしたが、やはりポイントを持ってるイケイケドンドンモードのときは内容がいい気がしなくもない、という結果になりました。

ポイントを持っているときは、NAGA度にして平均0.7、悪手率にして平均0.6% の変化が見られるので、意外に馬鹿にならない差かもしれません。

まとめ

予想はしていましたが、データの数が450程度ですと、やはり統計的な指標から確度をもって断言できる事柄はほとんどなく、
結果的には図をまとめただけみたいになってしまったのは残念です。

一方で、これは分析手(ぶんせきて)としては残念なことですが、打ち手(うちて)としては素晴らしいことかもしれません。

というのも、今回の検証でその日の体調や気分で打牌が大きく変わることはなく、仮に前の対局でラスを引いていようとも、厳しいポイント状況にあろうとも、メンタルをブレさずに淡々と打牌を続ける理想の打ち手の姿が浮かび上がったと言えなくもないからです。

これからもNAGA記録はとっていこうと思うので、今回書いたコードを元に、半年くらいを目処に内容を更新していきたいと思います。

【天鳳】安定段位の収束の可視化

安定段位のシミュレーションをしたので、その結果を載せます。

このテーマ自体はもう何番煎じか分からないような話題ですが、結果が分かりやすく可視化されているのは意外と見ない気がするので、その点で多少は需要があると思い記事にします。

安定段位について、シミュレーションから分かったことを結論から書くと、

・1000戦で1.5〜2.0段、5000戦で0.7〜0.8段 ほどブレる

・強者ほどブレ幅が大きい

・一般に上振れの方が下振れよりも幅が大きい

あたりです。

順番に見ていきます。

この手のシミュレーションは色々な仮定が敷けると思いますが、今回の検証の設定は以下です。


鳳凰卓想定


・真の実力(安定段位)として、1〜4位まで等差順列な着順分布を想定


・X回(X = 100, 400, 1000, 2000, 3000, 4000, 5000, 7000, 10000, 20000)ゲームを行い安定段位を計算する。これを各Xに対して10000回繰り返す。

真の実力:安定七段(.25 - .25 - .25 - .25)

塗りつぶした区間の上限が、10000回の試行のうち安定段位が上位2.5%を、
塗りつぶした区間の下限が、10000回の試行のうち安定段位が下位2.5%を意味します。

また、中央の線は10000回の試行の中央値で、期待される「真の実力」にほぼ一致することが分かります。

たとえ20000戦打ったとしても ±0.3段はブレるようです。

真の実力:安定八段(.265 - .255 - .245 - .235)

真の実力:安定九段(.278 - .259 - .241 - .222)



塗りつぶした区間の上限と下限、すなわち、上振れ2.5%点下振れ2.5%点が、期待される真の安定段位と比べてどれだけブレているのかを図にしたのが以下です。

つまり、各ゲーム数で最大どれくらいの上振れと下振れが起こるか、を示した図です。

安定七段

安定八段

安定九段

なんだか同じような図ばかりで面白くないですが、
後半3枚の図は、真の安定段位が高くなるほどグラフのY軸も大きくなっていることに注目してください。

これはつまり、実力がある人ほど上振れ/下振れのスケール感も大きくなっているということです。

以上3枚の上振れ/下振れ図の結果をまとめた表を載せておきます。

ngamesはゲーム数を、upper_lowerは上振れか下振れかを、
そしてantei7 〜 antei9の列の値が、真の実力が安定7・8・9段の人の、それぞれのゲーム数における、期待安定段位からのブレ幅の上振れ/ 下振れ2.5%点です。

自分の真の実力(は計りようがないですが)とゲーム数でどれくらいブレている可能性があるかの早見表のような感じでお使いください(※ ただしシミュレーションなのでブレがあります)。


表およびグラフより、真の安定段位が高い人は、低い人に比べてブレ幅が大きいことが分かります。

たとえば、1000戦を基準とすると、

真の実力が安定七段の人が、 + 1.52段の上振れ - 1.28段の下振れ を引きうるのに対し、

真の実力が安定九段の人は、 + 1.94段の上振れ - 1.58段の下振れ を引き得ます。


安定七段そこそこで天鳳位を目指している私(やその他多くの人たち)からすると、

デカ目の上振れは実力のあるものに訪れる

という少し悲しい事実が判明したことになります。


あと、常に上振れの方が下振れよりもスケールが大きくなっているのは、
「ラス回避」を連続で続けることで安定段位が天元突破することがあるからだと思われます。

短い期間であればラスなしの安定段位 + Inf段が現れることなどからも、上振れの方がスケールが大きくなることは体感的に分かるのではないでしょうか。


また、安定段位はどこまでいっても安定しないことが分かりました。

収束具合の傾きを見た目で判断する感じでは、2000 ~ 5000戦 くらいが実力を測る上でコスパがいい(?)ですかね。

実力七段なら、3000戦で ± 0.8段ほど、5000戦で ± 0.6段くらいです。

...ってまぁそんな打数ですら普通は打てないですけどね。
というか打たない方がいいと思います。

この記事は以上です。

【麻雀】運:実力 = 50 : 1 ? 着順の統計モデリング

麻雀は実力と運が勝敗を左右するゲームですが、運と実力の比率はどのくらいなのでしょうか?

7:3くらい? 8:2くらいでしょうか? あるいは100:1?

RとStanを用いて、この問いに答えてみたいと思います。

使用データ

データには、オンライン対戦麻雀【天鳳】の公開しているログ(https://tenhou.net/sc/raw/)による、プレイヤーの着順履歴を使用します。

サンプルサイズを絞るため、今回は2018年〜2020年の鳳凰卓4人打ち東風戦のデータを利用します。*1

さらにサンプルサイズを効率的に絞るため、 対戦履歴を「3年間で合計500ゲーム以上プレイしているプレイヤー4人による試合」に限ることとします。
すなわち、鳳凰卓常連 の4人による対戦のみに注目します。*2

整形したデータは以下のようになっています。

3年間で東風戦を500戦以上プレイした全226プレイヤーによる、計40657試合の履歴です。

モデルの説明

モデルの概要を説明します。

  1. 各プレイヤーの潜在的強さがあります。
    麻雀が完全実力ゲームであれば、この潜在的な強さに従って、この場合、1位から4位が「Aさん→Bさん→Cさん→Dさん」の順となります。

  2. しかし麻雀は完全実力ゲームではないので、運要素が発生します。
    運要素は、各プレイヤーに独立に、正規分布に従うものとして発生すると仮定します。

  3. 実際に各プレイヤーに正規分布乱数としての運要素が付加されます。

  4. こうして実力要素運要素の両方から、最終的な着順が決まります。
    この例では、「Bさん→Aさん→Dさん→Cさん」の順位になりました。

最後に、いくつか条件を付け加えます。
まず、各プレイヤーの強さも、何らかの正規分布にしたがっていると仮定します

このとき、「強さ」とは相対的なものなので、この正規分布の平均はゼロに固定できます。

また、先の例において、運要素、実力要素の数字のオーダーは不定です。
(運要素10に対して実力要素1と、運要素1に対して実力要素0.1、は同じことです)

ですので、仮定した運要素の正規分布標準偏差を1に固定します。
そして「各プレイヤーの強さの正規分布」の標準偏差をXとおき、これを推定します。

すなわち以下のようなイメージです。
麻雀には実力でどうこうできる標準偏差Xの正規分布と、実力ではどうにもならない運としての標準偏差1の正規分布、二つの要素から勝敗が決定する、ということです。

これにより「麻雀は運 : 実力 = 1 : X である」と推定できます。

なお、実際のStanコードは、松浦健太郎氏による書籍の、第10章のものを使用しました。

結果

モデルの収束には12時間ほどかかりました。
実力の分布である正規分布標準偏差Xは、中央値で0.0183と推定されました(下図)。

標準偏差X(s_muと宣言)の推定値

よって、麻雀の運と実力の比は、1:0.0183 ... すなわち、およそ50:1であると分かりました。

解釈

結果の解釈に際して、いくつか注意すべきと思われる点を列挙します。

「勝負ムラ」を「運要素」と捉えている(完全実力ゲームなど存在しない?)

将棋を例に挙げます。

将棋は運要素が極端に少ないことは、多くの人が同意する点であると思います。

ほとんど将棋のルールを知っているだけの私が、羽生善治さんと対戦しても、一度も勝てないでしょう。

しかし、藤井聡太さんと羽生善治さんならどうでしょうか?
今回使用したモデルに合わせて考えると、仮に藤井聡太さんの強さを49、羽生善治さんの強さを51(逆でも構いません)とすると、「完全実力ゲーム」なる将棋では、100回やって100回羽生善治さんが勝たねければなりません。

しかし、これはありそうもないことです。

この説明として、先に紹介した松浦健太郎さんの書籍では、「勝負ムラ」という単語が使われています。
つまり、その日(あるいはその対局)におけるパフォーマンスが、二人の強さ(49と51)を中心にして確率的に振る舞うと考えることで解決します。

これを運要素と捉えるかどうかは、何とも言い難い問題であるような気がしますが、
本モデルではこれを「運要素」として一纏めにしていることには注意しておきたいと思われます。

そもそもレベルの高いフィールドでの話をしている点

これは、先の将棋の例と繋がります。
本解析では、天鳳鳳凰卓という、一定の雀力がないと対戦権を得られないフィールドでの履歴を参照しました。

しかし実際には、特上卓・上級卓・一般卓という「鳳凰卓以下」のフィールドも存在します。

なにが言いたいかというと、レベルの高いフィールド(正確には、同程度のレベルの集合)では、「勝ったり負けたり」が発生しがちである...すなわち、「運要素」が大きくなるということです。

将棋においても、私が最上位層のプロ棋士に挑んでも一向に歯が立たないでしょうが、
最上位層同士の争いであれば、「勝ったり負けたり」が発生し、
「運要素(これを運要素と捉えるかは問題ですが)」が大きくなることが分かるでしょう。

この点で、本当はさまざまなレベルのプレイヤーが一纏めに対戦する場の履歴を解析すれば良かったのですが、
ネット麻雀ではレベルによる卓分けが明確になされている分、対処は難しかったです。

いろんなレベルの人たちが集まるフリー麻雀の対戦履歴が得られれば一番いいですね。

そのため、今回の解析結果の解釈としては、
「麻雀」というよりも「鳳凰卓東風戦」の運と実力比が50:1であると捉えておくのが無難です。

「個人ごとのパフォーマンスのブレの違い」を考慮していない

今回、運要素を全て「平均ゼロ・標準偏差1」に固定しましたが、実際には、プレイヤーごとにいくらかブレ幅に違いがあることが予想されます。

たとえば、将棋と同様に「その日その対局におけるコンディション」による違いや、
あるいは採用しがちな戦術によるブレ幅もあるでしょう(たとえば「トップラス麻雀」の人はブレ幅が大きいと言えます)

ここではそのようなプレイヤーごとの違いをモデル化していない(正確には、全て同じであるとしてモデル化している)ことにも、一応注意しておきたいところだと思います。

終わりに

以上です。
今回得られた結果から安定段位*3に関するシミュレーションもできたので、次はそれを書いてみます。

質問・コメント等あれば是非お願いします。

8/31 追記

思ったよりも多くの方に見られて緊張しているので、追記です。

このモデルのポイントは

運要素に正規分布を仮定する
プレイヤーの実力に正規分布を仮定する
プレイヤーはそれぞれの実力に応じた「下駄」を履いた状態で
正規分布乱数としての運要素が独立に付加され
最終的に着順が決まる

というものです

この4つ目の仮定に関しては、実は現実的ではないかもしれません。
なぜかというと、麻雀は「誰かが点棒を失えば、誰かが点棒を得る」からです。

すなわち、各プレイヤーに作用している運要素は、実際には独立ではないことが考えられます。
この点を改良したモデルを作れないかなぁと考えているので、できたら追記します。

モデルの適合の良さについては言及していなかったので、さらっと書いておきます。

推定されたパラメータ0.0183を用いて正規分布から発生させた色々な強さを持つ1000プレイヤーについて、総当たり、一人およそ2000対戦させ、
1000人分の安定段位を集計した結果が以下です。

男冥利さん(http://otokomyouri.com/toppage.aspx)というサイトからとってきた、実際の鳳東プレイヤーの安定段位の分布が以下です。

実力の標準偏差を0.05と0.1でシミュレートすると以下のようになります。
横軸の値に注目してください。


実力の幅を過大評価しているせいで、現実には誕生し得ないほどの安定段位を持った仮想プレイヤーが爆誕しています。
なお、安定段位の低い方の端は適合が悪そうですが、これは現実にはプレイヤーが鳳凰卓での対戦権を失ってしまい、2000ゲームという対戦数を重ねられないことが原因と考えられます。

以上から、簡易的な評価ですが、このモデルはまだ改良の余地がありそうですが、適合のほどは悪くないと思います。

Rコード

######################### 2019-2020年分の鳳東データの前処理 ##############################
# データの読み込み
files2018 <- list.files("data/2009.2.20_2021.8.3鳳凰卓scc", pattern = "scc2018", full.names = TRUE)
files2019 <- list.files("data/2009.2.20_2021.8.3鳳凰卓scc", pattern = "scc2019", full.names = TRUE)
files2020 <- list.files("data/2009.2.20_2021.8.3鳳凰卓scc", pattern = "scc2020", full.names = TRUE)

d_orig2018 <- do.call(rbind, lapply(files2018, read_delim, delim=" | ", quote="", 
                                col_names = c("X1","X2","X3","X4","X5","X6","X7","X8","X9","T1","T2","T3","T4")))
d_orig2019 <- do.call(rbind, lapply(files2019, read_delim, delim=" | ", quote="", 
                                col_names = c("X1","X2","X3","X4","X5","X6","X7","X8","X9","T1","T2","T3","T4")))
d_orig2020 <- do.call(rbind, lapply(files2020, read_delim, delim=" | ", quote="", 
                                col_names = c("X1","X2","X3","X4","X5","X6","X7","X8","X9","T1","T2","T3","T4")))

d_orig <- rbind(d_orig2018, d_orig2019, d_orig2020)

# 四鳳南喰赤 を選択
d <- filter(d_orig, X5 == "四鳳東喰赤-" | X5 == "四鳳東喰赤")[,10:13]

# 不要な文字列を削除
d$T1 <- str_remove(d$T1, pattern = "\\(\\+.*?\\)")
d$T2 <- str_remove(d$T2, pattern = "\\(\\+.*?\\)")
d$T2 <- str_remove(d$T2, pattern = "\\(\\-.*?\\)")
d$T2 <- str_remove(d$T2, pattern = "\\(0.0\\)")
d$T3 <- str_remove(d$T3, pattern = "\\(\\+.*?\\)")
d$T3 <- str_remove(d$T3, pattern = "\\(\\-.*?\\)")
d$T3 <- str_remove(d$T3, pattern = "\\(0.0\\)")
d$T4 <- str_remove(d$T4, pattern = "\\(\\-.*?\\)")
d$T1 <- str_remove(d$T1, "<br>")
d$T2 <- str_remove(d$T2, "<br>")
d$T3 <- str_remove(d$T3, "<br>")
d$T4 <- str_remove(d$T4, "<br>")

# 4,3,2,1着順に並び替え
d <- d[,c(4,3,2,1)]

#### データの把握 ####
# 全対戦数
nmatch <- nrow(d)

# 縦長データに変換
d_long <- pivot_longer(d, cols=everything(), names_to="order", names_prefix="T", values_to="player")
d_long$match <- rep(1:nrow(d), each=4)

# 全プレイヤーの総対戦数
M <- arrange(count(d_long, player), n)

colnames(d) <- c("forth", "third", "second", "first")

##### 対戦数少プレイヤーを省く #####

# 3年間で対戦数500以上のプレイヤー
sum(M$n>=500)   # 226player

# 対戦数500以上のプレイヤー名を抽出
P500 <- M$player[which(M$n>=500)]
res.p500 <- is.element(d_long$player, P500)   # 対戦数500以上のプレイヤーが含まれるか?

# 対戦数が500以上のプレイヤーが4人以上含まれる試合ナンバーを抽出
m500 <- data.frame(match = d_long$match[which(res.p500 == TRUE)])
tmp <- m500 %>% 
  group_by(match) %>% 
  filter(n()>3)
M500 <- unique(tmp$match)

# 対戦数500以上のプレイヤーが4人以上含まれる試合を選抜
d500 <- d[M500, ]

#' *対戦数500以上のプレイヤーが4人以上含まれる試合で解析*

# 個々人のID名に対応するナンバー
d500_long <- pivot_longer(d500, cols=everything(), names_to="order", names_prefix="T", values_to="player")
conv <- 1:length(unique(d500_long$player))
names(conv) <- unique(d500_long$player)

LtoW <- array(NA, dim=c(nrow(d500), ncol(d500)))
for(i in 1:nrow(d500)){
  for(j in 1:ncol(d500)){
    LtoW[i,j] <- conv[as.character(d500[i,j])]
  }
}

nmatch500 <- nrow(d500)
data <- list(nind=max(conv), nmatch=nmatch500, LtoW=LtoW)

#' *RUN*
params <- c("s_mu", "mu")
tonpu500 <- stan(file="analysis1.stan", data=data, seed=1234, pars=params,
                 thin=3, iter=6000, warmup=3000, chain=3)

Stanコード

data{
  int nind;
  int nmatch;
  int<lower=1, upper=nind> LtoW[nmatch, 4];
}

parameters{
  real<lower=0> s_mu;
  vector[nind] mu;
  ordered[4] performance[nmatch];
}

model{

  // 事前分布
  s_mu ~ normal(0, 0.3);

  for(m in 1:nmatch){
    for(k in 1:4){
      performance[m,k] ~ normal(mu[LtoW[m,k]], 1);
    }
  }

  mu ~ normal(0, s_mu);
}

*1:麻雀といえば東南戦という方も多いかもしれませんが、限られたプレイヤーで複数回対戦している東風戦の方が解析に都合が良かったため、こちらを使用しています

*2:オンライン対戦麻雀【天鳳】では、ポイントを賭けた段位制を導入しています。鳳凰は七段以上のプレイヤーが対戦を許されるフィールドであり、六段に落ちると鳳凰での対戦権を失います。すなわち、今回解析対象とする鳳凰卓常連は、一定以上のレベルが約束されたプレイヤーと言えます

*3:運要素の大きい麻雀において天鳳プレイヤーが唯一心の拠り所にしている麻雀の「長期成績」

【Stan】【農学】平均と標準偏差でガンマ分布を定義する:混合ガンマ分布による除草剤抵抗性雑草の割合推定

複数の確率分布を混合した分布(混合分布)で集団における各クラスタの割合推定をすることができます。

解析には2つのガンマ分布の混合分布を用いますが、その際、ガンマ分布のパラメータとして「平均・標準偏差」を用いるのが便利です。

今回は、農業雑草の研究をしている私のラボにおける仮想データから紹介します。

除草剤抵抗性雑草について

農業の現場において、作物の収量低下をもたらす雑草は、現在では除草剤処理により抑制することが多いです。

しかし、一つの除草剤の繰り返しの使用や、規定の薬量濃度に達しない低薬量散布などの要因により、 「除草剤抵抗性」形質を持つ個体が出現し、世界中で大きな農業被害をもたらしています。

ある個体が除草剤抵抗性の形質を有しているかを室内実験から判断するには、

除草剤抵抗性を有しないことが既知の(感受性の)対照群
ある水田圃場から採取した対象の個体群

のそれぞれに除草剤処理を施し、 t検定等を行い判断します。

すなわち、「抵抗性 / 感受性」の評価を一個体単位で行うのではなく、個体数に反復を取った水田圃場等の集団単位で行います。

これは、たとえばある個体が抵抗性(=抵抗性遺伝子をいくらか有している個体)であったとしても、
薬量処理濃度や、測定する評価指標(葉・根・重量etc.)によっては感受性と判断がつかないほどに生育不良を起こすことも多いためです。

データの分布

本題です。

今回は、上記のような個体群単位の情報が得られなかった場合を想定します。

たとえば、海外からの輸入穀物に雑草の種子がよく混入しているのですが、
これに含まれる除草剤抵抗性個体の割合を推定したい、とします。

ここにもはや圃場単位の情報は存在しないため、一個体ごとに抵抗性か感受性かを判断しなければなりません。


これらの輸入穀物混入種子を発芽させ、除草剤処理を施しました。
生育の評価指標として、「薬剤処理後の根の伸長量」を用いました。

200個体に対して実験を施し、得られた200行のデータの分布は以下のようになりました。

根の伸長量が小さい除草剤感受性型と、根の伸長量が大きい除草剤抵抗性型が混在することによる双峰型の分布になりました。

これを2つのガンマ分布の混合分布と見立て、除草剤抵抗性集団の割合を推定します。

なお、シミュレーションデータを生成した2つのガンマ分布の真値は以下のようにしました(コードはページ最後に載せています)。

パラメータはshapeパラメータとrateパラメータで指定しており、 図の二つの分布を、除草剤抵抗性個体(図の右の山)の割合aを0.7として混合しています。

モデルの記述

data{
  int N;          // サンプルサイズ
  vector[N] Y;    // 根の伸長量
}

parameters{
  real<lower=0, upper=1> a;   // 抵抗性個体の割合
  positive_ordered[2] shape;  // shapeパラメータ
  vector<lower=0>[2] rate;    // rateパラメータ
}

model{
  // 対数尤度の計算
  for(n in 1:N)
    target += log_sum_exp(
      log1m(a)   + gamma_lpdf(Y[n] | shape[1], rate[1]),  // 感受性個体
      log(a) + gamma_lpdf(Y[n] | shape[2], rate[2])       // 抵抗性個体
    );
}

感受性(図の左の山)のガンマ分布のパラメータを(shape1, rate1)
抵抗性(図の右の山)を(shape2, rate2)とし、
抵抗性の割合をaとしたとき、

ある一サンプルのデータが得られる尤度p

p = (1 - a) * Gamma(shape1, rate1) + a * Gamma(shape2, rate2)

であり、これの対数尤度が

log(p)

= log \{(1 - a) * Gamma(shape1, rate1) + a * Gamma(shape2, rate2)\}

= log [exp \{log(1-a) + log(Gamma(shape1, rate1)\} + exp \{log(a) + log(Gamma(shape2, rate2)\}]

ですので、
log_sum_exp関数の定義である

  log_sum_exp(x, y) = log (exp(x) + exp(y))

に従いmodelブロックで記述し、対数尤度を足し上げています。
なお、log1m(a) = log(1-a) です。

ラベルスイッチングを避けるため、正のパラメータであるshapeパラメータに shape[1] < shape[2] となる制約をpositive_ordered型で宣言することにより課しています。


これをRからキックすればよいわけですが、結果から言うと収束しません

トレースプロットを見ると、以下のようになっています。

f:id:IchijojiSagarimatsu:20210704174947p:plain

rate1が彼方に飛んでいってしまっているようです。

seed値や元データによってはshapeパラメータが彼方に飛んでいったりもします。

こういったあり得ない数値に飛ぶことを防止するため、
shapeパラメータとrateパラメータに範囲の制約を課せばいいのですが、
ガンマ分布のパラメータの解釈は少しわかりづらいという問題があります。

スクリプトを他者に見せる際にも、この値はどこから出てきたのかな?と思われそうでもあります。

そこで、より直感的なパラメータ、
すなわち平均と標準偏差を用いてガンマ分布を定義し直します

平均と標準偏差でガンマ分布を定義する

shapeパラメータα、rateパラメータβ のガンマ分布について、
平均mu標準偏差sigmaはそれぞれ
 mu = α / β  sigma = {\sqrt{α}}  /  β
で表されます。

これを変形すると、shapeパラメータとrateパラメータは、平均と標準偏差を用いて
 α = mu^2 / sigma^2  β = mu / sigma^2
となります。


これを用いてモデルを記述し、各パラメータに制約を入れます。

data{
  int N;          // サンプルサイズ
  vector[N] Y;    // 根の伸長量
}

parameters{
  real<lower=0, upper=1> a;           // 抵抗性個体の割合
  real<lower=0, upper=15> mu1;        // 感受性個体の平均
  real<lower=0, upper=100> mu2;       // 抵抗性個体の平均
  real<lower=0, upper=15> sigma1;     // 感受性個体の標準偏差
  real<lower=0, upper=100> sigma2;    // 抵抗性個体の標準偏差
}

transformed parameters{
  real shape1 = square(mu1 / sigma1);  // 感受性個体のshapeパラメータ
  real shape2 = square(mu2 / sigma2);  // 抵抗性個体のshapeパラメータ
  real rate1 = mu1 / square(sigma1);   // 感受性個体のrateパラメータ
  real rate2 = mu2 / square(sigma2);   // 抵抗性個体のrateパラメータ
}

model{
  // 対数尤度の計算
  for(n in 1:N)
    target += log_sum_exp(
      log1m(a)   + gamma_lpdf(Y[n] | shape1, rate1),
      log(a) + gamma_lpdf(Y[n] | shape2, rate2)
    );
}  

パラメータの制約なしではやはり収束しませんが、
parametersブロックに示したような制約を入れることで収束します。

ここでは、
「感受性個体の根の伸長量の平均は15mm以下であり、抵抗性個体は100mm以下、 またそれぞれのバラツキもそれより大きくなることはない」
という意味の制約を課しています。

shapeパラメータとrateパラメータを直接記述するよりも、
平均と標準偏差を用いた定義の仕方では、
より直感的な解釈が可能になっていると分かります

今回は単にモデルを収束させるための範囲制限でしたが、
たとえば平均と標準偏差を用いて事前分布を記述することも容易になります。

ちなみに、推定結果の事後分布は以下のようになります。

f:id:IchijojiSagarimatsu:20210704190035p:plain

すべてのパラメータはその95%区間に真値を含んでおり、良い精度で推定できています。
除草剤抵抗性個体の推定割合aは、事後分布の中央値で69%であると評価されました。

参考

StanとRでベイズ統計モデリング 松浦健太郎
 混合分布やlog_sum_exp関数について詳しく説明されています。

Rコード

最後に、使用したRコードを載せておきます。

library(rstan)
library(ggplot2)
set.seed(100)

# パラメータ真値の設定
a <- 0.7
N <- 200
shape1 <- 2
shape2 <- 5
rate1 <- 0.7
rate2 <- 0.2

 # 混合前の分布図示
curve(dgamma(x, shape1, rate1), 0, 60, ylab="", xlab="")
curve(dgamma(x, shape2, rate2), 0, 60, add=T)

y1 <- rgamma(N*(1-a), shape1, rate1)
y2 <- rgamma(N*a, shape2, rate2)

d <- data.frame(Y=c(y1, y2)) # シミュレーションデータセット

# データの把握
ggplot(d, aes(x=Y, y=..density..))+
  theme_test()+
  geom_histogram(fill="red", alpha=0.5)

# Stanのキック
data <- list(N=nrow(d), Y=d$Y)

# shapeパラメータ、rateパラメータによる記述
fit1 <- stan(file="params_default.stan", data=data, seed=123)
stan_trace(fit1)  # 収束しない

# 平均と標準偏差からパラメータを指定
fit2 <- stan(file="params_mean_sd.stan", data=data, seed=123)
stan_trace(fit2)
stan_hist(fit2)