確率計算

過去に書いたプログラムが、投稿したwebサービスから消えていたのでブログに残しておく。(2011年11月に作成)

ブラックシティとホワイトフォレストの最初に決まる13人の中にグレースがいる確率を求めるプログラム (だったと思う)。だだぢぢさんがTwitterで話題にしていて触発されて書いた。
各住人の重みのデータはktxadさんが公開されたものによる (URLは紛失)。
(追記2014/10/15: ktxadさんが公開されていたソースコードは現在消えていますが、dorara32さんがかわりに重みのデータを公開されています: http://note.chiebukuro.yahoo.co.jp/detail/n166836)

# encoding: utf-8
require "mathn" # 整数同士で割り算したら分数になるように

# 重みのついたデータから重複しないように何件かランダムに選ぶとき、
# 特定の1件が選ばれる確率を計算する

Candidate = Struct.new(:id, :weight)

CANDIDATES = [
  [ 0, 20], [ 1, 20], [ 2, 20], [ 3, 20], [ 4, 20], [ 5, 20], [ 6, 20],
  [ 7, 10], [ 8, 20], [ 9, 10], [10, 20], [11,  6], [12, 20], [13,  4],
  [14,  1], [15, 20], [16, 20], [17, 20], [18, 20], [19, 20], [20, 20],
  [21, 10], [22, 20], [23, 20], [24, 10], [25, 20], [26,  6], [27,  4],
  [28, 20], [29,  1]
].map{|(id,weight)| Candidate.new(id, weight) }

TARGET_ID = 29

CALCLATORS = [
  [:MonteCarlo, "モンテカルロ法による確率"],
  [:Naive,      "ナイーブな方法による確率"],
  [:Fast,       "高速化した方法による確率"],
]

def main
  [1, 2, 3, 10, 13].each do |num_choice|
    puts "選ぶ個数: #{num_choice}"
    CALCLATORS.each do |(class_name, description)|
      klass = ProbabilityCalclator.const_get(class_name)
      if class_name == :Naive and num_choice >= 5
        pr = "(時間がかかるため省略)"
      else
        pr = klass.new.calc(num_choice)
      end
      print "  #{description}: #{pr}"
      if pr.kind_of?(Rational)
        print " ≒ #{pr.to_f}"
      end
      puts
    end
  end
end

def sumup_weight(candidates)
  candidates.map {|x| x.weight }.inject(:+)
end

ProbabilityCalclator = Module.new

# モンテカルロ法
class ProbabilityCalclator::MonteCarlo
  TRY_TIMES = 10000

  def calc(num_choice)
    succeeded_count = 0
    TRY_TIMES.times do
      succeeded_count += 1 if try(num_choice)
    end
    succeeded_count.to_f / TRY_TIMES
  end

  def try(num_choice)
    candidates = CANDIDATES.dup
    num_choice.times do
      index = random_choose(candidates)
      candidate = candidates.delete_at(index)
      if candidate.id == TARGET_ID
        return true
      end
    end
    false
  end

  def random_choose(candidates)
    random_number = rand(sumup_weight(candidates))
    sum = 0
    candidates.each_with_index do |candidate, index|
      sum += candidate.weight
      if sum > random_number
        return index
      end
    end
    raise "unreachable"
  end
end

# ナイーブな方法
# CANDIDATES.size = 30, num_choice = 13 の場合
# 多分 29P13 回ぐらい探索しないといけなくて遅い
class ProbabilityCalclator::Naive
  def calc(num_choice)
    @visited = [false] * CANDIDATES.size
    @total_sum = sumup_weight(CANDIDATES)
    recur(num_choice)
  end

  # 残りnum個での確率を計算
  def recur(num)
    return 0 if num == 0
    result = 0
    CANDIDATES.each_with_index do |candidate, index|
      next if @visited[index] # 今すでに選択している候補はスキップ
      if candidate.id == TARGET_ID
        # TARGET_IDだった場合、確率は1。残りの分の探索はやらなくてOK
        pr = 1
      else
        # candidateを選び、まだ決定していない残りnum - 1個分での確率を計算する
        @visited[index] = true
        @total_sum -= candidate.weight
        pr = recur(num - 1)
        @total_sum += candidate.weight
        @visited[index] = false
      end
      result += (candidate.weight / @total_sum) * pr
    end
    result
  end
end

# 同じ重みの値をまとめて扱い、さらにメモ化することで高速化
class ProbabilityCalclator::Fast
  def calc(num_choice)
    # 重みの値をキーにしてその重みの候補があと何個残っているかを記憶する
    @weight_to_count = CANDIDATES.each_with_object(Hash.new(0)) {|candidate, r|
      r[candidate.weight] += 1
    }
    @total_sum = sumup_weight(CANDIDATES)
    @memo = Hash.new
    @weights = @weight_to_count.keys
    @target_weight =  CANDIDATES.find{|x| x.id == TARGET_ID }.weight
    recur(num_choice)
  end

  def recur(times)
    # メモ化
    # (@weight_to_countが同じ内容ならばrecur_internalの結果も同じだから
    #  記憶しておいて使いまわす)
    @memo[@weight_to_count.dup] ||= recur_internal(times)
  end

  def recur_internal(times)
    return 0 if times == 0
    result = 0
    @weights.each do |weight|
      count = @weight_to_count[weight] # 重みがweightの値であり未選択の候補の数

      if weight == @target_weight
        pr = 1 / count # TARGET_ID を選んだときの確率
        if count > 1
          # TARGET_IDとweightは同じだがTARGET_IDではないものを選んだときの確率
          pr += ((count - 1) / count) * choose(times, weight)
        end
      else
        pr = choose(times, weight)
      end
      result += pr * ((count * weight) / @total_sum)
    end
    result
  end

  # weightを選択した場合の残りtimes - 1個分での確率を計算
  def choose(times, weight)
    if @weight_to_count[weight] > 0
      @weight_to_count[weight] -= 1
      @total_sum -= weight
      pr = recur(times - 1)
      @total_sum += weight
      @weight_to_count[weight] += 1
      pr
    else
      0
    end
  end
end

if $0 == __FILE__
  main()
end

実行結果は次の通り。
なお「モンテカルロ法による確率」は不正確な結果で、「ナイーブな方法による確率」と「高速化した方法による確率」は完全に正確な結果である。

選ぶ個数: 1
  モンテカルロ法による確率: 0.0025
  ナイーブな方法による確率: 1/462 ≒ 0.0021645021645021645
  高速化した方法による確率: 1/462 ≒ 0.0021645021645021645
選ぶ個数: 2
  モンテカルロ法による確率: 0.0035
  ナイーブな方法による確率: 9729246465/2204009196532 ≒ 0.004414340230661892
  高速化した方法による確率: 9729246465/2204009196532 ≒ 0.004414340230661892
選ぶ個数: 3
  モンテカルロ法による確率: 0.0064
  ナイーブな方法による確率: 1568664824856411452425433501221/232169860937309598063431523052200 ≒ 0.006756539451431991
  高速化した方法による確率: 1568664824856411452425433501221/232169860937309598063431523052200 ≒ 0.006756539451431991
選ぶ個数: 10
  モンテカルロ法による確率: 0.0278
  ナイーブな方法による確率: (時間がかかるため省略)
  高速化した方法による確率: 65770250293419026078129194848669982508897034077533201092414085189931935552537702182324234499375672068027368327207661273376859495758981364818701358782794795202517/2475647499344346203202200274227924547084627755009408713100317804864210058755132452548664905544294615343129859802504836359196285761348970446409580976203759822039600 ≒ 0.026566888182117086
選ぶ個数: 13
  モンテカルロ法による確率: 0.0396
  ナイーブな方法による確率: (時間がかかるため省略)
  高速化した方法による確率: 72920965323202474895401221315730728325203050221435766480570752583389314598872312680973006197795827752084921849225068647216327771173023300970901986908137886930833925903506944736597859681/1932080335186385797665285107673085006540033135470821753828553162094595243823360809973187589984927249351524831714776478574108021855863498572634867939515180025821224902545585179805950634000 ≒ 0.03774220149917724

C++による実装も実は作っていたのでした。

#include <iostream>
#include <vector>
#include <map>

// 重みのついたデータから重複しないように何件かランダムに選ぶとき、
// 特定の1件が選ばれる確率を計算する

struct Candidate {
	int id;
	int weight;
};

const Candidate CANDIDATES[] = {
	{ 0, 20}, { 1, 20}, { 2, 20}, { 3, 20}, { 4, 20}, { 5, 20}, { 6, 20},
	{ 7, 10}, { 8, 20}, { 9, 10}, {10, 20}, {11,  6}, {12, 20}, {13,  4},
	{14,  1}, {15, 20}, {16, 20}, {17, 20}, {18, 20}, {19, 20}, {20, 20},
	{21, 10}, {22, 20}, {23, 20}, {24, 10}, {25, 20}, {26,  6}, {27,  4},
	{28, 20}, {29,  1}
};

const int NUM_CANDIDATES = sizeof(CANDIDATES) / sizeof(Candidate);

const int TARGET_ID = 29;

double calc(int num_choice);
void initialize();
double recur(int times);
double recur_internal(int times);
double choose(int times, int weight);

int main(void) {
	for (int num_choice = 1; num_choice <= NUM_CANDIDATES; num_choice ++) {
		std::cout << "num_choice = " << num_choice << ": " << calc(num_choice) << std::endl;
	}
}

typedef std::map<int,int> WeightToCount;

// 重みの値をキーにしてその重みの候補があと何個残っているかを記憶する
WeightToCount g_weight_to_count;
std::map<WeightToCount, double> g_memo;
int g_total_sum;

int g_target_weight;
std::vector<int> g_weights;

double calc(int num_choice) {
	initialize();
	return recur(num_choice);
}

void initialize() {
	g_weight_to_count.clear();
	g_weights.clear();
	g_memo.clear();
	g_total_sum = 0;
	for (int i = 0; i < (sizeof CANDIDATES / sizeof(Candidate)); i ++) {
		Candidate candidate = CANDIDATES[i];
		int weight = candidate.weight;
		if (g_weight_to_count.find(weight) != g_weight_to_count.end()) {
			g_weight_to_count[weight] += 1;
		} else {
			g_weight_to_count[weight] = 1;
			g_weights.push_back(weight);
		}
		g_total_sum += weight;
		if (candidate.id == TARGET_ID) {
			g_target_weight = weight;
		}
		
	}
}

double recur(int times) {
	// メモ化
	// (g_weight_to_countが同じ内容ならばrecur_internalの結果も同じだから
	//  記憶しておいて使いまわす)
	if (g_memo.find(g_weight_to_count) != g_memo.end()) {
		return g_memo[g_weight_to_count];
	} else {
		return g_memo[g_weight_to_count] = recur_internal(times);
	}
}

double recur_internal(int times) {
	if (times == 0) return 0;
	double result = 0;
	for (int i = 0; i < g_weights.size(); i ++) {
		int weight = g_weights[i]; // 重みがweightの値であり未選択の候補の数
		int count = g_weight_to_count[weight];
		double pr;
		if (weight == g_target_weight) {
			pr = 1.0 / count; // TARGET_ID を選んだときの確率
			if (count > 1) {
				// TARGET_IDとweightは同じだがTARGET_IDではないものを選んだときの確率
				pr += ((double)(count - 1) / count) * choose(times, weight);
			}
		} else {
			pr = choose(times, weight);
		}
		result += pr * (count * weight) / g_total_sum;
	}
	return result;
}

// weightを選択した場合の残りtimes - 1個分での確率を計算
double choose(int times, int weight) {
	if (g_weight_to_count[weight] > 0) {
		g_weight_to_count[weight] -= 1;
		g_total_sum -= weight;
		double pr = recur(times - 1);
		g_total_sum += weight;
		g_weight_to_count[weight] += 1;
		return pr;
	} else {
		return 0;
	}
}