2015-08-07 4 views
4

Какое наилучшее решение для нахождения вероятности прокатки суммы с n кубиками? Я решаю это, найдяВероятность получения конкретной суммы после прокатки n кубиков. Ruby

  1. средний.
  2. стандартное отклонение.
  3. z_score для чисел ниже x
  4. z_score для чисел выше x
  5. преобразования как для вероятностей
  6. вычитая одно из другого

Это то, что я сделал до сих пор.

# sides - number of sides on one die 
def get_mean(sides) 
    (1..sides).inject(:+)/sides.to_f 
end 

def get_variance(sides) 
    mean_of_squares = ((1..sides).inject {|sum, side| sum + side ** 2})/sides.to_f 
    square_mean = get_mean(sides) ** 2 

    mean_of_squares - square_mean 
end 

def get_sigma(variance) 
    variance ** 0.5 
end 

# x - the number of points in question 
def get_z_score(x, mean, sigma) 
    (x - mean)/sigma.to_f 
end 

# Converts z_score to probability 
def z_to_probability(z) 
    return 0 if z < -6.5 
    return 1 if z > 6.5 

    fact_k = 1 
    sum = 0 
    term = 1 
    k = 0 

    loop_stop = Math.exp(-23) 
    while term.abs > loop_stop do 
    term = 0.3989422804 * ((-1)**k) * (z**k)/(2*k+1)/(2**k) * (z**(k+1))/fact_k 
    sum += term 
    k += 1 
    fact_k *= k 
    end 

    sum += 0.5 
    1 - sum 
end 

# Calculate probability of getting 'х' total points by rolling 'n' dice with 'sides' number of sides. 
def probability_of_sum(x, n, sides=6) 

    mean = n * get_mean(sides) 
    variance = get_variance(sides) 
    sigma = get_sigma(n * variance) 

    # Rolling below the sum 
    z1 = get_z_score(x, mean, sigma) 
    prob_1 = z_to_probability(z1) 

    # Rolling above the sum 
    z2 = get_z_score(x+1, mean, sigma) 
    prob_2 = z_to_probability(z2) 

    prob_1 - prob_2 
end 

# Run probability for 100 dice 
puts probability_of_sum(400, 100) 

Что касается меня, когда я выбираю x = 200 вероятность равна 0. ли это правильное решение?

+0

Если я правильно понял .. нет. Если вы будете качать все 2 100, общее количество будет 200. Таким образом, есть вероятность, что это произойдет. – mtamhankar

+0

Вот что я тоже думал.Однако, согласно [68-95-99.7 (эмпирическому) правилу] (https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule), или правило 3-сигмы «Около ** 68% ** значений, полученных из нормального распределения, находятся в пределах ** 1 ** стандартное отклонение _σ_ от среднего значения, около ** 95% ** значений лежат в пределах ** 2 ** стандартных отклонений; и около ** 99,7% ** находятся в пределах ** 3 ** стандартных отклонений. " В этом случае с ** 100 ** 6-сторонними кубиками 'mean = 350' и' σ = 17'. Это означает, что ** 99,7% ** значений будет находиться в диапазоне ** 299 ** до ** 401 **. (350 +/- 17 * 3) – Pav31

+0

Также примечание: в алгоритме есть что-то принципиально неправильное, запустив что-то очень простое, например x = 12 и n = 2, которое должно быть 1/36, не работает. Я думаю, что @NeilSlater может быть связан с чем-то вроде дискретности. – mtamhankar

ответ

1

Вот моя окончательная версия.

  1. Изменены Коррекции сумм в get_z_score к x-0.5 и x+0.5 соответственно для более точного результата.
  2. Добавлен return 0 if x < n || x > n * sides для покрытия случаев, где сумма меньше количества кубиков и выше, чем количество костей, умноженное на количество сторон.
  3. тесты Benchmark Добавленные с результатом

Основные функции

# sides - number of sides on one die 
def get_mean(sides) 
    (1..sides).inject(:+)/sides.to_f 
end 

def get_variance(sides) 
    mean_of_squares = ((1..sides).inject {|sum, side| sum + side ** 2})/sides.to_f 
    square_mean = get_mean(sides) ** 2 

    mean_of_squares - square_mean 
end 

def get_sigma(variance) 
    variance ** 0.5 
end 

# x - the number of points in question 
def get_z_score(x, mean, sigma) 
    (x - mean)/sigma.to_f 
end 

# Converts z_score to probability 
def z_to_probability(z) 
    return 0 if z < -6.5 
    return 1 if z > 6.5 

    fact_k = 1 
    sum = 0 
    term = 1 
    k = 0 

    loop_stop = Math.exp(-23) 
    while term.abs > loop_stop do 
    term = 0.3989422804 * ((-1)**k) * (z**k)/(2*k+1)/(2**k) * (z**(k+1))/fact_k 
    sum += term 
    k += 1 
    fact_k *= k 
    end 

    sum += 0.5 
    1 - sum 
end 

# Calculate probability of getting 'х' total points by rolling 'n' dice with 'sides' number of sides. 
def probability_of_sum(x, n, sides=6) 
    return 0 if x < n || x > n * sides 

    mean = n * get_mean(sides) 
    variance = get_variance(sides) 
    sigma = get_sigma(n * variance) 

    # Rolling below the sum 
    z1 = get_z_score(x-0.5, mean, sigma) 
    prob_1 = z_to_probability(z1) 

    # Rolling above the sum 
    z2 = get_z_score(x+0.5, mean, sigma) 
    prob_2 = z_to_probability(z2) 

    prob_1 - prob_2 
end 

Benchmark

require 'benchmark' 

Benchmark.bm do |x| 
    x.report { @prob = probability_of_sum(350, 100).to_f } 
    puts "\tWith x = 350 and n = 100:" 
    puts "\tProbability: #{@prob}" 
end 
puts 

Benchmark.bm do |x| 
    x.report { @prob = probability_of_sum(400, 100).to_f } 
    puts "\tWith x = 400 and n = 100:" 
    puts "\tProbability: #{@prob}" 
end 
puts 

Benchmark.bm do |x| 
    x.report { @prob = probability_of_sum(1000, 300).to_f } 
    puts "\tWith x = 1000 and n = 300:" 
    puts "\tProbability: #{@prob}" 
end 

Результат

 user  system  total  real 
    0.000000 0.000000 0.000000 ( 0.000049) 
    With x = 350 and n = 100: 
    Probability: 0.023356331366255034 

     user  system  total  real 
    0.000000 0.000000 0.000000 ( 0.000049) 
    With x = 400 and n = 100: 
    Probability: 0.00032186531055478085 

     user  system  total  real 
    0.000000 0.000000 0.000000 ( 0.000032) 
    With x = 1000 and n = 300: 
    Probability: 0.003232390001131513 
+0

Привет, спасибо за решение! Не могли бы вы объяснить, что значит значение 0.3989422804 и почему loop_stop = Math.exp (-23)? – AlexMrKlim

+0

Это было взято из [здесь] (http://stackoverflow.com/questions/31875909/z-score-to-probability-and-vice-verse-in-ruby/31876457#31876457), которое было взято из [здесь] (http://stackoverflow.com/questions/16194730/seeking-a-statistical-javascript-function-to-return-p-value-from-az-score/16197404#16197404). Также вы можете попробовать найти _z-score для формулы вероятности_, ​​чтобы лучше понять. – Pav31

0

Я также решил это с помощью метода Monte Carlo, и результаты относительно близки.

# x - sum of points to find probability for 
# n - number of dice 
# trials - number of trials 
def monte_carlo(x, n, trials=10000) 
    pos = 0 

    trials.times do 
    sum = n.times.inject(0) { |sum| sum + rand(1..6) } 
    pos += 1 if sum == x 
    end 

    pos/trials.to_f 
end 

puts monte_carlo(300, 100, 30000) 
6

Добавление результат двух независимых распределений вероятности является такой же, как convolving оба распределения. Если распределения дискретны, то это дискретная свертка.

Так что, если один кубик представлен как:

probs_1d6 = Array.new(6) { Rational(1,6) } 

Затем 2d6 можно рассчитать следующим образом:

probs_2d6 = [] 
probs_1d6.each_with_index do |prob_a,i_a| 
    probs_1d6.each_with_index do |prob_b,i_b| 
    probs_2d6[i_a + i_b] = (probs_2d6[i_a + i_b] || 0) + prob_a * prob_b 
    end 
end 

probs_2d6 
# => [(1/36), (1/18), (1/12), (1/9), (5/36), (1/6), 
#  (5/36), (1/9), (1/12), (1/18), (1/36)] 

Хотя это п-квадрат для сторон кости, и вполне логично комбинация может уменьшить это, что делает это, как правило, менее гибким для более сложных настроек. Самое приятное в этом подходе - вы можете продолжать добавлять больше кубиков и делать другие более экзотические комбинации. Например, чтобы получить 4d6, вы можете сверлить два результата для 2d6. Использование рациональных решений позволяет вам решать проблемы с плавающей точкой.

Я пропустил одну деталь, вам нужно сохранить начальное смещение (+1 для нормальной шестисторонней матрицы) и добавить его вместе, чтобы узнать, с чем связаны вероятности.

Я сделал более сложную версию этой логики в gem games_dice, в плавающей точке, а не в Rational, которая может справиться с несколькими другими комбинациями костей.

Вот основной переписана вашего метода с использованием вышеуказанного подхода в наивном способе (просто комбинируя эффекты кости по одному за раз):

def probability_of_sum(x, n, sides=6) 
    return 0 if x < n 
    single_die_probs = Array.new(sides) { Rational(1,sides) } 

    combined_probs = [1] # Represents chance of getting 0 when rolling 0 dice :-) 

    # This is not the most efficient way to do this, but easier to understand 
    n.times do 
    start_probs = combined_probs 
    combined_probs = [] 
    start_probs.each_with_index do |prob_a,i_a| 
     single_die_probs.each_with_index do |prob_b,i_b| 
      combined_probs[i_a + i_b] = (combined_probs[i_a + i_b] || 0) + prob_a * prob_b 
     end 
    end 
    end 

    combined_probs[ x - n ] || 0 
end 

puts probability_of_sum(400, 100).to_f 
# => 0.0003172139126369326 

Примечания метод фактически вычисляет распределение полной вероятности 100 -to-600, поэтому вам нужно только один раз вызвать его и сохранить массив (плюс смещение +100) один раз, и вы можете делать другие полезные вещи, такие как вероятность получить больше определенного числа. Все с полной точностью благодаря использованию номеров Rational в Ruby.

Поскольку в вашей ситуации у вас есть только один тип кубиков, мы можем избежать использования Rational до конца, работая только с целыми числами (по существу, считая комбинированные значения) и делим на общее количество комбинаций (сторон мощность числа рулонов). Это намного быстрее, и возвращает значения для 100х костей в рамках второй:

def probability_of_sum(x, n, sides=6) 
    return 0 if x < n 
    combined_probs = [1] # Represents chance of getting 0 when rolling 0 dice :-) 

    n.times do 
    start_probs = combined_probs 
    combined_probs = [] 
    start_probs.each_with_index do |prob_a,i_a| 
     sides.times do |i_b| 
      combined_probs[i_a + i_b] = (combined_probs[i_a + i_b] || 0) + prob_a 
     end 
    end 
    end 

    Rational(combined_probs[ x - n ] || 0, sides ** n) 
end 
+1

Просто незначительное разъяснение. Плотность вероятности или функция массы для суммы независимых переменных не похожа на свертку функций плотности или массы, это точно свертка. –

6

Существует точное решение с контрастной суммой биномиальных коэффициентов. Я написал это в нескольких местах (по Quora и MSE), и вы можете найти его в другом месте, хотя есть некоторые ошибочные версии. Будьте осторожны, если вы это реализуете, вам может потребоваться отменить биномиальные коэффициенты, которые намного больше конечного результата, и если вы используете арифметику с плавающей запятой, вы можете потерять слишком много точности.

Предложение Neil Slater использовать динамическое программирование для вычисления свертки является хорошим. Он медленнее суммирования биномиальных коэффициентов, но достаточно надежный. Вы можете ускорить его несколькими способами, например, с помощью возведения в степень возведения в квадрат и с помощью быстрого преобразования Фурье, но многие люди найдут, что они будут излишними.

Чтобы исправить ваш метод, вы должны использовать (простую) корректуру непрерывности в нормальном приближении и ограничить контекст, в котором у вас достаточно кубиков, и вы оцениваете достаточно далеко от максимума и минимума, которые вы ожидаете от нормального приближения будет хорошим, либо в абсолютном, либо в относительном смысле. Коррекция непрерывности состоит в том, что вы заменяете счетчик n интервалом от n-1/2 до n + 1/2.

Точный подсчет количества способов прокатки 200 - это 7745954278770349845682110174816333221135826585518841002760, поэтому вероятность равна делению на 6^100, что составляет около 1,188563 x 10^-20.

Нормальное приближение с простой коррекцией непрерывности - Phi ((200,5-350)/sqrt (3500/12)) - Phi ((199,5-350)/sqrt (3500/12)) = 4,2 x 10^19. Это точно в абсолютном выражении, так как оно очень близко к 0, но оно отключено в 35 раз, поэтому оно не велико в относительном выражении. Нормальное приближение дает лучшее относительное приближение ближе к центру.