概要

『あつまれどうぶつの森』におけるカブの売値（たぬき商店での平日の値段）は9～660ベルの間で変動するとされ、特に最高値660ベル、最安値9ベルが出るのは非常に珍しいとされている。
そこで、ある週において、最高値および最安値が出現する確率を計算してみる。

前提

求める確率は、ゲームプレイ週数を無限大に近づけた際に収束する値とする。
以下のカブ価決定アルゴリズムの解析結果（2020年7月13日時点）が正しいものとする。
https://gist.github.com/Treeki/85be14d297c80c8b3c0a76375743325b
ゲーム内の時刻を現実のものから変える行為（いわゆる時間操作）は行わないものとする。
擬似乱数に規則性を見いだし、特定の値になるように調整する行為（いわゆる乱数調整）は行わないものとする。

カブ価変動パターン

カブ価の変動パターンは大きく以下の4種類に分けられる。

減少型: 単調減少し、買値を超えることはない。
波型: 何度か上下を繰り返し、買値の近辺の値をとり続ける。
跳ね大型: 単調減少ののち増加に転じ、大きなピークを迎え、減少する。
跳ね小型: 単調減少ののち増加に転じ、小さなピークを迎え、減少する。

最高値は「3. 跳ね大型」でのみ、最安値は「4. 跳ね小型」でのみ発生する。
そのため、まずは各型の確率を求めてみる。

各型になる確率は、一つ前の週の型によってのみ決定され（マルコフ連鎖）、以下の表のようになっている。（解析結果のコード135～209行目）

前週\今週	減少型	波型	跳ね大型	跳ね小型
減少型	5%	25%	45%	25%
波型	15%	20%	30%	35%
跳ね大型	20%	50%	5%	25%
跳ね小型	15%	45%	25%	15%

この表から、下記のサイトの方法で、ゲームプレイ週数を無限大に近づけた際に各型が収束する確率（定常分布）を求めてみる。
（本当は収束することの証明が必要であるが、ここでは省略する。）
https://withcation.com/2020/01/23/%E3%83%9E%E3%83%AB%E3%82%B3%E3%83%95%E9%80%A3%E9%8E%96%E3%81%AE%E5%AE%9A%E5%B8%B8%E5%88%86%E5%B8%83%E3%82%92%E8%A7%A3%E3%81%8F/

減少型、波型、跳ね大型、跳ね小型の確率がそれぞれ $a, b, c, d$ になるとすると、以下の方程式が成り立つ。

\begin{pmatrix}
a & b & c & d
\end{pmatrix}
\begin{pmatrix}
0.05 & 0.25 & 0.45 & 0.25 \\
0.15 & 0.20 & 0.30 & 0.35 \\
0.20 & 0.50 & 0.05 & 0.25 \\
0.15 & 0.45 & 0.25 & 0.15 \\
\end{pmatrix}
=
\begin{pmatrix}
a & b & c & d
\end{pmatrix}

移項して、

\begin{pmatrix}
a & b & c & d
\end{pmatrix}
\begin{pmatrix}
-0.95 & 0.25 & 0.45 & 0.25 \\
0.15 & -0.80 & 0.30 & 0.35 \\
0.20 & 0.50 & -0.95 & 0.25 \\
0.15 & 0.45 & 0.25 & -0.85 \\
\end{pmatrix}
=
\begin{pmatrix}
0 & 0 & 0 & 0
\end{pmatrix}

一見、4種の文字を使った4式の連立方程式を行列表記したようにみえるが、実は任意の3式を変形すると残り1式になってしまうため、式が一つ足りない（4x4の行列の行列式は0になり、逆行列が求められない）。
そのため、 $ a+b+c+d=1 $ を任意の列（ここでは1列目）に埋め込む。

\begin{pmatrix}
a & b & c & d
\end{pmatrix}
\begin{pmatrix}
1 & 0.25 & 0.45 & 0.25 \\
1 & -0.80 & 0.30 & 0.35 \\
1 & 0.50 & -0.95 & 0.25 \\
1 & 0.45 & 0.25 & -0.85 \\
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0
\end{pmatrix}
\\
\begin{pmatrix}
a & b & c & d
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
1 & 0.25 & 0.45 & 0.25 \\
1 & -0.80 & 0.30 & 0.35 \\
1 & 0.50 & -0.95 & 0.25 \\
1 & 0.45 & 0.25 & -0.85 \\
\end{pmatrix}^{-1}
=
\begin{pmatrix}
0.1476 & 0.3463 & 0.2474 & 0.2588
\end{pmatrix}

各型の確率が求められた。Pythonのコードを以下に示す。

import numpy as np

A = np.array(
    [[0.05, 0.25, 0.45, 0.25],
    [0.15, 0.20, 0.30, 0.35],
    [0.20, 0.50, 0.05, 0.25],
    [0.15, 0.45, 0.25, 0.15]]
    )

A -= np.identity(4)  # 移項
A[:, 0] = 1  # 条件a+b+c+d=1を入れる
B = np.array([1, 0, 0, 0])

print(B @ np.linalg.inv(A))  # [0.1476074  0.34627733 0.24736279 0.25875248]

最高値（660ベル）となる確率

最高値（660ベル）になるためには、「3. 跳ね大型」になること（確率0.2474）に加え、以下の条件が必要である。

A: その島の買値（日曜午前にウリから買える値段）は90～110ベルであるが、これが110ベルである。
B: ピーク時には買値の2～6倍になるが、これが6倍になる。

Aについてはランダムであり（コード123行目のrandint）、確率は$\frac{1}{21}$である。
Bについては、まず倍率として2.0～6.0のfloat型の数値を生成するが、ピーク時の値段は（買値）×（倍率）の結果の小数点以下を切り上げたものであるため、ちょうど6.0倍でなくてもよい可能性がある。
倍率は、uint32_t型（ $0$ ～ $2^{32}$ の整数値をとる）の乱数をもとに決められ、乱数の値が大きいほど大きくなる。
そこで、ピーク時の値段が660ベルになるuint32_t型乱数のボーダーを二分探索を用いて調べることで、条件Bを満たす確率を求める。

#include <bits/stdc++.h>
using namespace std;

// uint32_t型の乱数をもとに、float型の数を返す関数。乱数が大きいほどbに近い数となる。
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

// 小数点以下を切り上げる関数
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

int main() {
  uint32_t ll = 0;  //uint32_tの下限値
  uint32_t rr = pow(2, 32) - 1;  //uint32_tの上限値
  int32_t basePrice = 110;  // 買値
  uint64_t l, r, m;
  int32_t peakPrice;

  l = 0;
  r = pow(2, 32) - 1;
  m = (l + r) / 2;

  // 二分探索
  while (r - l > 1) {
    peakPrice = intceil(randfloat(2.0, 6.0, m) * basePrice);  // ピーク時の値段
    if (peakPrice < 660) {l = m;}
    else {r = m;}
    m = (l + r) / 2;
  }

  cout << 0 << " から " << rr << " の値をとる乱数が " << m << " より大きい場合、最高値となる。" << endl;
  cout << "その確率は " << (double)(rr - m) / (double)(rr) << " である。" << endl;
}

出力

0 から 4294967295 の値をとる乱数が 4285206015 より大きい場合、最高値となる。
その確率は 0.00227273 である。

求めた確率すべてを掛け合わせると、 $0.2474 \times \frac{1}{21} \times 0.002272 = 2.68 \times 10^{-5}$ となる。
0.00268%と低い確率だが、ゲームの販売本数は2020年5月7日時点で1177万本であるため、カブ価をチェックしているのがそのうち1%だったとしても、毎週世界中で3人くらいは最高値660ベルを確認している計算になる。

それでは、次に最安値の確率も計算してみる。

最安値（9ベル）となる確率

最安値（9ベル）になるためには、「4. 跳ね小型」になること（確率0.2588）に加え、以下の条件が必要である。

A: その島の買値（日曜午前にウリから買える値段）は90～110ベルであるが、これが90ベルである。
B: 「初めて増加に転じる時期」が {月曜午前, 月曜午後, ... , 木曜午後} の8つのうち、月曜午前または木曜午後である。
※なお、Bで月曜午前の場合と木曜午後の場合の確率は等しく、プログラムもほぼ同じであるため、以降のC,Dの条件は木曜午後の場合で説明する。
C: 以上の条件の下で、月曜午前のカブ価は36～81ベルの可能性があるが、これが36ベルである。
D: 以上の条件の下で、木曜午後のカブ価は9～65ベルの可能性があるが、これが9ベルである。

A, Bについては整数の乱数で選ばれるため、それぞれ $\frac{1}{21}, \frac{1}{4}$ の確率となる。
Cについては、まず倍率として0.4～0.9のfloat型の数値を生成し、買値の90ベルに掛けた結果が月曜午前のカブ価となる。
小数点以下は切り上げであるため、36ベルとなるためには、倍率は0.4（誤差の関係上、わずかに上回る分には問題ない）でなければならない。
前述したように倍率は、uint32_t型（ $0$ ～ $2^{32}$ の整数値をとる）の乱数をもとに決められるが、ここでは乱数の値が小さいほど倍率が大きくなるプログラムとなっている（コード317行目）。
やはり、カブ価が36ベルとなる乱数のボーダーを二分探索を用いて調べることで、条件Cを満たす確率を求めるのがよさそうだ。

#include <bits/stdc++.h>
using namespace std;

// uint32_t型の乱数をもとに、float型の数を返す関数。乱数が大きいほどbに近い数となる。
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

// 小数点以下を切り上げる関数
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

int main() {
  uint32_t ll = 0;  //uint32_tの下限値
  uint32_t rr = pow(2, 32) - 1;  //uint32_tの上限値
  int32_t basePrice = 90;  // 買値
  uint64_t l, r, m;
  int32_t monAmPrice;

  l = 0;
  r = pow(2, 32) - 1;
  m = (l + r) / 2;
  
  // 二分探索
  while (r - l > 1) {
    monAmPrice = intceil(randfloat(0.9, 0.4, m) * basePrice);  // 月曜午前の値段
    if (monAmPrice > 36) {l = m;}
    else {r = m;}
    m = (l + r) / 2;
  }

  cout << ll << " から " << rr << " の値をとる乱数が " << m << " より大きい場合、条件Cが満たされる。" << endl;
  cout << "その確率は " << (double)(rr - m) / (double)(rr) << " である。" << endl;
}

出力

0 から 4294967295 の値をとる乱数が 4294966783 より大きい場合、条件Cが満たされる。
その確率は 1.19209e-07 である。

小数点以下を切り上げるという性質上、最高値のときの条件Bより、今回のほうが確率は大幅に低くなっている。

最後に、条件「D: 木曜午後のカブ価は9～65ベルの可能性があるが、これが9ベルである。」を満たす確率を考えてみる。
月曜午前時点で0.4である倍率は、月曜午後～木曜午後までの6期間において単調減少する。
この際、各期間において、倍率の下落幅が0.03～0.05の中から決められ、更新される。
90ベルが9ベルになるためには、倍率は $0.1 = 0.4 - 0.05 \times 6$ とならなければならず、6期間連続で下落幅がほぼ最大(0.05)になる場合を考えることになる。
しかし、カブ価の小数点以下切り上げ時に「0.99999を足してから小数を切り捨てる」という処理を行っているため、例えば9.0001ベルは10ベルに切り上げられるが、9.0000001ベルは9ベルとなる、という現象が発生する。そのため、下落幅は6回すべてが最大値でなくてもよい可能性がある。
それでも条件Cでは二分探索をすることで乱数の閾値を求められたが、それは乱数発生が1回だけであったからできたことであり、今回は6回の乱数が関係してくるため単純にはできない。
そこで、以下のような方法で確率を求める。

D-1. 乱数が最高値の $2^{32}-1$ である場合の倍率下落幅（0.03に最も近い（注））を $A$ とする。6回中5回の下落幅が $A$ であるという条件の下で、残り1回の下落幅がいくつ以上であったら9ベルになるか、という閾値を求める。その閾値を $B$ とする。
D-2. $B$ から $A$ の下落幅について、乱数がどの範囲であった場合その下落幅になるか、を求める。
D-3. 6回の下落幅それぞれを $B$ から $A$ の間で全探索し、9ベルになるものについては、そうなる場合の数（D-2で求めた幅を6つ掛け合わせたもの）を足し合わせる。
D-4. その値を $(2^{32})^6$ で割ったものが、条件Dを満たす確率である。

（注）0.05ではなく0.03であるのは、解析では「0.02を引いてから、0～0.03を引く」という操作をしており、それに倣ったため。

下落幅については10進数で小数点以下を表示するより、float型の仮数部を見たほうが扱いやすいので、そう表記する。
仮数部の16進数表現を得るのに、以下の記事を参考にした。
https://qiita.com/nia_tn1012/items/d26f0fc993895a09b30b
また、今までは10進数表記していた乱数の値も、これからは16進数表記することとする。

まずはD-1を二分探索で行う。

#include <bits/stdc++.h>
using namespace std;

// uint32_t型の乱数をもとに、float型の数を返す関数。乱数が大きいほどbに近い数となる。
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

// 小数点以下を切り上げる関数
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

// float型の仮数部を返す関数
int printb(float ff) {
    union { float f; int i; } a;
    int i;
    a.f = ff;

    return a.i & 0x7FFFFF;
  }

int main() {
  uint32_t ll = 0;  //uint32_tの下限値
  uint32_t rr = pow(2, 32) - 1;  //uint32_tの上限値
  int32_t basePrice = 90;  // 買値
  uint64_t l, r, m;

  float declineRate;
  int32_t thuPmPrice;
  l = 0;
  r = pow(2, 32) - 1;
  m = (l + r) / 2;

  // 二分探索
  while (r - l > 1) {
    declineRate = (0.4 - 0.03 * 6 - randfloat(0, 0.02, rr) * 5 - randfloat(0, 0.02, m));
    thuPmPrice = intceil(declineRate * basePrice);  // 木曜午後の値段
    if (thuPmPrice > 9) {l = m;}
    else {r = m;}
    m = (l + r) / 2;
  }
  cout << hex << "下落幅Bの仮数部は " << printb(randfloat(0, 0.02, m+1)) \
  << " であり、その下落幅になる中で最小の乱数の値は " << m+1 << " である。" << endl;
  cout << hex << "下落幅Aの仮数部は " << printb(randfloat(0, 0.02, rr)) << " である。" << endl;
}

出力

下落幅Bの仮数部は 23d6db であり、その下落幅になる中で最小の乱数の値は ffffb600 である。
下落幅Aの仮数部は 23d709 である。

よって、 23d6db ～ 23d709 を探索すればよいが、念のため 23d6da ～ 23d709 の48通りを調べる。（計算順序によっては、丸め誤差の関係で仮数部が 23d6da の下落幅でも9ベルになる可能性があるため）

次に、D-2からD-4を一気に行った。計算時間は10分ほどかかった。

#include <bits/stdc++.h>
using namespace std;

// uint32_t型の乱数をもとに、float型の数を返す関数。乱数が大きいほどbに近い数となる。
float randfloat(float a, float b, uint32_t u32)
  {
    uint32_t val = 0x3F800000 | (u32 >> 9);
    float fval = *(float *)(&val);
    return a + ((fval - 1.0f) * (b - a));
  }

// 小数点以下を切り上げる関数
int intceil(float val)
  {
    return (int)(val + 0.99999f);
  }

// float型の仮数部を返す関数
int printb(float ff) {
    union { float f; int i; } a;
    int i;
    a.f = ff;
    return a.i & 0x7FFFFF;
  }

// 倍率の下落のしかたを全探索する再帰関数
int64_t fullSearch(float declineRateNow, int step, int basePrice, vector<float> declineRates, int c) {
  int64_t ans = 0;
  int32_t thuPmPrice;
  if (step == 6) {
    thuPmPrice = intceil(declineRateNow * basePrice); // 木曜午後の値段が9ベルかを確認
    if (thuPmPrice == 9) {return 1;}
    else {return 0;}
  }
  else {
    declineRateNow -= 0.03;
    for (int i=0; i<c; i++) {
      ans += fullSearch(declineRateNow - declineRates.at(i), step + 1, basePrice, declineRates, c);
    }
  }
  return ans;
}

int main() {
  uint32_t m = pow(2, 32) - 1;
  float declineRate = randfloat(0, 0.02, m);
  float declineRate_b = randfloat(0, 0.02, m);
  uint32_t l, r;
  int c = 0;
  vector<float> declineRates(50);
  
  // D-2
  r = m;
  cout << "No. | 乱数下限 | 乱数上限 | 幅 | 下落幅の仮数部" << endl;
  while (true) {
    declineRate_b = declineRate;
    declineRate = randfloat(0, 0.02, m);
    if (printb(declineRate) < printb(declineRate_b)) {
      l = m + 1;
      cout << hex << setfill('0') << setw(2) << c << " | " << l << " | " << r << " | " \
      << r-l+1 << " | " << printb(declineRate_b) << endl;
      declineRates.at(c) = declineRate_b;
      r = m;
      c++;
      if (printb(declineRate) < 0x23d6da) {break;}
    }
    m--;
  }

  // D-3
  int32_t basePrice = 90;
  int64_t count9bell = fullSearch(0.4, 0, basePrice, declineRates, c);

  // D-4
  // (2^9)^6を掛けてから(2^32)^6で割る代わりに、2で23*6回割っている
  double prob_d = count9bell;
  for (int i=0; i<23*6; i++) {
    prob_d /= 2;
  }
  cout << dec << "37^6 通りある「倍率の下落のしかた」のうち、カブ価が9ベルとなるのは " << count9bell << " 通りである。" << endl;
  cout << "その場合、条件Dが満たされ、その確率は " << prob_d << " である。" << endl;
}

最初に、D-2として、仮数部が 23d6da ～ 23d709 の48通りである下落幅について、乱数がどの範囲であった場合その下落幅になるか、を求めた。以下の出力が得られた。

No. | 乱数下限 | 乱数上限 | 幅 | 下落幅の仮数部
00 | fffffe00 | ffffffff | 200 | 23d709
01 | fffffc00 | fffffdff | 200 | 23d707
02 | fffffa00 | fffffbff | 200 | 23d706
03 | fffff800 | fffff9ff | 200 | 23d705
04 | fffff600 | fffff7ff | 200 | 23d704
05 | fffff400 | fffff5ff | 200 | 23d702
06 | fffff200 | fffff3ff | 200 | 23d701
07 | fffff000 | fffff1ff | 200 | 23d700
08 | ffffee00 | ffffefff | 200 | 23d6fe
09 | ffffec00 | ffffedff | 200 | 23d6fd
0a | ffffea00 | ffffebff | 200 | 23d6fc
0b | ffffe800 | ffffe9ff | 200 | 23d6fb
0c | ffffe600 | ffffe7ff | 200 | 23d6f9
0d | ffffe400 | ffffe5ff | 200 | 23d6f8
0e | ffffe200 | ffffe3ff | 200 | 23d6f7
0f | ffffe000 | ffffe1ff | 200 | 23d6f6
10 | ffffde00 | ffffdfff | 200 | 23d6f4
11 | ffffdc00 | ffffddff | 200 | 23d6f3
12 | ffffda00 | ffffdbff | 200 | 23d6f2
13 | ffffd800 | ffffd9ff | 200 | 23d6f0
14 | ffffd600 | ffffd7ff | 200 | 23d6ef
15 | ffffd400 | ffffd5ff | 200 | 23d6ee
16 | ffffd200 | ffffd3ff | 200 | 23d6ed
17 | ffffd000 | ffffd1ff | 200 | 23d6eb
18 | ffffce00 | ffffcfff | 200 | 23d6ea
19 | ffffcc00 | ffffcdff | 200 | 23d6e9
1a | ffffca00 | ffffcbff | 200 | 23d6e7
1b | ffffc800 | ffffc9ff | 200 | 23d6e6
1c | ffffc600 | ffffc7ff | 200 | 23d6e5
1d | ffffc400 | ffffc5ff | 200 | 23d6e4
1e | ffffc200 | ffffc3ff | 200 | 23d6e2
1f | ffffc000 | ffffc1ff | 200 | 23d6e1
20 | ffffbe00 | ffffbfff | 200 | 23d6e0
21 | ffffbc00 | ffffbdff | 200 | 23d6de
22 | ffffba00 | ffffbbff | 200 | 23d6dd
23 | ffffb800 | ffffb9ff | 200 | 23d6dc
24 | ffffb600 | ffffb7ff | 200 | 23d6db

表の数値はすべて16進数である。
例えば、最後の行を見ると、（ 00000000 ～ ffffffff の値をとる乱数が） ffffb600 から ffffb7ff までの 0x200 通りの値をとった場合、下落幅の仮数部は 23d6db になる、ということである。
この表より、以下のことがわかる。

乱数がどんな値であっても、下落幅が 23d708, 23d703,... などの仮数部になることはない。（48通り調べているのに、表の行数が 0x25 すなわち37行なので、そのような下落幅が11通りあることがわかる。）
なりうる下落幅について、乱数の範囲の幅はすべて 0x200 すなわち $2^9 = 512$ である。（実はこれはrandfloat関数のコードからもわかる）

よって、D-3, D-4については、6回の下落幅それぞれをこの37通りで全探索し、下落のしかた $37^6$ パターンにおいてカブ価が9ベルとなった回数を数え、その回数に $(2^9)^6$ を掛けてから $(2^{32})^6$ で割ったものが、条件「D: 木曜午後のカブ価は9～65ベルの可能性があるが、これが9ベルである。」を満たす確率、ということになる。
これを計算した結果、以下の出力が得られた。

37^6 通りある「倍率の下落のしかた」のうち、カブ価が9ベルとなるのは 8689156 通りである。
その場合、条件Dが満たされ、その確率は 2.49367e-035 である。

変動パターン（跳ね小型）の確率と条件A～Dの確率すべてを掛け合わせると、
$0.2588 \times \frac{1}{21} \times \frac{1}{4} \times (1.192 \times 10^{-7}) \times (2.494 \times 10^{-35}) = 9.16 \times 10^{-45}$ となる。
これは天文学的な低確率である。
仮に、地球に住む77億人全員がこのゲームを持ち、138億年前の宇宙誕生から今に至るまで毎週カブ価を確認したとしても、1回でも「1カブ9ベル」を見られる確率は $5.08 \times 10^{-23}$ なのである。
なお、これは「前提」に記載したカブ価決定アルゴリズムの解析結果のソースコードが、完全にゲーム内部のものと同じである場合の確率である。小数点以下を切り上げるintceil関数で足す値が0.99999でなく0.999999であれば最安値の確率はもっと低くなるわけであり、またfloat型をdouble型にするだけでも確率は変わるだろう。

まとめ

カブ価が最高値の660ベルとなる確率は $2.68 \times 10^{-5}$ であり、毎週世界のどこかでは出ているといえる。
カブ価が最安値の9ベルとなる確率は $9.16 \times 10^{-45}$ であり、まず起こり得ない。

参考

あつまれどうぶつの森カブ価予測ツール | hyperWiki
カブの買値や売値を入力すると、カブ価の予測値や、カブ価が各変動パターンである確率が表示される。

AtCoder コードテスト
競技プログラミングのサイトAtCoderのコードテストページ。C++の標準出力に日本語を使用しても、文字化けされずに出力される。

更新履歴

2020/07/14 最安値の確率を修正、誤字修正（誤: 下落率 → 正: 下落幅）