HAKUTAI Tech Notes

IT関連、数学のことなどを主として思いつきで色々書き溜めていきます

【matter.js】曲がりくねった棒を簡単に生成する

matter.jsで曲がりくねった棒状のBodyを作ろうとすると思った以上に大変である。
Bodies.fromVertices()Svg.pathToVertices()で座標点を指定して作ると、凹面が含まれる形状の場合は凸包形(凹面が埋められる)に変換されてしまうためdecomp.jsを使う必要があり、あまり複雑な形状だと思い通りに機能してくれないことが多い。さらに、これらの関数ではドーナツ形のように空洞を含む形状はサポートされていない。
そこで、もっと単純な条件の入力だけで等幅の曲がりくねった棒状のBodyを生成できる関数を作成してみた。

棒の生成の考え方

関数的には棒の骨格となる中心線の座標データと棒の幅だけを引数として渡す。
やはりBodies.fromVertices()を使うが、棒全体の座標を与えて一気に作るのではなく、棒の節点(屈折部分)と節点の間に四角形を1つずつ生成してそれらを最後にくっつけて1本の棒にするという方針になる。
中心線の各節点から腕(棒の幅と曲がりの角度に応じて長さが変わる)を左右に伸ばして、棒の1部品となる四角形の頂点の座標を決める。 なお、ここで言う棒の左側・右側とは、中心線の進行方向を見た時の左側・右側という意味で使っている。

中心線の各線分の傾きの計算

中心線のn、n+1番目の節点を$S_n(x_n, y_n)$, $S_{n+1}(x_{n+1}, y_{n+1})$とした時、線分$S_nS_{n+1}$がx軸の正の向きとなす角$\theta_n$を計算する。
線分$S_nS_{n+1}$の傾きを$m$とすると、両節点のx, y座標の値から$\theta_n$は下記のように求められる。

$$ m = \frac{y_{n+1} - y_n}{x_{n+1} - x_n} = \tan \theta_n $$

$$ \theta_n = \tan^{-1}m $$

これと同等の計算は、javascriptのMath.atan2() 関数を利用すればいい。Math.atan2() は、原点 (0, 0) から任意の1点 (x, y) までの半直線とx軸の正の向きのなす角をラジアン単位で返す。原点以外の2点間で計算したい場合は、2点のx, y座標の差分をそれぞれ引数に指定すればいい。引数はy, xの順で渡す必要がある。

なお、数学的には角の向きは反時計回りを正とするが、matter.jsではy軸が下向きなので角の向きは時計回りが正となる点に注意する必要がある。

節点から左右に伸びる腕の長さの計算

節点$S_n$, $S_{n+1}$間に四角形(図の黄色い部分)を作る時、 四角形の両端の辺(図の$P_1P_2$, $P_3P_4$)が中心線の屈折角の二等分線となるようにしてやると、棒の幅に基づいて自然な感じで屈折部を作ることができる。
$\angle S_{n-1}S_nS_{n+1}$は線分$P_1P_2$で等分されるので、まずは右側の角$\alpha_n$を求める。

$$ \alpha_n = \frac{180^{\circ} - (\theta_n + \theta_{n-1})}{2} = 90^{\circ} - \frac{\theta_n + \theta_{n-1}}{2} $$ であるが、$\theta_n$と$\theta_{n-1}$の回転方向を考慮すると図の場合は$\theta_n \gt 0$, $\theta_{n-1} \lt 0$なので、

$$ 90^{\circ} - \frac{\theta_n - \theta_{n-1}}{2} $$ になる。屈折のパターンに応じて$\theta_{n-1}$, $\theta_n$の値が変わってくるが、いずれの場合もこの式で$\alpha_n$を計算することができる。


次に、節点から左右に伸ばす腕の長さ、つまり線分$S_nP_1$または$S_nP_2$の長さ$l_n$を求める。棒の右側の$P_1$から反対側へ下ろした垂線の足を$H_n$とし、線分$P_1H_n$と$P_1P_2$のなす角を$\beta_n$とすると、図形的な関係から$\beta_n$は次のように求められる。

$$ \beta_n = 90^{\circ} - \alpha_n = \frac{\theta_n - \theta_{n-1}}{2} $$

線分$P_1H_n$の長さ、つまり棒の幅を$t$とすると、これはあらかじめ引数で与えられるもので既知であるから、$l_n$は下記のように求められる。

$$ l_n = \frac{1}{2} \left| \frac{t}{\cos \beta_n} \right| $$

四角形を構成する4頂点の座標の計算

四角形の形状と位置を決定するために、4頂点$P_1$, $P_2$, $P_3$, $P_4$の座標を求める。色々やり方はあると思うが、ここでは座標点の回転移動を利用して求めることにする。なお、四角形の両端の辺について、線分$P_1P_2$を始点側の辺、線分$P_3P_4$を終点側の辺と呼ぶことにする。

まずは始点側の$P_1$, $P_2$の座標を計算する。 回転移動の対象とする点を$Q_s(x_s, y_s)$とする。$Q_s$は線分$S_nS_{n+1}$上に点$S_n$からの距離が$l_n$となるよう取れば計算しやすい。この時$Q_s$の座標は下記の計算で求められる。

$$ \begin{cases} x_s = x_n + l_n \cos \theta_n\\ y_s = y_n + l_n \sin \theta_n \end{cases} $$

後は、節点$S_n$周りに$Q_s$を$+\alpha_n$回転して右側の$P_1$の座標が求まり、$P_1$を180°回転して$P_2$が求まる。

終点側の$P_3$,$P_4$の座標の計算も同様に計算する。 回転移動の対象とする点を$Q_e(x_e, y_e)$として、線分$S_nS_{n+1}$上に点$S_{n+1}$からの距離が$l_{n+1}$となるように取る。この時$Q_e$の座標は下記の計算で求められる。

$$ \begin{cases} x_e = x_{n+1} - l_{n+1} \cos \theta_n\\ y_e = y_{n+1} - l_{n+1} \sin \theta_n \end{cases} $$

この後はやはり節点$S_{n+1}$周りに$Q_s$を回転して$P_3$,$ P_4$を求めるのだが、四角形を構成する都合上$P_1$→$P_2$→$P_3$→$P_4$の順で決定していかなければならないので先に左側の$P_3$を計算する。 $P_3$の座標は$Q_e$を節点$S_{n+1}$周りに$Q_s$を(180° − $\alpha_{n+1}$)回転し、$P_4$は$P_3$を180°回転すれば求まる。

座標の回転移動について

ここで、行列を用いた座標の回転移動について少しまとめておく。
一般に、座標の回転移動を表す行列は下記のように表す。

$$ \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \\ \end{pmatrix} $$

ただし、この行列は原点周りの回転を表したものであり、原点以外の任意の点周りの回転を行う場合は平行移動も組み合わせる必要がある。
画像処理などでは、点の平行移動:Translationと線形変換(回転:Rotation、拡大縮小:Scaling、せん断:Skew)をまとめてアフィン変換と呼び、3x3の行列で表すのが一般的である。せん断については今回は省くが、平行移動量$T_x$, $T_y$、回転角度$\theta$、拡大率$S_x$, $S_y$を表す行列は下記のように表される。

$$ \begin{pmatrix} S_x\cos \theta & -\sin \theta & T_x \\ \sin \theta & S_y\cos \theta & T_y \\ 0 & 0 & 1 \end{pmatrix} $$

特に今回は平行移動と回転移動しか考慮しないので$S_x$, $S_y$は1である。
任意の点$A(x_a, y_a)$の周りに点$B(x, y)$を$\theta$回転して点$B’(x’, y’)$に移す、つまり$\overrightarrow{AB}$を$\overrightarrow{AB’}$に変換することを考えるとき、次のように3段階の変換に分けて適用する必要がある。


① $\overrightarrow{AB}$の始点$A$を原点と重なる点$A_0$まで平行移動する(点$B$の移動先を点$B_0$とする)

$$ \begin{pmatrix} 1 & 0 & -x_a \\ 0 & 1 & -y_a \\ 0 & 0 & 1 \end{pmatrix} $$

② $\overrightarrow{A_0B_0}$を$\theta$回転する(点$B_0$の移動先を点$B_0'$とする)

$$ \begin{pmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

③$\overrightarrow{A_0B_0’}$の始点$A_0$を元の点$A$に戻るように平行移動する

$$ \begin{pmatrix} 1 & 0 & x_a \\ 0 & 1 & y_a \\ 0 & 0 & 1 \end{pmatrix} $$

これらの変換行列を点Bに適用(変換行列はオリジナルの座標に順次左側から掛け合わせていく)して点B’に変換することを式で表すと、

$$ \begin{pmatrix} x'\\ y'\\ 1 \end{pmatrix}= \begin{pmatrix} 1 & 0 & x_a \\ 0 & 1 & y_a \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & -x_a \\ 0 & 1 & -y_a \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x\\ y \\ 1 \end{pmatrix} $$

これを展開して整理すると、

$$ \begin{cases} x' = x\cos \theta - y\sin \theta - x_a\cos \theta + y_a\sin \theta + x_a \\ y' = x\sin \theta + y\cos \theta - x_a\sin \theta - y_a\cos \theta + y_a \end{cases} $$

四角形の中心座標の計算

4頂点$P_1$, $P_2$, $P_3$, $P_4$の配列をBodies.fromVertices()に渡せば棒の部品となる四角形を生成できる。しかし、4頂点の正確な中心座標を計算してfromVertices()の引数に指定しないと、本来生成したい位置からずれた所に四角形が生成されてしまいきれいな棒状にならない場合がある。
そこで、4頂点の中心座標を求めるために、一旦4頂点の座標値を元にMatter.Verticesオブジェクトを作る。このVerticesオブジェクトをVertices.centre()の引数として渡せば一発で中心座標を返してくれる。
ただし、Vertices.centre()で中心を計算するためには$P_1$→$P_2$→$P_3$→$P_4$と辿った時にぐるっと一周する位置関係でなければならず、クロスする位置関係になっていると正しく中心を計算することができない。

このようなクロスする位置関係になってしまうのは下記の2パターンである。

$$ \begin{cases} \frac{\theta_n - \theta_{n-1}}{2} \lt -90^{\circ} \rightarrow \alpha_n \lt 0^{\circ} \cdots (1)\\ \frac{\theta_n - \theta_{n-1}}{2} \gt 90^{\circ} \rightarrow \alpha_n \gt 180^{\circ} \cdots (2) \end{cases} $$

(1)(2)になるのは、下の図のように隣り合う角$\theta_n, \theta_{n-1}$が共に鈍角かつ異符号の時である。

(1)の場合は$\alpha_n$に180°を足し、(2)の場合は$\alpha_n$から180°を引くという補正を行えばクロスすることはなくなる。

ソースコード

ここまでの考え方に示した内容をコードに書き起こした。

/**
 * 中心線に沿った棒状のBodyを生成する
 */
function createMeandering(sections, thickness, shouldClose = false, options ={}) {
  const parts = [];
  const sectionAngles = [];

  // パスを閉じる場合、先頭の節点を末尾に追加する
  if (shouldClose) {
    sections.push(sections[0]);
  }

  for (let j = 0; j < sections.length; j++) {
    if (j !== sections.length - 1) {
      // 2頂点の線分とx軸のなす角度を計算する
      sectionAngles.push(Math.atan2(sections[j + 1].y - sections[j].y, sections[j + 1].x - sections[j].x));
    }
  }

  if (shouldClose) {
    const beginningAngle = sectionAngles[0];
    const endAngle = sectionAngles[sectionAngles.length - 1];
    sectionAngles.unshift(endAngle);
    sectionAngles.push(beginningAngle);
  } else {
    sectionAngles.unshift(sectionAngles[0]);
    sectionAngles.push(sectionAngles[sectionAngles.length - 1]);
  }

  for (let j = 0; j < sections.length; j++) {
    if (j !== sections.length - 1) {
      const vertices = [];

      // 始点側の2頂点を生成する
      const armAngle1 = (sectionAngles[j + 1] - sectionAngles[j]) / 2;
      const armLength1 = Math.abs(thickness / (2 * Math.cos(armAngle1)));
      const startRadius = { x: sections[j].x + armLength1 * Math.cos(sectionAngles[j + 1]), y: sections[j].y + armLength1 * Math.sin(sectionAngles[j + 1]) };
      // 始点側の回転角度を計算する
      const rotAngle1 = calcRotAngle(-armAngle1);
      // 始点側の動径ベクトルを回転して2頂点の座標を計算する
      vertices.push(rotate(sections[j], startRadius, rotAngle1));
      vertices.push(rotate(sections[j], vertices[0], Math.PI));

      // 終点側の2頂点を生成する
      const armAngle2 = (sectionAngles[j + 2] - sectionAngles[j + 1]) / 2;
      const armLength2 = Math.abs(thickness / (2 * Math.cos(armAngle2)));
      const endRadius = { x: sections[j + 1].x - armLength2 * Math.cos(sectionAngles[j + 1]), y: sections[j + 1].y - armLength2 * Math.sin(sectionAngles[j + 1]) };
      // 終点側の回転角度を計算する
      const rotAngle2 = calcRotAngle(armAngle2);
      // 終点側の動径ベクトルを回転して2頂点の座標を計算する
      vertices.push(rotate(sections[j + 1], endRadius, rotAngle2));
      vertices.push(rotate(sections[j + 1], vertices[2], Math.PI));

      // 4頂点で作られるVerticesオブジェクトの中心座標を求める
      const tetragon = Vertices.create(vertices, Matter.Body);
      const center = Vertices.centre(tetragon);
      // 四角形を生成してpartsに追加する
      parts.push(Bodies.fromVertices(center.x, center.y, vertices, options));
    }
  }

  // 結合前のパーツの配列を返却する
  return parts;
}

/**
 * 任意の点の周りに座標点を回転する
 */
function rotate(origin, target, rotAngle) {
  // 回転移動対象の点を、原点周りに平行移動 → 原点周りの移動 → 元の任意点周りに戻す
  const x = target.x * Math.cos(rotAngle) - target.y * Math.sin(rotAngle) + origin.x - origin.x * Math.cos(rotAngle) + origin.y * Math.sin(rotAngle);
  const y = target.x * Math.sin(rotAngle) + target.y * Math.cos(rotAngle) + origin.y - origin.x * Math.sin(rotAngle) - origin.y * Math.cos(rotAngle);

  return { x, y };
}

/**
 * 回転角度を計算する
 */
function calcRotAngle(armAngle) {
  let rotAngle = Math.PI / 2 + armAngle;
  // 回転角度が0°より小、180°より大になる場合の補正
  if (armAngle > Math.PI / 2) {
    rotAngle -= Math.PI;
  } else if (armAngle < -Math.PI / 2) {
    rotAngle += Math.PI;
  }

  return rotAngle;
}
createMeandeling()
中心線の座標データの配列 [object]
棒の幅 number
パスを閉じるフラグ boolean
Matter.Bodyのオプション object

メインの処理を行う関数。
指定が必須な中心線の座標データの配列と棒の幅の他に、パスを閉じるフラグとしてtrueを指定することで棒の始点と終点を自動的に結んで閉じたような形状を作れるようにした。
なお、中心線の屈折角が0°になる区間を含む場合はうまく生成できないので注意が必要。

なお、最終的に1本の棒状とする場合はcreateMeandering()から返るBodyの配列を別途結合してやる必要がある。個々の四角形を結合してから返すようにすると、そのBodyに対してさらに別のBodyを結合した時にやはり凸包形に変換されてしまうのでこのような返却の形にした。

rotate()
回転の中心点の座標 object
回転対象の点の座標 object
回転角度(ラジアン) number

任意の点の周りに座標点を回転する関数。

calcRotateAngle()
節点から伸びる腕の角度(ラジアン) number

長方形の4頂点を計算するために中心線上の点を回転する角度を計算する関数。
回転角度が0°より小、180°より大になる場合の補正もここで行う。

使用例

新潟県

const mainland = [
  { x: 466.9, y: 59.2 }, { x: 445.5, y: 105.8 }, { x: 439.7, y: 159.7 }, { x: 397.1, y: 206.4 }, { x: 324.3, y: 245.8 },
  { x: 283.6, y: 327.7 }, { x: 225.5, y: 392 }, { x: 181.2, y: 428 }, { x: 59.7, y: 477.1 }, { x: 77.8, y: 494.1 },
  { x: 87.8, y: 534.1 }, { x: 118.8, y: 498.8 }, { x: 142.1, y: 503 }, { x: 146, y: 522.7 }, { x: 199.3, y: 513.4 },
  { x: 225.5, y: 471.9 }, { x: 257.4, y: 466.7 }, { x: 267.1, y: 497.8 }, { x: 283.6, y: 510.3 }, { x: 286.5, y: 544.5 },
  { x: 316.6, y: 536.2 }, { x: 335, y: 517.5 }, { x: 349.6, y: 481.2 }, { x: 372.8, y: 460.5 }, { x: 392.2, y: 484.3 },
  { x: 405.8, y: 483.3 }, { x: 400, y: 422.1 }, { x: 390.3, y: 408.6 }, { x: 398.4, y: 356.5 }, { x: 435.8, y: 349.6 },
  { x: 476.6, y: 335 }, { x: 468.8, y: 301.8 }, { x: 507.6, y: 253.1 }, { x: 486.2, y: 232.3 }, { x: 499.8, y: 157.7 },
  { x: 527.9, y: 147.3 }, { x: 538.6, y: 127.6 }, { x: 506.6, y: 106.9 }, { x: 501.8, y: 71.6 }
];

const sado = [
  { x: 260.3, y: 184 }, { x: 235.9, y: 236.9 }, { x: 185.9, y: 259.2 }, { x: 207.9, y: 213.3 }, { x: 198.4, y: 205.7 },
  { x: 188.3, y: 214 }, { x: 192.5, y: 185.9 }, { x: 237.7, y: 122.9 }, { x: 255, y: 119.7 }, { x: 241.9, y: 159.8 },
  { x: 234.7, y: 187.2 }
];

const options = {
  render: {
    fillStyle: '#6b8e23'
  }
}

const mainlandParts = createMeandering(mainland, 10, true, options);
const sadoParts = createMeandering(sado, 10, true, options);
const niigata = Body.create({ parts: [...mainlandParts, ...sadoParts], isStatic: false });

Composite.add(engine.world, niigata);

// 床
Composite.add(engine.world, Bodies.rectangle(400, 585, 800, 30, { isStatic: true }));


図形の入れ子

// 正多角形の節点座標を計算する関数
function polygonSections(center, sides, radius) {
  const sections = [
    { x: center.x, y: center.y - radius }
  ];
  const rotAngle = 2 * Math.PI / sides;
  for (let i = 0; i < sides - 1; i++) {
    sections.push(rotate(center, sections[i], rotAngle));
  }

  return sections;
}

// 正二十角形
const icosagonSections = polygonSections({ x: 300, y: 200 }, 20, 200); 
// 正六角形
const hexagonSections = polygonSections({ x: 300, y: 200 }, 6, 150); 
// 星形
const starSections = [
  { x: 300, y: 100 },{ x: 358.78, y: 280.9 },{ x: 204.89, y: 169.1 },{ x: 395.11, y: 169.1 },{ x: 241.22, y: 280.9 }
];
// 床
const floorSections = [
  { x: 0, y: 400 },{ x: 100, y: 440 },{ x: 200, y: 470 },{ x: 300, y: 490 },{ x: 400, y: 500 },
  { x: 500, y: 490 },{ x: 600, y: 470 },{ x: 700, y: 440 },{ x: 800, y: 400 },
];

const options = {
  inertia: 1e7,
  density: 1,
}

const icosagon = createMeandering(icosagonSections, 30, true, { ...options, render: { fillStyle: '#ffd700' }});
const hexagon = createMeandering(hexagonSections, 15, true, { ...options, render: { fillStyle: '#0000ff' } });
const star = createMeandering(starSections, 10, true, { ...options, render: { fillStyle: '#ff0000' } });
const floor = createMeandering(floorSections, 30, false, { ...options, render: { fillStyle: '#4d4d4d' } });

Composite.add(engine.world, Body.create({ parts: icosagon }));
Composite.add(engine.world, Body.create({ parts: hexagon }));
Composite.add(engine.world, Body.create({ parts: star }));
Composite.add(engine.world, Body.create({ parts: floor, isStatic: true }));

確率計算だけでカタールW杯優勝国を予想する

先日、サッカーのカタールW杯に出場する日本代表のメンバーが発表され、いよいよ開催の気運が高まってきた。ということでW杯に絡めて何かやりたいなと思い、今大会でどこが優勝しそうかを予想するシミュレーションをやってみた。予想とは言っても真面目にチーム毎の戦力分析を行うとかいったものではなく、単に各チームの勝利確率を計算してそれだけを根拠に大会のレギュレーションに沿って対戦を進めていき、どこが優勝しそうかをシミュレーションしてみるというお遊びのものです。
シミュレーションの肝となる勝利確率の計算はイロレーティング(Elo Rating)の値をベースにした。



イロレーティングから勝利確率を計算する

イロレーティングとは

対戦型の競技において相対評価によって双方の実力を表すために使われる指標(レート)である。元々はチェスの強さを評価するために考案されたものだそうだ。
ある競技において参加プレイヤーの平均のレートを主に1500と定め、その値の大小により実力を評価することができる。試合が行われる度に2プレイヤー(チーム)のレート差と結果に応じて双方のレートが更新される。
サッカー競技におけるイロレーティングは、各国の国際Aマッチの結果を用いてその試合の重み(最高値がW杯決勝の60、最低値が親善試合の20)、ホームチームのアドバンテージ、得失点差を考慮して計算されている。いわゆる強豪国であればレートは2000を超えているところが多い。
各国の最新のレートは、「World Football Elo Ratings」というサービスで確認することができる。

www.eloratings.net

なお、サッカーの世界ランキングとしてFIFAランキングがあるが、2018年に計算方法が改訂されイロレーティングをベースに計算されるようになった。これにより、南米と欧州びいきの不公平な指標でしかなかった「大陸連盟間の強さ」は考慮されなくなった。

勝利確率の計算

対戦する2チームのレート差から勝利確率を計算することができる。

チームA、Bの勝利確率を$P_A$, $P_B$、レートを$R_A$, $R_B$とすると、

$$ P_A = \frac{1}{10^{\frac{R_B - R_A}{400}}+1} $$

で計算される。
チームBの勝利確率は余事象の確率から$P_B = 1 - P_A$である。

World Football Elo Ratings - About

引き分けについて

イロレーティングの対象は引き分けのない競技であることを前提としており、サッカーなどのように引き分けがあり得る競技の場合、引き分けは0.5勝0.5敗という扱いでレートが計算されている。このため、対戦両チームのレートから引き分けの確率だけを逆算して求めるといったことができない。したがって、なんらかの方法で引き分けの確率を別に求め(後述)、あらかじめレート差から計算しておいた両チームの勝利確率に対して引き分けを考慮した補正を行う必要がある。
引き分け確率を$P_d$として、$P_d$を考慮した時のチームA、Bの勝利確率$P_A'$, $P_B'$を下記のように計算した。(下の図は$P_A$ = 65%, $P_B$ = 35%, $P_d$ = 20%とした場合の例)

$$ \begin{align} P_A'& = P_A - \frac{P_d}{2}\\ P_B'& = P_B - \frac{P_d}{2} = 1-P_A - \frac{P_d}{2} \end{align} $$

なお、「$\frac{P_d}{2}$ ≧ 弱い方の勝利確率」となってしまう場合、「$P_d$ = 弱い方の勝利確率」とした。

シミュレーションの条件

実際のW杯における対戦の組み合わせに沿って、対戦国同士のレート差から算出した勝利確率・引き分け確率を用いて抽選を行い結果を決定していく。 本物の試合であれば実力差以外にも対戦の順序や選手の士気、コンディションなどが対戦結果に影響してくると思われるが、このシミュレーションでは対戦結果を左右するものはレート差から計算される確率だけである。また、他チームの対戦結果などの外的要因の影響は一切考慮しない独立なものとする。
また、途中の対戦結果を受けて各チームのレートが変動することは考えず、シミュレーション開始時点のレートを最後まで使う。

グループリーグの考え方

  • 勝利確率に基づき勝敗を決定する
  • 引き分けを考慮する
  • 勝:3、分:1、負:0の勝ち点を計上していき各グループの成績上位2チームが決勝トーナメントへ進む


実際のレギュレーションにおいて、グループリーグ内の順位を決定する項目の優先度は下記の通りである。

1. 勝点
2. 得失点差
3. 総得点
4. 直接対決の結果
5. フェアプレーポイント・反則ポイント
6. くじ引き

しかし、このシミュレーションはレート差しか見ておらず、チームごとの得失点などは考慮できないので下記の通りにしている。

1. 勝点
2. 直接対決の結果
3. くじ引き


なお、同グループ内で勝点が並ぶパターンは以下の3通りが考えられる。

(1) 2チームが同じ勝点で並ぶ
➡︎ 直接対決の結果で判断し、引き分けだったらくじ引き
(2) 3チームが同じ勝点で並ぶ
➡︎ 残りの1チームは全勝or全分けor全敗になっており、当該3チーム同士の直接対決の結果は全分けor三すくみ状態なのでくじ引き
(3) 全チーム同じ勝点になる
➡︎ リーグ戦6試合が全て引き分けになるか全チーム1勝1分1敗であり、直接対決の結果が同じなのでくじ引き

トーナメント戦の考え方

  • 勝率確率に基づき勝敗を決定する
  • 引き分けを考慮する
  • 引き分けの場合は勝敗50%ずつの抽選(PK戦のつもり)で決着を付ける

A〜Hの各グループリーグを1位、2位通過した16チームをトーナメントに割り振る方法も実際のレギュレーションに従う。


引き分け確率の求め方

前述の通りイロレーティングでは引き分けの確率を直接算出することができないので、なんらかの方法で別途求めなければならない。
そこで、実際に過去に行われたW杯10大会分(1982年スペイン大会〜2018年ロシア大会)のデータを利用し、対戦チームのレート差と引き分けの試合数の傾向から引き分け確率を求める式を推定した。幸いWorld Football Elo Ratingsには現在のレートだけでなく、過去に世界中で行われた国際親善試合や公式戦の結果も膨大に記録されているので非常にありがたい。
まず、過去10大会分の全試合のレート差の絶対値を幅50ずつの階級に分け、各階級毎の引き分けの割合(各階級の引き分けの試合数 / 各階級の総試合数)を算出した。さらに、各階級の中央値と引き分けの割合の分布から最小二乗法により近似直線の式を求め、レート差$\Delta R$から引き分け確率$R_d$を計算するための式とした。

$$ R_d = -0.00045\left| \Delta R \right| + 0.316 $$

なお、$450 < \left| \Delta R \right| \leq 500$の階級の引き分けの割合が飛び上がってしまっているが、この階級に含まれる試合数は他と比べて少ないことを考慮して近似式を計算する際には除外した。
「レート差が小さい、つまり力が拮抗しているほど引き分けになりやすい」という感覚的にもなんとなく良さそうだ。
ちなみに、この式は単純な一次関数なのでレート差が702以上の場合は引き分け確率が0以下になってしまうが、W杯出場国同士ではまずあり得ないレート差だと思われるのでとりあえずはヨシ。

シミュレーションの実行

各グループリーグの対戦〜トーナメントの決勝戦までを1回の試行として、全部で1000回繰り返した結果を集計する。
大会のシミュレーションを行うコードはPython 3.10.5で作成・実行した。
ソースコードは下記のGithubで公開している。
world_cup_winner_simulator


シミュレーションに用いた参加32ヶ国のレートの一覧は下記になる。括弧内はイロレーティングでの順位 である。※2022年11月1日時点

A 🇶🇦カタール 1664(50位)
🇪🇨エクアドル 1840(18位)
🇸🇳セネガル 1687(43位)
🇳🇱オランダ 2040(4位)
B 🏴󠁧󠁢󠁥󠁮󠁧󠁿イングランド 1920(14位)
🇮🇷イラン 1817(21位)
🇺🇸アメリカ 1798(24位)
🏴󠁧󠁢󠁷󠁬󠁳󠁿ウェールズ 1790(26位)
C 🇦🇷アルゼンチン 2141(2位)
🇸🇦サウジアラビア 1632(56位)
🇲🇽メキシコ 1813(22位)
🇵🇱ポーランド 1809(23位)
D 🇫🇷フランス 2005(6位)
🇦🇺オーストラリア 1719(39位)
🇩🇰デンマーク 1971(9位)
🇹🇳チュニジア 1687(43位)
E 🇪🇸スペイン 2045(3位)
🇨🇷コスタリカ 1736(35位)
🇩🇪ドイツ 1960(10位)
🇯🇵日本 1798(24位)
F 🇧🇪ベルギー 2025(5位)
🇨🇦カナダ 1770(29位)
🇲🇦モロッコ 1753(32位)
🇭🇷クロアチア 1922(13位)
G 🇧🇷ブラジル 2169(1位)
🇷🇸セルビア 1892(16位)
🇨🇭スイス 1929(12位)
🇨🇲カメルーン 1613(61位)
H 🇵🇹ポルトガル 2004(7位)
🇬🇭ガーナ 1540(74位)
🇺🇾ウルグアイ 1936(11位)
🇰🇷韓国 1783(28位)

結果

ローカルPC上でカタールW杯を1000回開催するのに要した時間は13.24秒だった。
なお、今回のシミュレーションで試行1回当たりの引き分け試合数の平均は14.7試合だった。実際の過去10大会における引き分け試合数の平均が1大会当たり14.6試合であることを考えるとかなり良い線を行っているのではないだろうか。

結果の集計方法

各国毎に最終到達ランクが「グループリーグ敗退、ベスト16、ベスト8、ベスト4、準優勝、優勝」になった回数をそれぞれカウントして度数分布で表す。
さらに、「グループリーグ敗退→0ポイント、ベスト16→1ポイント、ベスト8→2ポイント、ベスト4→3ポイント、準優勝→4ポイント、優勝→5ポイント」としてポイントの期待値(ポイント* 各ランクの到達回数 / 全試行回数)を計算する。期待値が5に近いほど今大会で好成績を残す可能性が高く、(期待値の意味から言えば厳密には異なるが)優勝に近い国であるといえる。

結果発表

では期待値の昇順で32位から発表していこう。


◆ 第32位 ◆

🇬🇭 ガーナ (期待値:0.030 / レート:1540)
リーグ 敗退 972
ベスト16 26
ベスト8 2
ベスト4 0
準優勝 0
優勝 0

◆ 第31位 ◆

🇨🇲 カメルーン (期待値:0.039 / レート:1613)
リーグ 敗退 966
ベスト16 30
ベスト8 3
ベスト4 1
準優勝 0
優勝 0

◆ 第30位 ◆

🇹🇳 チュニジア (期待値:0.119 / レート:1687)
リーグ 敗退 895
ベスト16 93
ベスト8 10
ベスト4 2
準優勝 0
優勝 0

◆ 第29位 ◆

🇸🇦 サウジアラビア (期待値:0.144 / レート:1632)
リーグ 敗退 880
ベスト16 100
ベスト8 16
ベスト4 4
準優勝 0
優勝 0

◆ 第28位 ◆

🇦🇺 オーストラリア (期待値:0.163 / レート:1719)
リーグ 敗退 858
ベスト16 124
ベスト8 15
ベスト4 3
準優勝 0
優勝 0

◆ 第27位 ◆

🇨🇷 コスタリカ (期待値:0.184 / レート:1736)
リーグ 敗退 860
ベスト16 108
ベスト8 22
ベスト4 8
準優勝 2
優勝 0

◆ 第26位 ◆

🇶🇦 カタール (期待値:0.218 / レート:1664)
リーグ 敗退 834
ベスト16 123
ベスト8 34
ベスト4 9
準優勝 0
優勝 0

◆ 第25位 ◆

🇲🇦 モロッコ (期待値:0.230 / レート:1753)
リーグ 敗退 818
ベスト16 142
ベスト8 33
ベスト4 6
準優勝 1
優勝 0

◆ 第24位 ◆

🇨🇦 カナダ (期待値:0.321 / レート:1770)
リーグ 敗退 757
ベスト16 179
ベスト8 55
ベスト4 5
準優勝 3
優勝 1

◆ 第23位 ◆

🇸🇳 セネガル (期待値:0.329 / レート:1687)
リーグ 敗退 765
ベスト16 156
ベスト8 67
ベスト4 10
準優勝 1
優勝 1

◆ 第22位 ◆

🇰🇷 韓国 (期待値:0.353 / レート:1783)
リーグ 敗退 714
ベスト16 236
ベスト8 36
ベスト4 11
準優勝 3
優勝 0

◆ 第21位 ◆

🇯🇵 日本 (期待値:0.354 / レート:1798)
リーグ 敗退 731
ベスト16 197
ベスト8 59
ベスト4 13
準優勝 0
優勝 0

◆ 第20位 ◆

🇺🇸 アメリカ (期待値:0.532 / レート:1798)
リーグ 敗退 631
ベスト16 247
ベスト8 91
ベスト4 23
準優勝 6
優勝 2

◆ 第19位 ◆

🏴󠁧󠁢󠁷󠁬󠁳󠁿 ウェールズ (期待値:0.580 / レート:1790)
リーグ 敗退 591
ベスト16 277
ベスト8 98
ベスト4 29
準優勝 5
優勝 0

◆ 第18位 ◆

🇵🇱 ポーランド (期待値:0.636 / レート:1809)
リーグ 敗退 566
ベスト16 305
ベスト8 72
ベスト4 43
準優勝 12
優勝 2

◆ 第17位 ◆

🇮🇷 イラン (期待値:0.670 / レート:1817)
リーグ 敗退 539
ベスト16 300
ベスト8 121
ベスト4 33
準優勝 6
優勝 1

◆ 第16位 ◆

🇲🇽 メキシコ (期待値:0.671 / レート:1813)
リーグ 敗退 532
ベスト16 334
ベスト8 81
ベスト4 41
準優勝 8
優勝 4

◆ 第15位 ◆

🇷🇸 セルビア (期待値:0.731 / レート:1892)
リーグ 敗退 570
ベスト16 246
ベスト8 108
ベスト4 44
準優勝 23
優勝 9

◆ 第14位 ◆

🇨🇭 スイス (期待値:1.028 / レート:1929)
リーグ 敗退 420
ベスト16 315
ベスト8 151
ベスト4 61
準優勝 37
優勝 16

◆ 第13位 ◆

🇪🇨 エクアドル (期待値:1.126 / レート:1840)
リーグ 敗退 356
ベスト16 306
ベスト8 240
ベスト4 60
準優勝 30
優勝 8

◆ 第12位 ◆

🇭🇷 クロアチア (期待値:1.181 / レート:1922)
リーグ 敗退 302
ベスト16 382
ベスト8 207
ベスト4 64
準優勝 32
優勝 13

◆ 第11位 ◆

🇩🇪 ドイツ (期待値:1.244 / レート:1960)
リーグ 敗退 304
ベスト16 359
ベスト8 208
ベスト4 67
準優勝 42
優勝 20

◆ 第10位 ◆

🇺🇾 ウルグアイ (期待値:1.293 / レート:1936)
リーグ 敗退 197
ベスト16 512
ベスト8 163
ベスト4 71
準優勝 43
優勝 14

◆ 第9位 ◆

🏴󠁧󠁢󠁥󠁮󠁧󠁿 イングランド (期待値:1.435 / レート:1920)
リーグ 敗退 239
ベスト16 352
ベスト8 236
ベスト4 104
準優勝 46
優勝 23

◆ 第8位 ◆

🇩🇰 デンマーク (期待値:1.704 / レート:1971)
リーグ 敗退 146
ベスト16 411
ベスト8 182
ベスト4 150
準優勝 76
優勝 35

◆ 第7位 ◆

🇵🇹 ポルトガル (期待値:1.724 / レート:2004)
リーグ 敗退 117
ベスト16 453
ベスト8 198
ベスト4 110
準優勝 65
優勝 57

◆ 第6位 ◆

🇧🇪 ベルギー (期待値:1.946 / レート:2025)
リーグ 敗退 123
ベスト16 338
ベスト8 240
ベスト4 131
準優勝 105
優勝 63

◆ 第5位 ◆

🇫🇷 フランス (期待値:1.960 / レート:2005)
リーグ 敗退 101
ベスト16 393
ベスト8 181
ベスト4 161
準優勝 98
優勝 66

◆ 第4位 ◆

🇪🇸 スペイン (期待値:2.018 / レート:2045)
リーグ 敗退 105
ベスト16 295
ベスト8 313
ベスト4 132
準優勝 74
優勝 81

◆ 第3位 ◆

🇳🇱 オランダ (期待値:2.244 / レート:2040)
リーグ 敗退 45
ベスト16 239
ベスト8 390
ベスト4 166
準優勝 73
優勝 87

◆ 第2位 ◆

🇦🇷 アルゼンチン (期待値:2.845 / レート:2141)
リーグ 敗退 22
ベスト16 240
ベスト8 166
ベスト4 245
準優勝 97
優勝 230

◆ 第1位 ◆

🇧🇷 ブラジル (期待値:2.948 / レート:2169)
リーグ 敗退 44
ベスト16 182
ベスト8 202
ベスト4 193
準優勝 112
優勝 267


ということで、優勝に最も近い国はブラジルでした!!
🎉おめでとうございます🥳

おわりに

成績の期待値で並べてみた結果、1〜8位までの上位集団については多少の前後はあるものの概ねイロレーティング 通りの並び順になった。期待値が飛び抜けている1位ブラジル、2位アルゼンチンには是非とも20年ぶりに欧州以外の国として優勝を掴み取ってもらいたい。
ところで、我が日本代表の期待値は残念ながら0.354(グループリーグ 敗退濃厚)という結果となり、改めて厳しいグループであることを認識させられてしまった…。
しかし、これを言ったらこの記事の意味を根底から覆してしまうようなものだが、数字だけで試合結果が決まってしまうのであればスポーツ観戦なんて面白くもなんともなくなってしまう。スポーツくじだって成立しないでしょう。
なので、日本代表には是非ともレーティングの順位や前評判を払拭するような大番狂わせを期待したい。


12/20追記
開幕からあっという間に日程が進んでいき、遂にアルゼンチンが36年ぶりの優勝に輝いてW杯の幕が閉じた。シミュレーションで優勝と予想したブラジルはベスト8で終わってしまったものの、見事に20年ぶりに南米勢がW杯を手にすることができました!!!
また、大変厳しいだろうという予想に反して日本がドイツとスペインに逆転勝ちを収めグループを1位通過した快進撃を筆頭に、確率の数字だけでは到底論じることができないような大波乱や熱戦が多く、夢のような1ヶ月でした。