A Method for Estimating the Upper of Russian Roulette

Analyze Problem

  • In essence, the RR(Russian Roulette) is a special situation of the Bernoulli experiment, namely the Geometric distribution (repeat X times experiments with probability P). Therefore, the expectation of Geometric distribution is E = 1/p and we can use the expectation as the estimator.
  • However, the expectation is just an average number of a function, whereas we need an estimator to estimate the upper.
  • Therefore, I decide to use a Rejection-Region-based method to estimate the maximum recursion depth of RR. Also, I implemented and visualized this method with Python and the code will append in the end.

Build Mathematical Model

  • In statistics, an event that occurs with a very small probability α to happen could be seen as an impossible event. More specifically, such a very small probability α is named Significance Level in Unilateral Hypothesis testing.
  • Given that x is the current recursion depth and p is the success probability, we have:

P(continue_RR) = pow(p, x).

  • If we set the P(continue_RR) to a very small probability α (Significance Level). Then, we gain the equation:

x = log(p, α), p is the base

  • So, if given a p, we can estimate x now. For example: (p = 0.9, α = 0.001) >> x ≈ 66 (ceil)

Click here to view model code

Python
import numpy as np
from matplotlib import pyplot as plt
import math
from sys import float_info
import os

def decimalPlacesNum(num)->int:
    temp = int(num)
    dpn = 0
    while(temp != num):
        num *= 10
        temp = int(num)
        dpn += 1
    return dpn

# Parameters
SHOW_MAX_ITER_TIMES = 100
SIGNIFICANCE_LEVEL = 0.001
PLOT_TITLE = "Russian Roulette"
SAVE_PLOT = True
prob_RR = 0.9

# Program Begins
if __name__ == "__main__":
    plt.ylim(0,SHOW_MAX_ITER_TIMES)
    plt.xlim(0,1.0)
    
    stride = min(0.01, 1.0 / (10 ** decimalPlacesNum(prob_RR)))
    x = np.arange(float_info.epsilon,1, stride)
    y = list( map(lambda e : math.ceil(math.log(SIGNIFICANCE_LEVEL, e)), x) )
    
    idx_nearest_prob = min(np.searchsorted(x, prob_RR, side = "right"), len(x)-1)
    nearest_prob = x[idx_nearest_prob]
    expected_iter_times = y[idx_nearest_prob]
    
    # Clamp Point Position Range & Draw
    point_position = (nearest_prob, min(expected_iter_times, SHOW_MAX_ITER_TIMES))
    plt.scatter([point_position[0]],[point_position[1]],color="red")
    # Auxiliary lines & Annotation
    plt.plot([point_position[0], point_position[0]],[0,point_position[1]],linestyle="--")
    plt.plot([0, point_position[0]],[point_position[1],point_position[1]],linestyle="--")
    
    if (prob_RR >= 1.0): annotation = "(1,∞)"
    else: annotation = "({},{})".format(nearest_prob - float_info.epsilon, expected_iter_times)
    
    plt.annotate(annotation,
                 xy = (point_position[0], point_position[1]), xycoords = "data",
                 xytext = (-len(annotation)*(len(annotation)//3), 10), textcoords = "offset points",
                 fontsize = 12)
    
    plt.title(PLOT_TITLE)
    plt.xlabel("Probobility")
    plt.ylabel("Theoretical Upper Iteration Times")
    plt.plot(x,y)
    if SAVE_PLOT: 
        print("[TIPS]: Plot has been saved in {}".format(os.path.join(os.path.abspath(os.path.curdir), PLOT_TITLE)))
        plt.savefig(PLOT_TITLE, dpi=300)
    plt.show()

Leave a Reply

Your email address will not be published. Required fields are marked *