range.rb

# -*- coding: utf-8 -*-
=begin
  This script was written by Takashi SUGA on July 2017 - December 2022
  You may use and/or modify this file according to the license described in the MIT LICENSE.txt file
  https://raw.githubusercontent.com/suchowan/watson-api-client/master/LICENSE.

  See also https://suchowan.seesaa.net/article/202212article_25.html
=end

require 'pp'
require 'when_exe'
require 'combination_generator'

# map メソッド
class CombinationGenerator
  def to_a
    combinations = []
    self.each do |combination|
      combinations << combination
    end
    combinations
  end
end

# 藤井聡太 対戦相手のレートと勝敗 (2022年12月24日時点)
# レートは http://kishibetsu.com/rating.html による
Results = %w(
  1236○ 1436○ 1344○ 1340○ 1661○ 1544○ 1417○ 1557○ 1336○ 1549○
  1594○ 1526○ 1792○ 1503○ 1508○ 1441○ 1591○ 1536○ 1705○ 1772○
  1608○ 1505○ 1611○ 1585○ 1607○ 1572○ 1419○ 1761○ 1712○ 1779●
  1469○ 1617○ 1701● 1279○ 1489○ 1477○ 1856● 1666○ 1383○ 1649○
  1540○ 1885● 1540● 1688○ 1651● 1495○ 1376○ 1464○ 1466● 1535○
  1576○ 1650○ 1514○ 1472○ 1634○ 1396○ 1483○ 1670○ 1641● 1563○
  1798● 1689○ 1735○ 1742● 1476○ 1651● 1469○ 1747○ 1819○ 1493○
  1516○ 1543○ 1446○ 1592○ 1509○ 1829○ 1821○ 1616○ 1551○ 1546○
  1737○ 1796○ 1542● 1534○ 1711○ 1628○ 1678○ 1659○ 1478○ 1588○
  1604● 1777○ 1722○ 1764● 1474○ 1806● 1436○ 1528○ 1707○ 1627○

  1305○ 1222○ 1405○ 1728○ 1822● 1790● 1584○ 1680○ 1655○ 1553○
  1550○ 1677○ 1606○ 1603○ 1750○ 1825● 1542○ 1668○ 1481○ 1648○
  1565○ 1382○ 1787○ 1806○ 1456○ 1516○ 1736● 1518○ 1728○ 1926○
  1747○ 1604○ 1704○ 1462○ 1766● 1535○ 1628○ 1704○ 1686○ 1727○
  1633○ 1595○ 1908● 1688● 1838○ 1783● 1230○ 1496○ 1524○ 1830○
  1527○ 1789○ 1334○ 1804○ 1807● 1798● 1742○ 1937● 1544○ 1767●
  1620○ 1661○ 1716● 1711○ 1500○ 1604○ 1777○ 1904● 1520○ 1761○
  1881○ 1583○ 1668○ 1557○ 1796○ 1897● 1637○ 1592○ 1587○ 1816○
  1607○ 1529○ 1543○ 1847○ 1807○ 1573○ 1664○ 1589○ 1512○ 1843●
  1845○ 1825○ 1505○ 1604○ 1740○ 1607● 1488○ 1590○ 1556○ 1770○

  1859○ 1830○ 1855○ 1770○ 1938○ 1950○ 1728● 1650○ 1622○ 1932○
  1724○ 1933○ 1821○ 1740○ 1850○ 1655○ 1935● 1819○ 1952○ 1827○
  1518○ 1746● 1649○ 1804○ 1728○ 1812○ 1637○ 1906● 1949○ 1831●
  1791○ 1927● 1798● 1780○ 1667○ 1963● 1764○ 1799○ 1585○ 1761○
  1441○ 1596○ 1603○ 1550○ 1666○ 1957○ 1653○ 1504○ 1948○ 1751○
  1828○ 1640○ 1672○ 1828○ 1729○ 1709● 1730○ 1699○ 1764○ 1921○
  1744● 1952○ 1675○ 1951○ 1733○ 1856○ 1943● 1945○ 1742○ 1804○
  1942○ 1936○ 1930○ 1682○ 1559○ 1925● 1775○ 1936○ 1933○ 1770○
  1933○ 1928● 1939○ 1921○ 1673○ 1762○ 1939○ 1937● 1861● 1812○
  1780○ 1818○ 1696○ 1820○ 1938○ 1681● 1714○ 1937○ 1939○ 1901○

  1927○ 1790○ 1922○ 1681○ 1797○ 1911● 1879● 1805○ 1925○ 1764● 
  1683○ 1912● 1909○ 1904○ 1675○ 1898○ 1710○ 1787○ 1744● 1779○ 
  1764○ 1935● 1953○ 1719○ 1878● 1956○ 1556○ 1895○ 1939○ 1890○ 
  1716○ 1787○ 1794● 1913○ 1560○ 1740○ 1918○ 1811○ 1910○ 1783○ 
  1902○ 1813● 1800○ 1896○ 1834○ 1831○ 1724○ 1801● 1806○ 1832○ 
  1829○ 1815○ 1814○ 1828● 1765○ 1841○ 1806○ 1814○ 1811○)

class Rating

  # 負け数マージン
  LossMargin = 0

  # 簡略計算閾値
  LossLimit = 7

  # レートの種
  RateSeed = 2000

  # 勝率キャッシュ
  Rates = {}

  # 組み合わせ
  Combinations = []

  def initialize(period)

    # 期間と、その期間内の結果
    @period = period
    results =Results[@period]

    # 対戦相手のレート
    @opponents = results.map(&:to_f)

    # 平均
    @average = @opponents.inject(&:+) / @opponents.size

    # メジアン
    d, m = (@opponents.size - 1).divmod(2)
    sort = @opponents.sort
    @median = (sort[d] + sort[d+m]) / 2.0

    # 負け数
    @total_loss = results.select {|result| result =~ /●/ }.length
    return if @total_loss > LossLimit

    # 組み合わせ
    laptime0 = Time.now.to_f
    Combinations[@total_loss] ||= (0..(@total_loss+LossMargin)).to_a.map {|loss|
      CombinationGenerator.new(loss, (0...@opponents.length).to_a).to_a
    }
    p [@period, 'laptime0-1', Time.now.to_f - laptime0]
  end

  def estimate
    # 点推定(平衡レート)
    laptime1 = Time.now.to_f
    p [@period, @total_loss, @average, 'point_estimation', estimate_point.round(2)]
    laptime2 = Time.now.to_f
    p [@period, 'laptime1-2', laptime2 - laptime1]
    results = [@total_loss, @average, @median, estimate_point]
    return results if @total_loss > LossLimit

    # 危険率 2.5% の場合の区間推定
    [ 0.025 ].each do |threshold|
      results += estimate_range(threshold)
      p [@period, threshold, results[-2..-1].map {|rate| rate.round(2)}] #=>
    end
    laptime3 = Time.now.to_f
    p [@period, 'laptime2-3', laptime3 - laptime2]
    return results

  rescue => e
    pp e.message
    pp e.backtrace
    return results

=begin
    # レートと負け数の関係
    losses = {}
    [1600,1700,1800,1900,2000,2100,2200].each do |rate|
      losses[rate] = estimate_loss(rate)
    end
    (1..(@total_loss+LossMargin)).each do |loss|
      losses[When::Ephemeris.root(1800, loss) {|rate| estimate_loss(rate)}] = loss.to_f
    end
    losses.keys.sort.each do |rate|
      p [@period, rate.round(2), losses[rate].round(2)]
    end
    laptime4 = Time.now.to_f
    p [@period, 'laptime3-4', laptime4 - laptime3]

    # 最尤値での負け数ごとの確率
    p (0..@total_loss+LossMargin).to_a.map {|loss|
      total_probability(estimate_point, loss).round(4)
    }
    p [@period, 'laptime4-5', Time.now.to_f - laptime4]
=end
  end

  # 点推定
  def estimate_point
    @estimated ||= When::Ephemeris.root(RateSeed, 0) {|rate| change_rate(rate)}
  end

  # 期待負け数
  def estimate_loss(rate)
    (0...@opponents.length).to_a.map {|i|
      game_probability(@opponents[i]-rate)
    }.inject(:+)
  end

  # レーティング法によるレートの変化
  def change_rate(rate)
    estimate_loss(rate) - @total_loss
  end

  # 区間推定
  def estimate_range(threshold=0.05)
    [[threshold,@total_loss],[1-threshold,@total_loss-1]].map {|th,loss|
      When::Ephemeris.root(@estimated || RateSeed, th) {|rate| range_probability(rate, loss)}
    }
  end

  # 指定の負け数以下の発生確率
  def range_probability(rate, loss_limit)
    (0..loss_limit).to_a.map {|loss|
      total_probability(rate, loss)
    }.sort.inject(:+)
  end

  # 指定の負け数の発生確率
  def total_probability(rate, loss)
    Combinations[@total_loss][loss].map {|losses|
      series_probability(rate, losses)
    }.sort.inject(:+)
  end

  # 指定の勝敗組み合わせの発生確率
  def series_probability(rate, losses)
    (0...@opponents.length).to_a.map {|i|
      r = rate - @opponents[i]
      game_probability(losses.include?(i) ? -r : r)
    }.inject(:*)
  end

  # 一局の勝率
  def game_probability(r)
    Rates[r] ||= 1 / (1 + 10.0 ** (-r / 400.0))
  end
end

open('fujii-rate.csv', 'w') do |csv|
(Results.length-29).times do |start|
    period = start...(start+30)
    rating = Rating.new(period)
    csv.puts((["#{start+1}..#{start+30}"] + rating.estimate.map {|item| item.to_s}).join(','))
  end
end