Calculating the constant C in Dota 2 pseudo-random distribution

http://www.playdota.com/forums/showthread.php?t=7993

Maybe this thread will help you. A "numerical" approach is mentioned but unfortunately I cannot find the source of it. If you want to get the actual values of C you have to test them for various %-based skills. Maybe all the percentages are not represented so you might have to infer the other values for C, but it's actually pointless if there are no skills/items corresponding.

For example, if you want to have the actual C for P(E)=15%, you take Coup de Grâce, PA's ultimate and you measure the number of time you crit right after a crit divided by the number of attacks you made after a crit. You should find a value close to 3.22% like stated in the wiki. In the playdota thread it is mentioned they have done a million attacks.

I'm not sure whether the "theoretical C" is relevant here if you are able to do some testing you should find the actual C and that is the only one that matters.

If you know about Markov chains you can basically solve the P(E) for any C value. You need to think about a bunch of states {1,2,3...}, which correspond to how many "stacks" of C you've accumulated.

So you have a matrix P where the nth row and kth column holds the probability that you go from n stacks to k stacks. However, the only valid transitions are:

• From n stacks to 1 stack when you proc and it resets.
• From n stacks to n+1 stacks when you don't proc.

The matrix is of finite size because once you reach n = ⌈1/C⌉ (rounded up) stacks, the probability of going back to 1 stack is 100%.

Now with your A matrix constructed, you can compute what is known as a stationary distribution, where aX = a, where a is a vector and X is a matrix. If you have Matlab you can solve for a by taking the eigenvector of the eigenvalue 1. This eigenvector represents the proportion of time spent in each state, so by dividing it by the sum of all the values, you obtain a long term proportion of time spent in each state.

Once you know the probability of being in each state you can multiply it by the probability of a proc in that state and sum it across states to get the average proc rate. Then you just try different C values until you achieve desired P(E).

Below is some Matlab code I put together, the values are not exactly the same as the wiki page by they are sufficiently similar that I have not put more thought into it.

function [C, prob] = pseudorand(target) % Takes in a target p and outputs the C value
    C = 0.30210;                    % Initial points for iterating
    prob = 0.5;
while abs(target-prob)>10^-8        % Iterate and adjust C until desire p is 
                                      achieved.
    if C<0.1
        C = C + (target-prob)/8;    % Be more careful around small C values
    else
        C = C + (target-prob)/2;
    end


    P = zeros(ceil(1/C));           % Determine the size of the matrix.

    for n = 1:size(P,1)             % Set values of first column to n*C, rounding
        if n*C < 1                    down to 1 on the final row.
            P(n,1) = n*C;
        else
            P(n,1) = 1;
        end
    end                            

    for n = 1:size(P,1)-1
        P(n,n+1) = 1-P(n,1);        % Populating entries of nth row and (n+1)th 
    end                               column with value 1-n*C (once again with 
                                      rounding on the last row).

    [v, d] = eig(P');               % Obtain eigenvectors and eigenvalues.
    stationary = v(:,1)/sum(v(:,1));% Take first eigenvector and normalise.

    prob = stationary'*P(:,1);      % Find P(E) for corresponding C value.
end

EDIT (response to comment below): I'm sorry but I don't quite understand what you mean by "testing values". This process rests on the assumption that the pseudo-random probability satisfies p = n × C. The assumption is that the Warcraft 3 developers wrongly calculated the C values which lead to a inconsistency between expected (in the tooltip) and observed. So the true values can either be found by taking the observed (through parsing) probability to calculate the corresponding C value using my algorithm or you could somehow acquire the table used for C values in the source code, which is from where I assume the wiki got them from.

The other solution would to be parse in a way that discriminates between the number of misses experienced, for example you can calculate the frequency of procs after 2,3,4 non-procs then fit a straight line across the points to determine the linear growth of the C value. This method is probably the most cumbersome but provides some good granularity and also gives some indictation as to whether or not the probability is indeed p = n × C (ignoring the point where n × C goes above 1, the values should pretty much sit in a straight line).

user72955's iterative approach is fundamentally along the right lines, but you could run into problems if you ported that Matlab code to other systems with different floating point precision, or if you wanted the numbers calculated to the limits of precision of your programming environment.

In general, when you are testing for convergence on a system where floating point numbers are inexactly represent by binary formats, you will not want to make the comparison between expected and actual with an arbitrary error value as in that code where it tests abs(target-prob)>10^-8. Testing this way is prone to getting caught in an infinite loop if the iteration converges to a value that is not as accurate as your are guessing it will be. In the code I give below, I test the outputs of consecutive iterations for convergence to 0 as my terminating condition, rather than latest output vs expect value, to avoid this issue. Because the calculation of P from C is deterministic, and precision the number format is finite, it will always converge to a point where consecutive outputs are identical and as close as possible to the expected value, without the possibility of getting caught in a cycle or failing to converge due to binary representation inaccuracies.

The second issue I have with that solution is the incrementing of C by adding an arbitrary percentage of the output error. This is not a linear gain system where additive negative feedback will quickly and accurately converge, as the poster realized when they had to put in the divide-by-8 fudge factor for smaller increments for low values of C. If more precision were asked for, more fudge factors would be needed to keep the system convergent. It's better not to make assumptions about the relation of the output to the input, and adjust your input by an interval halving approach.

Here is a solution in C# code that will converge to the limit of the precision of the decimal format type. The 'm' after numbers just tells the compiler they are decimal format literals, as opposed to floats or doubles.

    public decimal CfromP( decimal p )
    {
        decimal Cupper = p;
        decimal Clower = 0m;
        decimal Cmid;
        decimal p1;
        decimal p2 = 1m;
        while(true)
        {
            Cmid = ( Cupper + Clower ) / 2m;
            p1 = PfromC( Cmid );
            if ( Math.Abs( p1 - p2 ) <= 0m ) break;

            if ( p1 > p )
            {
                Cupper = Cmid;
            }
            else
            {
                Clower = Cmid;
            }

            p2 = p1;
        }

        return Cmid;
    }

    private decimal PfromC( decimal C )
    {
        decimal pProcOnN = 0m;
        decimal pProcByN = 0m;
        decimal sumNpProcOnN = 0m;

        int maxFails = (int)Math.Ceiling( 1m / C );
        for (int N = 1; N <= maxFails; ++N)
        {
            pProcOnN = Math.Min( 1m, N * C ) * (1m - pProcByN);
            pProcByN += pProcOnN;
            sumNpProcOnN += N * pProcOnN;
        }

        return ( 1m / sumNpProcOnN );
    }

CfromP() is the main function you call, and PfromC is a helper function that forward computes P (the easy direction to go).

It outputs this:

C(0.05) = 0.003801658303553139101756466
C(0.10) = 0.014745844781072675877050816
C(0.15) = 0.032220914373087674975117359
C(0.20) = 0.055704042949781851858398652
C(0.25) = 0.084744091852316990275274806
C(0.30) = 0.118949192725403987583755553
C(0.35) = 0.157983098125747077557540462
C(0.40) = 0.201547413607754017070679639
C(0.45) = 0.249306998440163189714677100
C(0.50) = 0.302103025348741965169160432
C(0.55) = 0.360397850933168697104686803
C(0.60) = 0.422649730810374235490851220
C(0.65) = 0.481125478337229174401911323
C(0.70) = 0.571428571428571428571428572
C(0.75) = 0.666666666666666666666666667
C(0.80) = 0.750000000000000000000000000
C(0.85) = 0.823529411764705882352941177
C(0.90) = 0.888888888888888888888888889
C(0.95) = 0.947368421052631578947368421

This isn't adding any new thought to this discussion, but here's the above C# code translated to Python. I also added a function to test a value of C against many iterations with random rolls, and the variance between runs easily supports the theoretical value for C.

#!/usr/bin/env python
import math
import random

def PfromC(C):
    if not isinstance(C, float):
        C = float(C)

    pProcOnN = 0.0
    pProcByN = 0.0
    sumNpProcOnN = 0.0

    maxFails = int(math.ceil(1.0 / C))

    for N in range(1, maxFails + 1):
        pProcOnN = min(1.0, N * C) * (1.0 - pProcByN)
        pProcByN += pProcOnN
        sumNpProcOnN += N * pProcOnN

    return (1.0 / sumNpProcOnN)

def CfromP(p):
    if not isinstance(p, float):
        p = float(p) # double precision

    Cupper = p
    Clower = 0.0
    Cmid = 0.0
    p2 = 1.0

    while(True):
        Cmid = (Cupper + Clower) / 2.0
        p1 = PfromC(Cmid)

        if math.fabs(p1 - p2) <= 0:
            break

        if p1 > p:
            Cupper = Cmid
        else:
            Clower = Cmid

        p2 = p1

    return Cmid

def testPfromC(C, iterations):
    if not isinstance(C, float):
        C = float(C)

    nsum = 0

    for i in range(iterations):
        n = 1

        while (C * n) < 1.0 and random.random() > (C * n):
            n += 1

        nsum += n

    av_iters = float(nsum) / iterations
    return 100.0 / av_iters

if __name__ == "__main__":
    iterations = 100000

    for i in range(1, 100):
        p = i / 100.0
        C = CfromP(p)
        print "%.2ft%ft%f" % (p, C, testPfromC(C, iterations))

The above code gives the following values for p, theoretical C, and experimental p:

0.01    0.000156    1.002238
0.02    0.000620    2.002579
0.03    0.001386    3.001511
0.04    0.002449    4.002288
0.05    0.003802    4.994758
0.06    0.005440    6.005638
0.07    0.007359    7.003934
0.08    0.009552    8.010426
0.09    0.012016    8.987512
0.10    0.014746    10.014611
0.11    0.017736    11.004199
0.12    0.020983    12.005177
0.13    0.024482    13.024971
0.14    0.028230    13.975791
0.15    0.032221    14.959288
0.16    0.036452    15.999053
0.17    0.040920    17.009059
0.18    0.045620    17.972132
0.19    0.050549    18.976664
0.20    0.055704    19.954822
0.21    0.061081    20.943636
0.22    0.066676    22.012710
0.23    0.072488    22.979420
0.24    0.078511    24.014735
0.25    0.084744    25.031101
0.26    0.091183    25.959664
0.27    0.097826    27.102915
0.28    0.104670    28.023057
0.29    0.111712    29.000554
0.30    0.118949    29.959944
0.31    0.126379    31.002753
0.32    0.134001    31.939392
0.33    0.141805    32.999597
0.34    0.149810    34.025757
0.35    0.157983    34.908278
0.36    0.166329    35.980153
0.37    0.174909    37.053368
0.38    0.183625    38.058557
0.39    0.192486    38.999431
0.40    0.201547    40.008482
0.41    0.210920    41.034227
0.42    0.220365    42.099236
0.43    0.229899    42.969354
0.44    0.239540    44.056551
0.45    0.249307    45.176913
0.46    0.259872    45.906516
0.47    0.270453    46.863211
0.48    0.281008    48.012061
0.49    0.291552    48.983830
0.50    0.302103    49.938326
0.51    0.312677    51.091310
0.52    0.323291    51.903294
0.53    0.334120    53.118028
0.54    0.347370    54.012304
0.55    0.360398    54.951094
0.56    0.373217    56.077072
0.57    0.385840    56.917788
0.58    0.398278    57.943470
0.59    0.410545    58.905304
0.60    0.422650    59.988002
0.61    0.434604    61.122079
0.62    0.446419    62.023197
0.63    0.458104    63.044547
0.64    0.469670    63.999181
0.65    0.481125    65.143609
0.66    0.492481    66.020110
0.67    0.507463    66.997186
0.68    0.529412    68.085570
0.69    0.550725    68.980741
0.70    0.571429    70.037330
0.71    0.591549    71.011632
0.72    0.611111    72.145387
0.73    0.630137    73.091401
0.74    0.648649    73.969421
0.75    0.666667    74.791519
0.76    0.684211    75.818460
0.77    0.701299    77.138471
0.78    0.717949    77.965415
0.79    0.734177    79.041386
0.80    0.750000    80.042903
0.81    0.765432    80.993302
0.82    0.780488    81.967885
0.83    0.795181    82.964831
0.84    0.809524    84.091559
0.85    0.823529    85.107832
0.86    0.837209    86.080003
0.87    0.850575    87.164960
0.88    0.863636    88.069258
0.89    0.876404    88.907065
0.90    0.888889    90.034934
0.91    0.901099    91.098742
0.92    0.913043    91.925283
0.93    0.924731    93.009413
0.94    0.936170    94.065413
0.95    0.947368    95.073302
0.96    0.958333    95.967448
0.97    0.969072    96.931159
0.98    0.979592    97.973900
0.99    0.989899    98.999119

The event of interest is an attack proc (a critical hit, a bash, etc) that has a certain probability of occurring for every attempt. The pseudo-random probability c is a downward adjusted probability to the fixed probability p, i.e. c < p.

I have a probability of c to proc on my first attempt. If I do not proc, and this happens with probability 1-c, then the probability of me proc-ing on my second attempt will increase by c to 2c. Eventually, nc >= 1, and the integer n is the maximum number of attempts made for a proc to register. I will be guaranteed to register a proc after n attempts, i.e. nc=1.

Let X be the random variable that represents the number of attempts before a proc is registered for a pseudo-random probability, c. The probability mass function of X is given as follows:

P(1) = c

P(2) = (1-c) x 2c

⋮

P(n-1) = (1-c) x (1-2c) x … x (1-(n-2)c) x (n-1)c

P(n) = (1-c) x (1-2c) x … x (1-(n-1)c)

where the maximum number of attempts n = ⌈1/C⌉.

The expected value E(X) gives the average number of attempts made before proc-ing.

E(X) = c + 2^2(1-c)c + … + (n-1)^2(1-c)…(1-(n-2)c)(n-1)c + n(1-c)…(1-(n-1)c)

Equate E(X) with 1/p, to get the equality

c + 2^2(1-c)c + … + (n-1)^2(1-c)…(1-(n-2)c)(n-1)c + n(1-c)…(1-(n-1)c) = 1/p.

The LHS of the equality is a polynomial of degree n-1.

The intuition behind the RHS 1/p is as follows: suppose p = 1/4, on average I need to make 4 (1/p = 1/(1/4)) attempts before the event of interest occurs.

Given a pseudo-random probability c, we can find its counter part p by taking the reciprocal of the polynomial evaluated at c. Unfortunately, I do not know the closed-form solution of p in terms of c. If p is given, then I suggest using a numerical solver to obtain c.

Implementation code in Python

import numpy as np
import pandas as pd

def c_to_p(c):
    
    ev = c
    
    prob = c
    
    n = int(np.ceil(1/c))
    
    cum_prod = 1
    
    for x in range(2, n):

        cum_prod *= 1-(x-1)*c
        prob_x = cum_prod*(x*c)
        prob += prob_x
        ev += x*prob_x
        
             
    prob_x = 1-prob
    ev += n*prob_x
        
        
    return 1/ev

# read the tabulated figure from the source and compare it with p computed by the formula

df = pd.read_html('https://dota2.gamepedia.com/Random_distribution')[0]
df['Self computed p'] = df['C'].apply(c_to_p)

Output (from Jupyter)