Glicko2

altre · Apr 12, 2020 · a2043e7 · a2043e7
1 parent 33d9504
commit a2043e7
Show file tree

Hide file tree

Showing 3 changed files with 94 additions and 7 deletions.
diff --git a/README.md b/README.md
@@ -4,3 +4,5 @@
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://altre.github.io/SkillRating.jl/dev)
 [![Build Status](https://travis-ci.com/altre/SkillRating.jl.svg?branch=master)](https://travis-ci.com/altre/SkillRating.jl)
 [![Coveralls](https://coveralls.io/repos/github/altre/SkillRating.jl/badge.svg?branch=master)](https://coveralls.io/github/altre/SkillRating.jl?branch=master)
+
+Implementation of glicko, glicko2 and elo skill rating algorithms.
diff --git a/src/SkillRating.jl b/src/SkillRating.jl
@@ -18,7 +18,10 @@ module SkillRating
     * `sjs`: Game outcomes from the perspective of the "player": 1 - won, 0.5 - draw, 0 - lost.
     
     Calculate a new glicko rating for games within a rating period or for a single game.
-    See: http://www.glicko.net/glicko/glicko.pdf
+    
+    See: http://www.glicko.net/glicko/glicko.pdf and Glickman, Mark E., "Parameter estimation 
+    in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
+    (http://www.glicko.net/research/glicko.pdf)
     """
     glicko_rating(r, RD, rjs::Vector, RDjs::Vector, sjs::Vector) = r + q/(1/RD^2 + 1/d2(r, rjs, RDjs)) * sum(g.(RDjs) .* (sjs .- Es.(r, rjs, RDjs)))
 
@@ -34,7 +37,10 @@ module SkillRating
     * `sjs`: Game outcomes from the perspective of the "player": 1 - won, 0.5 - draw, 0 - lost.
 
     Calculate a new glicko rating deviation for games within a rating period  or for a single game.
-    See: http://www.glicko.net/glicko/glicko.pdf
+    
+    See: http://www.glicko.net/glicko/glicko.pdf and 
+    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
+    (http://www.glicko.net/research/glicko.pdf)
     """
     glicko_RD(r, RD, rjs::Vector, RDjs::Vector, sjs::Vector) = sqrt(1/(1/RD^2 + 1/d2(r, rjs, RDjs)))
 
@@ -50,15 +56,21 @@ module SkillRating
     * `maxRD`: Rating deviation of a player who does not play.
 
     Calculate rating deviation after the end of a rating period.
-    See: http://www.glicko.net/glicko/glicko.pdf
+    
+    See: http://www.glicko.net/glicko/glicko.pdf and 
+    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
+    (http://www.glicko.net/research/glicko.pdf)
     """
     glicko_RD_increase(RD; c = 63.2, maxRD=350) = min(sqrt(RD^2+c^2), maxRD)
 
     """
         glicko_c(rating_perods_to_max, maxrd, medianrd)
 
     Calculate c, the increase in uncertainty between rating periods.
-    See: http://www.glicko.net/glicko/glicko.pdf
+    
+    See: http://www.glicko.net/glicko/glicko.pdf and 
+    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
+    (http://www.glicko.net/research/glicko.pdf)
     """
     glicko_c(rating_perods_to_max, maxrd, medianrd) = sqrt((maxrd^2 - medianrd^2)/rating_perods_to_max)
 
@@ -70,7 +82,10 @@ module SkillRating
     * `RB`: Rating of player B
 
     Calculate the expected outcome using glicko.
-    See: http://www.glicko.net/glicko/glicko.pdf
+    
+    See: http://www.glicko.net/glicko/glicko.pdf and 
+    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
+    (http://www.glicko.net/research/glicko.pdf)
     """
     glicko_expected(r1, RD1, r2, RD2) = 1 / (1 + 10^(-g(sqrt(RD1^2 + RD2^2)) * (r1 - r2)/400))
 
@@ -99,5 +114,74 @@ module SkillRating
     """
     elo_expected(RA, RB) = 1 / (1 + 10 ^ ((RB - RA) / 400))
 
-    export glicko_rating, glicko_RD, glicko_RD_increase, glicko_c, glicko_expected, elo_rating, elo_expected
+    Eg2(μ, μj, ϕj) = 1/(1+exp(-gg2(ϕj)*(μ - μj)))
+    gg2(ϕ) = 1/sqrt(1 + 3*(ϕ^2)/π^2)
+    compute_v(μ, μjs, ϕjs) = 1/sum(gg2.(ϕjs).^2 .* Eg2.(μ, μjs, ϕjs) .* (1 .- Eg2.(μ, μjs, ϕjs)))
+    compute_Δ(μ, μjs, ϕjs, sjs, v) = v * sum(gg2.(ϕjs) .* (sjs .- Eg2.(μ, μjs, ϕjs)))
+    ϕincrease(ϕ, σ) = sqrt(ϕ^2 + σ^2)
+
+    """
+        glicko2(μ::Real, ϕ::Real, σ::Real, μjs::Vector{Real}, ϕjs::Vector{Real}, sjs::Vector{Real}; τ=0.05, ϵ=10^-6)
+    
+    # Arguments:
+    * `μ`: Rating of player
+    * `ϕ`: Rating deviation of player
+    * `σ`: Rating volatility
+    * `μjs`: Ratings of oppenents within rating period
+    * `μjs`: Rating deviations within rating period
+    * `sjs`: Game outcomes from the perspective of the "player": 1 - won, 0.5 - draw, 0 - lost.
+    * `τ=0.5`: smaller values constrain the change in volatility over time, reasonable range between 0.3 and 1.2
+    * `ϵ=10^-6`: Convergence bound
+
+    Returns tuple `(μ, ϕ, σ)`: glicko2 rating, rating deviation and rating volatility for games in one rating period.
+
+    See: http://www.glicko.net/glicko/glicko2.pdf and 
+    Glickman, Mark E., "Dynamic paired comparison models with stochastic variances" (2001), 
+    Journal of Applied Statistics, 28, 673-689. 
+    (http://www.glicko.net/research/dpcmsv.pdf)
+    """
+    function glicko2(μ::Real, ϕ::Real, σ::Real, μjs::Vector{R}, ϕjs::Vector{R}, sjs::Vector{S}; τ=0.5, ϵ=10^-6) where {R <: Real, S <: Real}
+        v = compute_v(μ, μjs, ϕjs)
+        Δ = compute_Δ(μ, μjs, ϕjs, sjs, v) 
+        a = log(σ^2)
+        f(x) = exp(x)*(Δ^2 - ϕ^2  - v - exp(x)) / (2 * (ϕ^2 + v + exp(x))^2) - (x - a)/τ^2
+        A = a
+        if Δ^2 > ϕ^2 + v 
+            B = log(Δ^2 - ϕ^2  - v)
+        else
+            k = 1
+            while (f(a - k*τ) < 0)
+                k += 1
+            end
+            B = a - k*τ
+        end
+        fA = f(A)
+        fB = f(B)
+        while abs(B - A) > ϵ
+            C = A + (A-B) * fA/ (fB-fA)
+            fC = f(C)
+            if fC * fB < 0
+                A = B
+                fA = fB
+            else
+                fA /= 2
+            end
+            B = C
+            fB = fC
+        end
+        σnew = exp(A/2)
+        ϕincreased = sqrt(ϕ^2 + σnew^2)
+        ϕnew = 1/sqrt(1/ϕincreased^2 + 1/v)
+        μnew = μ + ϕnew^2 * sum(gg2.(ϕjs) .* (sjs .- Eg2.(μ, μjs, ϕjs)))
+        μnew, ϕnew, σnew
+    end
+
+    """
+        glicko1_to_glicko2(r::real, RD::Real)
+
+    Convert ratings `r` and rating deviations `RD` from glicko1 (Elo-like) scale to glicko 2 scale (μ and ϕ).
+    """
+    glicko1_to_glicko2(r::Real, RD::Real) = ((r-1500)/173.7178, RD/173.7178)
+
+    export glicko_rating, glicko_RD, glicko_RD_increase, glicko_c, glicko_expected, elo_rating, elo_expected, glicko2
 end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -5,7 +5,8 @@ using SkillRating
     @test round(glicko_rating(1500, 200, [1400, 1550, 1700], [30, 100, 300], [1, 0, 0])) == 1464
     @test round(glicko_RD(1500, 200, [1400, 1550, 1700], [30, 100, 300], [1, 0, 0]), digits = 1) == 151.4
     @test round(glicko_expected(1400, 80, 1500, 150), digits = 3) == 0.376
+    @test round.(glicko2(0, 1.1513, 0.06, [-0.5756, 0.2878, 1.1513], [0.1727, 0.5756, 1.7269], [1, 0, 0]; τ = 0.5), digits=4) == (-0.2069, 0.8722, 0.06)
 
     @test round(elo_rating(1613, [1609, 1477, 1388, 1586, 1720], [0,0.5,1,1,0])) == 1601
     @test round(elo_expected(1613, 1609), digits=2) == 0.51
-end
+end