大富豪の格差を計算してみる
日常の数学シリーズ。
大富豪（地域によっては大貧民とも呼ばれるトランプゲームの一種→Wikipedia:大富豪)で、カードの交換後に大富豪と大貧民で持っている２の枚数の期待値がどれくらいかを計算してみた*1。
- fugo.scm
; 人数 players: 4 ; カードの総数 all-cards: 52 ; プレイヤーの持つカードの数 player-cards: 13 ; 注目するカードの数 interesting-cards: 4 ; 移動するカードの数 exchange: 2 (define (comb n m) (cond ((or (< n 0) (< m 0)) (error "Negative numbers passed")) ((< n m) (error "First argument must not be smaller than second argument")) (#t (let loop ((n n) (m m) (r 1)) (cond ((< m 0) (error "Not an integer")) ((= m 0) r) (#t (loop (- n 1) (- m 1) (* r (/ n m))))))))) (define (probability-of-interesting-cards expected interesting-cards player-cards all-cards) (/ (* (comb interesting-cards expected) (comb (- all-cards interesting-cards) (- player-cards expected))) (comb all-cards player-cards))) (define standard-density (let ((p probability-of-interesting-cards)) (list (cons 0 (p 0 4 13 52)) (cons 1 (p 1 4 13 52)) (cons 2 (p 2 4 13 52)) (cons 3 (p 3 4 13 52)) (cons 4 (p 4 4 13 52))))) (define (filter-density carfn cdrfn density) (map (lambda (coord) (cons (carfn (car coord)) (cdrfn (cdr coord)))) density)) (define (filter-density-car f density) (filter-density f (lambda (x) x) density)) (define (filter-density-cdr f density) (filter-density (lambda (x) x) f density)) (define (tune-density density) (filter-density-cdr exact->inexact density)) (define (center seq) (let loop ((p1 seq) (p2 seq)) (if (or (null? p2) (null? (cdr p2))) p1 (loop (cdr p1) (cdr (cdr p2)))))) (define (copy-list seq) (apply list seq)) (define (merge seq1 seq2 cmpfn) (let loop ((seq1 seq1) (seq2 seq2) (r ())) (cond ((null? seq1) (append (reverse! r) seq2)) ((null? seq2) (append (reverse! r) seq1)) ((cmpfn (car seq1) (car seq2)) (loop (cdr seq1) seq2 (cons (car seq1) r))) (#t (loop seq1 (cdr seq2) (cons (car seq2) r)))))) (define (merge-sort! seq cmpfn) (let sort! ((seq seq)) (if (or (null? seq) (null? (cdr seq))) seq (let ((c (center (cons #f seq)))) (let ((seq2 (cdr c))) (set-cdr! c ()) (merge (sort! seq) (sort! seq2) cmpfn)))))) (define (merge-sort seq cmpfn) (merge-sort! (copy-list seq) cmpfn)) (define (normalize-density density) (define (sort-density density) (merge-sort density (lambda (x y) (< (car x) (car y))))) (let loop-1 ((rest (sort-density density)) (r ())) (if (null? rest) (reverse! r) (let loop-2 ((x (car rest)) (xs (cdr rest))) (if (null? xs) (loop-1 () (cons x r)) (let ((y (car xs)) (xs (cdr xs))) (if (= (car x) (car y)) (loop-2 (cons (car x) (+ (cdr x) (cdr y))) xs) (loop-1 (cons y xs) (cons x r))))))))) (define (shift-density density offset min max) (normalize-density (filter-density-car (lambda (x) (let ((x (+ x offset))) (cond ((< x min) min) ((> x max) max) (#t x)))) density))) (define density-looser (shift-density standard-density -2 0 4)) (define density-winner (let ((p probability-of-interesting-cards)) (normalize-density (list (cons 0 (* (p 0 4 13 52) (p 0 4 13 39))) (cons 1 (* (p 0 4 13 52) (p 1 4 13 39))) (cons 2 (* (p 0 4 13 52) (p 2 4 13 39))) (cons 2 (* (p 0 4 13 52) (p 3 4 13 39))) (cons 2 (* (p 0 4 13 52) (p 4 4 13 39))) (cons 1 (* (p 1 4 13 52) (p 0 3 13 39))) (cons 2 (* (p 1 4 13 52) (p 1 3 13 39))) (cons 3 (* (p 1 4 13 52) (p 2 3 13 39))) (cons 3 (* (p 1 4 13 52) (p 3 3 13 39))) (cons 2 (* (p 2 4 13 52) (p 0 2 13 39))) (cons 3 (* (p 2 4 13 52) (p 1 2 13 39))) (cons 4 (* (p 2 4 13 52) (p 2 2 13 39))) (cons 3 (* (p 3 4 13 52) (p 0 1 13 39))) (cons 4 (* (p 3 4 13 52) (p 1 1 13 39))) (cons 4 (* (p 4 4 13 52) (p 0 0 13 39))))))) (display (tune-density density-winner)) (newline) (display (tune-density density-looser)) (newline) ; (tune-density density-winner) ; ((0 . 0.0552220888355342) ; (1 . 0.249699879951981) ; (2 . 0.420264105642257) ; (3 . 0.235966386554622) ; (4 . 0.0388475390156062)) ; (tune-density density-looser) ; ((0 . 0.956158463385354) ; (1 . 0.0412004801920768) ; (2 . 0.00264105642256903))
というわけで、参加人数4人の場合、大富豪は2を2枚持っている確率が42%、1枚も持っていない確率は5%しかない。対して、大貧民の方は95%の確率で２が1枚もない状態からスタートする。分かっていたことだけど、やっぱり貧富の差は大きい。
同じことをPythonで書いた場合（Python Windows CE port 使用）
- fugo.py
# -*- coding: utf-8 -*- from itertools import groupby from operator import itemgetter from rational import Rational # 人数 players: 4 # カードの総数 all-cards: 52 # プレイヤーの持つカードの数 player-cards: 13 # 注目するカードの数 interesting-cards: 4 # 移動するカードの数 exchange: 2 def comb(n, m): if n < 0 or m < 0: raise "Negative numbers passed" elif n < m: raise "First argument must not be smaller than second argument" else: r = Rational(1) while m > 0: r = r * Rational(n, m) n -= 1 m -= 1 if m < 0: raise "Not an integer" return r def probability_of_interesting_cards(expected, interesting_cards, player_cards, all_cards): return comb(interesting_cards, expected) * \ comb(all_cards - interesting_cards, player_cards - expected) / \ comb(all_cards, player_cards) def standard_density(): p = probability_of_interesting_cards return [(0, p(0, 4, 13, 52)), (1, p(1, 4, 13, 52)), (2, p(2, 4, 13, 52)), (3, p(3, 4, 13, 52)), (4, p(4, 4, 13, 52))] standard_density = standard_density() def filter_density(carfn, cdrfn, density): return [(carfn(x), cdrfn(xs)) for x, xs in density] def filter_density_car(f, density): return filter_density(f, lambda x: x, density) def filter_density_cdr(f, density): return filter_density(lambda x: x, f, density) def tune_density(density): return filter_density_cdr(float, density) def normalize_density(density): def sort_density(density): return sorted(density, key=itemgetter(0)) r = [] for key, it in groupby(sort_density(density), itemgetter(0)): values = map(itemgetter(1), it) r.append((key, reduce(lambda x, y: x + y, values, Rational(0)))) return r def shift_density(density, offset, min, max): def f(x): x += offset if x < min: return min elif x > max: return max else: return x density = filter_density_car(f, density) return normalize_density(density) def density_looser(): density = shift_density(standard_density, -2, 0, 4) return normalize_density(density) density_looser = density_looser() def density_winner(): p = probability_of_interesting_cards density = [ (0, p(0, 4, 13, 52) * p(0, 4, 13, 39)), (1, p(0, 4, 13, 52) * p(1, 4, 13, 39)), (2, p(0, 4, 13, 52) * p(2, 4, 13, 39)), (2, p(0, 4, 13, 52) * p(3, 4, 13, 39)), (2, p(0, 4, 13, 52) * p(4, 4, 13, 39)), (1, p(1, 4, 13, 52) * p(0, 3, 13, 39)), (2, p(1, 4, 13, 52) * p(1, 3, 13, 39)), (3, p(1, 4, 13, 52) * p(2, 3, 13, 39)), (3, p(1, 4, 13, 52) * p(3, 3, 13, 39)), (2, p(2, 4, 13, 52) * p(0, 2, 13, 39)), (3, p(2, 4, 13, 52) * p(1, 2, 13, 39)), (4, p(2, 4, 13, 52) * p(2, 2, 13, 39)), (3, p(3, 4, 13, 52) * p(0, 1, 13, 39)), (4, p(3, 4, 13, 52) * p(1, 1, 13, 39)), (4, p(4, 4, 13, 52) * p(0, 0, 13, 39)), ] return normalize_density(density) density_winner = density_winner() if __name__ == '__main__': print tune_density(density_winner) print tune_density(density_looser) # tune_density(density_winner) # [(0, 0.055222088835534214), # (1, 0.24969987995198079), # (2, 0.42026410564225691), # (3, 0.23596638655462185), # (4, 0.038847539015606246)] # tune_density(density_looser) # [(0, 0.95615846338535415), # (1, 0.041200480192076831), # (2, 0.0026410564225690276)]
rational.pyは小町算 #3 - odz bufferのものを一部修正して使わせてもらった。
- rational.py
def gcd(x, y): while y != 0: x, y = y, x % y return x def lcm(x, y): return x * y / gcd(x, y) class Rational(object): def __init__(self, numerator, denominator=1): self.numerator = int(numerator) self.denominator = int(denominator) assert denominator != 0 c = gcd(self.numerator, self.denominator) if c != 1: self.numerator /= c self.denominator /= c if self.denominator < 0: self.numerator = -self.numerator self.denominator = -self.denominator def __neg__(self): return Rational(-self.numerator, self.denominator) def __add__(self, other): denominator = lcm(self.denominator, other.denominator) a = denominator / self.denominator b = denominator / other.denominator return Rational(self.numerator * a + other.numerator * b, denominator) def __sub__(self, other): return self + (-other) def __mul__(self, other): return Rational(self.numerator * other.numerator, self.denominator * other.denominator) def __div__(self, other): return self * Rational(other.denominator, other.numerator) def __float__(self): return float(self.numerator) / self.denominator def __int__(self): return self.numerator / self.denominator def __eq__(self, other): if isinstance(other, Rational): return self.numerator == other.numerator and \ self.denominator == other.denominator elif isinstance(other, int): return self.denominator == 1 and self.numerator == other else: return float(self) == other def __str__(self): if self.denominator == 1: return str(self.numerator) else: return '%d/%d' % (self.numerator, self.denominator) def __repr__(self): return '%d/%d' % (self.numerator, self.denominator)
ちなみに、ここまで全部Ad[es]上でやってます。
