This notebook is designed to be executed in order. Some of the functions remember their previous values, so if you execute the notebook out of order it might be good to clear the functions or restart the kernel. There are optional examples included, which can be executed again or opened.

These are the supposed probabilities of winning a match of length n if the true rating difference is d points on FIBS.

fp[d_, n_] := fp[d, n] = 
        N[1/(1 + 10^(-d Sqrt[n]/2000)), 20]

Example: The estimated probability of winning a 25 point match if one is 50 points stronger.

fp[50, 25]

This gives True if a uniformly distributed on [0,1] random variable is less than p, and False otherwise. I believe that Mathematica's pseudo random number generator is good enough for these purposes, though I have not tested it in detail.

win[p_] := {Random[] < p}[[1]]

This fixes the level of precision of x to n digits after the decimal point. Otherwise, in repeated calculations Mathematica loses precision erratically far faster than it ought to, probably because it loses track of how many digits are accurate when there are a lot of zeros at the end of the window of precision.

fixprecision[x_, n_] := N[IntegerPart[x 10^n]/10^n, n]

The following take in the previous rating, the opponent's rating, the length of the match, and whether one won (True) or not (False).

fnewrat[oldrat_, opprat_, n_, wins_] := 
      fixprecision[
        oldrat + 
          4 Sqrt[n] 
            If[wins, 
              1 - fp[oldrat - opprat, n], -fp[oldrat - opprat, n]], 
        20]

The following is the force law with win probability function sp and rating snewrat of a player with strength rat1, rating rat2, and opponent strength opprat, for an n point match. I have not included other probability function and rating formulas, so sp should always be fp and snewrat should always be fnewrat.

drift[sp_,  snewrat_, rat1_, rat2_, opprat_, 
        n_] := drift[sp,  snewrat, rat1, rat2, opprat, n] = 
        sp[rat1 - opprat, n] 
            snewrat[rat2, opprat, n, True] + (1 - sp[rat1 - opprat, n]) 
            snewrat[rat2, opprat, n, False] - rat2

The following is the drift per experience point.

driftperexp[sp_,  snewrat_, rat1_, rat2_, opprat_, n_] := 
      drift[sp,  snewrat, rat1, rat2, opprat, n]/n

The following is the drift per experience point of a player with playing strength 0, rating delta, and opponent strength and rating opprat.

dperexp[sp_, snewrat_, delta_, opprat_, 
        n_] := dperexp[sp, snewrat, delta, opprat, n] = 
        driftperexp[sp,  snewrat, 0, delta, opprat, n]

The following is the variance per experience point. It should be the total stakes^2 p (1-p). The playing strength is assumed to be 0, and the oppenent has strength opprat.

varperexp[sp_, snewrat_, opprat_, 
        n_] := varperexp[sp, snewrat, opprat, 
          n] = ((snewrat[0, 0, n, True] - 
                    snewrat[0, 0, n, False])^2 sp[opprat, 
                n] (1 - sp[opprat, n]))/n

This is the k drift/variance.

ffactor[sp_, snewrat_, currat_, truerat_, opprat_, 
        n_] := ffactor[sp, snewrat, currat, truerat, opprat, n] = 
        2 dperexp[sp, snewrat, currat - truerat, opprat - truerat, n]/
            varperexp[sp, snewrat, truerat - opprat, n]

The stable distribution has probability density function pdf(x). The differential equation it satisfies is f'/f = ffactor. ctimespf is off by a constant factor, computed in normfibs.

ctimespdf[sp_, snewrat_, x_, opprat_, n_] := 
      E^NIntegrate[
          ffactor[sp, snewrat, t, 0, opprat, n], {t, 0, x}]

normfibs[opprat_, matchlength_] := normfibs[opprat, matchlength] = 
        NIntegrate[
          ctimespdf[fp, fnewrat, x, opprat, matchlength], {x, -2000, 
            2000}]

pdffibs[x_, opprat_, n_] := 
      ctimespdf[fp, fnewrat, x, opprat, n]/normfibs[opprat, n]

To find the probability that one is between rat1 and rat2, numerically integrate pdffibs[x,_,_] from x=rat1 to x=rat2.

The following gives the nth moment of the stable distribution on FIBS assuming opponents stronger by opprat and matches of length k. fibsstddev computes the standard deviation, with a miniscule error used to avoid a complaint by Mathematica when a numerical integral produces a value close to 0.

fibsnthmoment[opprat_, k_, n_] := 
      NIntegrate[x^n pdffibs[x, opprat, k], {x, -1000, 1000}]

fibsstddev[opprat_, k_] := 
      Sqrt[fibsnthmoment[opprat, k, 2] - 
          fibsnthmoment[opprat + 0.01, k, 1]^2]

Example: The standard deviation of the stable distribution when one plays 5-point matches against opponents who are correctly rated 400 points stronger.

fibsstddev[400, 5]

For the purpose of finding percentiles of the distributions, the following is Newton's Method. iternm iterates this n times on the initial guess of x.

NewtonsMethod[f_, df_, x_, targval_] := 
      x - (f[x] - targval)/df[x]

iternm[f_, df_, x_, targval_, n_] := 
NestList[
        Function[t, NewtonsMethod[f, df, t, targval]], x, n]

Opponent's rating is opprat. The match length is n. fracabove is the desired fraction of the stable distribution above the output. k is the number of iterations of Newton's Method to use. If the last two values in the output are equal (usually for k>=5) then Newton's Method has converged and you can trust the last value given.

fibspercentile[opprat_, n_, fracabove_, k_] := 
      iternm[Function[t, NIntegrate[pdffibs[x, opprat, n], {x, t, 1000}]], 
        Function[t, -pdffibs[t, opprat, n]], 0, fracabove, k]

The following computes the 90th percentile of the stable distribution. 10% of the time, one should be overrated by at least 53.4739 points if one only plays 5-point matches against people who are correctly rated at one's true rating.

fibspercentile[0, 5, 1/10, 5]

Another example. This computes the 30th percentile when one plays 1 point matches against players who are correctly rated at 100 points stronger.

fibspercentile[100, 1, 7/10, 5]

When one plays people who are stronger or weaker, the mode of the stable distribution is one's correct rating, but not the median.

fibspercentile[100, 25, 1/2, 5]