 > #start of original code
 > norm = require("gaussian")
 > norm = norm(0,1)
-> pdf = norm.pdf
-> cdf = norm.cdf
-> icdf = norm.ppf    # inverse CDF
-> # Update rules for approximate marginals for the win and draw cases,
-> # respectively.
-> Vwin=(t, e)->
->   return pdf(t-e) / cdf(t-e)
-> Wwin=(t, e)->
->   return Vwin(t, e) * (Vwin(t, e) + t - e)
-> Vdraw=(t, e)->
->   return (pdf(-e-t) - pdf(e-t)) / (cdf(e-t) - cdf(-e-t))
-> Wdraw=(t, e)->
->   return pow(Vdraw(t, e) , 2) + ((e-t) * pdf(e-t) + (e+t) * pdf(e+t)) / (cdf(e-t) - cdf(-e-t))
-> class Gaussian
->   ###
->   Object representing a gaussian distribution.  Create as:
->     new Gaussian({mu=..., sigma=...})
->       or
->     new Gaussian({pi=..., tau=...})
->       or
->     new Gaussian()    # gives 0 mean, infinite sigma
->   ###
->   constructor:(parms={})->
+51c54
+<     if "pi" in parms
+---
 >     if parms.pi != undefined
->       @pi = parms["pi"]
->       @tau = parms["tau"]
+54c57
+<     else if "mu" in parms
+---
 >     else if parms.mu != undefined
->       @pi = pow(parms["sigma"] , -2)
->       @tau = @pi * parms["mu"]
->     else
->       @pi = 0
->       @tau = 0
+60,68c63
+<   repr:()->
+<     return "N(pi={0.pi},tau={0.tau})".format(this)
+<   toString:()->
+<     if @pi == 0.0
+<       return "N(mu=0,sigma=inf)"
+<     else
+<       sigma = sqrt(1/@pi)
+<       mu = @tau / @pi
+<       return "N(mu={0:.3f},sigma={1:.3f})".format(mu, sigma)
+---
 >     if isNaN(@pi) || isNaN(@tau) then throw new Error "Gaussian parms can not be NaN"
->   MuSigma:()->
->     ### Return the value of this object as a (mu, sigma) tuple. ###
->     if @pi == 0.0
+72c67
+<       return 0, float("inf")
+---
 >       return [0, Infinity]
->     else
+74c69
+<       return @tau / @pi, sqrt(1/@pi)
+---
 >       return [@tau / @pi, sqrt(1/@pi)]
->   mul:(other)->
->     return new Gaussian({"pi":@pi+other.pi, "tau":@tau+other.tau})
->   div:(other)->
->     return new Gaussian({"pi":@pi-other.pi, "tau":@tau-other.tau})
-> class Variable
->   ### A variable node in the factor graph. ###
->   constructor:()->
->     @value = new Gaussian()
->     @factors = {}
->   AttachFactor:(factor)->
->     @factors[factor] = new Gaussian()
->   UpdateMessage:(factor, message)->
->     old_message = @factors[factor]
->     @value = @value.div(old_message).mul(message)
->     @factors[factor] = message
->   UpdateValue:(factor, value)->
->     old_message = @factors[factor]
->     @factors[factor] = value.mul(old_message).div(@value)
->     @value = value
->   GetMessage:(factor)->
->     return @factors[factor]
-> class Factor
->   ### Base class for a factor node in the factor graph. ###
->   constructor:(variables)->
+98a94
 >     @id = genId()
->     @variables = variables
->     for v in variables
->       v.AttachFactor(this)
+101a98,99
 >   toString:()->"Factor_#{@id}"
 > 
-> # The following Factor classes implement the five update equations
-> # from Table 1 of the Herbrich et al. paper.
-> class PriorFactor extends Factor
->   ### Connects to a single variable, pushing a fixed (Gaussian) value
->   to that variable. ###
->   constructor:(variable, param)->
->     super([variable])
->     @param = param
->   Start:()->
->     @variables[0].UpdateValue(this, @param)
-> class LikelihoodFactor extends Factor
->   ### Connects two variables, the value of one being the mean of the
->   message sent to the other. ###
->   constructor:(mean_variable, value_variable, variance)->
->     super([mean_variable, value_variable])
->     @mean = mean_variable
->     @value = value_variable
->     @variance = variance
->   UpdateValue:()->
->     ### Update the value after a change in the mean (going "down" in
->     the TrueSkill factor graph. ###
->     y = @mean.value
->     fy = @mean.GetMessage(this)
->     a = 1.0 / (1.0 + @variance * (y.pi - fy.pi))
->     @value.UpdateMessage(this, new Gaussian({"pi":a*(y.pi - fy.pi), "tau":a*(y.tau - fy.tau)}))
->   UpdateMean:()->
->     ### Update the mean after a change in the value (going "up" in
->     the TrueSkill factor graph. ###
->     # Note this is the same as UpdateValue, with @mean and
->     # @value interchanged.
->     x = @value.value
->     fx = @value.GetMessage(this)
->     a = 1.0 / (1.0 + @variance * (x.pi - fx.pi))
->     @mean.UpdateMessage(this, new Gaussian({"pi":a*(x.pi - fx.pi),"tau":a*(x.tau - fx.tau)}))
-> class SumFactor extends Factor
->   ### A factor that connects a sum variable with 1 or more terms,
->   which are summed after being multiplied by fixed (real)
->   coefficients. ###
->   constructor:(sum_variable, terms_variables, coeffs)->
+141c139
+<     assert len(terms_variables) == len(coeffs)
+---
 >     if len(terms_variables) != len(coeffs) then throw new Error "assert error"
->     @sum = sum_variable
->     @terms = terms_variables
->     @coeffs = coeffs
+145c143
+<     super([sum_variable] + terms_variables)
+---
 >     super([sum_variable].concat(terms_variables))
->   _InternalUpdate:(variable, y, fy, a)->
+147c145,146
+<     new_pi = 1.0 / (sum(pow(a[j],2) / (y[j].pi - fy[j].pi) for j in range(len(a))))
+---
 >     new_pi = map(range(len(a)), (j)->pow(a[j],2) / (y[j].pi - fy[j].pi))
 >     new_pi = 1.0 / sum(new_pi)
->     new_tau = new_pi * sum(a[j] * (y[j].tau - fy[j].tau) / (y[j].pi - fy[j].pi) for j in range(len(a)))
->     variable.UpdateMessage(this, new Gaussian({"pi":new_pi, "tau":new_tau}))
+149a149
 >     return
->   UpdateSum:()->
->     ### Update the sum value ("down" in the factor graph). ###
+152,153c152,153
+<     y = [t.value for t in @terms]
+<     fy = [t.GetMessage(this) for t in @terms]
+---
 >     y = (t.value for t in @terms)
 >     fy = (t.GetMessage(this) for t in @terms)
->     a = @coeffs
->     @_InternalUpdate(@sum, y, fy, a)
+155a156
 >     return
->   UpdateTerm:(index)->
->     ### Update one of the term values ("up" in the factor graph). ###
->     # Swap the coefficients around to make the term we want to update
->     # be the 'sum' of the other terms and the factor's sum, eg.,
->     # change:
->     #
->     #    x = y_1 + y_2 + y_3
->     #
->     # to
->     #
->     #    y_2 = x - y_1 - y_3
->     #
->     # then use the same update equation as for UpdateSum.
->     b = @coeffs
+170,171c171,172
+<     a = [-b[i] / b[index] for i in range(len(b)) if i != index]
+<     a.insert(index, 1.0 / b[index])
+---
 >     a = ((-b[i] / b[index]) for i in range(len(b)) when i != index)
 >     a.splice(index, 0, 1.0 / b[index])
->     v = @terms.slice(0)
->     v[index] = @sum
+174,175c175,176
+<     y = [i.value for i in v]
+<     fy = [i.GetMessage(this) for i in v]
+---
 >     y = (i.value for i in v)
 >     fy = (i.GetMessage(this) for i in v)
->     @_InternalUpdate(@terms[index], y, fy, a)
-> class TruncateFactor extends Factor
->   ### A factor for (approximately) truncating the team difference
->   distribution based on a win or a draw (the choice of which is
->   determined by the functions you pass as V and W). ###
->   constructor:(variable, V, W, epsilon)->
->     super([variable])
->     @var = variable
->     @V = V
->     @W = W
->     @epsilon = epsilon
->   Update:()->
->     x = @var.value
->     fx = @var.GetMessage(this)
->     c = x.pi - fx.pi
->     d = x.tau - fx.tau
->     sqrt_c = sqrt(c)
+193,195c194,196
+<     args = (d / sqrt_c, @epsilon * sqrt_c)
+<     V = @V(*args)
+<     W = @W(*args)
+---
 >     args = [d / sqrt_c, @epsilon * sqrt_c]
 >     V = @V.apply(this,Array.prototype.slice.call(args, 0))
 >     W = @W.apply(this,Array.prototype.slice.call(args, 0))
->     new_val = new Gaussian({"pi":c / (1.0 - W), "tau":(d + sqrt_c * V) / (1.0 - W)})
->     @var.UpdateValue(this, new_val)
-> DrawProbability=(epsilon, beta, total_players=2)->
->   ### Compute the draw probability given the draw margin (epsilon). ###
->   return 2 * cdf(epsilon / (sqrt(total_players) * beta)) - 1
-> DrawMargin=(p, beta, total_players=2)->
->   ### Compute the draw margin (epsilon) given the draw probability. ###
->   return icdf((p+1.0)/2) * sqrt(total_players) * beta
-> INITIAL_MU = 25.0
-> INITIAL_SIGMA = INITIAL_MU / 3.0
+205a207,209
 > BETA = null
 > EPSILON = null
 > GAMMA = null
-> SetParameters=(beta=null, epsilon=null, draw_probability=null, gamma=null)->
->   ###
->   Sets three global parameters used in the TrueSkill algorithm.
->   beta is a measure of how random the game is.  You can think of it as
->   the difference in skill (mean) needed for the better player to have
->   an ~80% chance of winning.  A high value means the game is more
->   random (I need to be *much* better than you to consistently overcome
->   the randomness of the game and beat you 80% of the time); a low
->   value is less random (a slight edge in skill is enough to win
->   consistently).  The default value of beta is half of INITIAL_SIGMA
->   (the value suggested by the Herbrich et al. paper).
->   epsilon is a measure of how common draws are.  Instead of specifying
->   epsilon directly you can pass draw_probability instead (a number
->   from 0 to 1, saying what fraction of games end in draws), and
->   epsilon will be determined from that.  The default epsilon
->   corresponds to a draw probability of 0.1 (10%).  (You should pass a
->   value for either epsilon or draw_probability, not both.)
->   gamma is a small amount by which a player's uncertainty (sigma) is
->   increased prior to the start of each game.  This allows us to
->   account for skills that vary over time; the effect of old games
->   on the estimate will slowly disappear unless reinforced by evidence
->   from new games.
->   ###
->   if beta is null
->     BETA = INITIAL_SIGMA / 2.0
->   else
->     BETA = beta
->   if epsilon is null
->     if draw_probability is null
->       draw_probability = 0.10
->     EPSILON = DrawMargin(draw_probability, BETA)
->   else
->     EPSILON = epsilon
->   if gamma is null
->     GAMMA = INITIAL_SIGMA / 100.0
->   else
->     GAMMA = gamma
-> SetParameters()
-> AdjustPlayers=(players)->
->   ###
->   Adjust the skills of a list of players.
->   'players' is a list of player objects, for all the players who
->   participated in a single game.  A 'player object' is any object with
->   a "skill" attribute (a (mu, sigma) tuple) and a "rank" attribute.
->   Lower ranks are better; the lowest rank is the overall winner of the
->   game.  Equal ranks mean that the two players drew.
->   This function updates all the "skill" attributes of the player
->   objects to reflect the outcome of the game.  The input list is not
->   altered.
->   ###
->   players = players.slice(0)
->   # Sort players by rank, the factor graph will connect adjacent team
->   # performance variables.
+229d232
+<   global BETA, EPSILON, GAMMA
+260c263
+<   players.sort(key=lambda p: p.rank)
+---
 >   players.sort((a,b)->a.rank-b.rank)
->   # Create all the variable nodes in the graph.  "Teams" are each a
->   # single player; there's a one-to-one correspondence between players
->   # and teams.  (It would be straightforward to make multiplayer
->   # teams, but it's not needed for my current purposes.)
+265,268c268,271
+<   ss = [new Variable() for p in players]
+<   ps = [new Variable() for p in players]
+<   ts = [new Variable() for p in players]
+<   ds = [new Variable() for p in players.slice(0,-1)]
+---
 >   ss = (new Variable() for p in players)
 >   ps = (new Variable() for p in players)
 >   ts = (new Variable() for p in players)
 >   ds = (new Variable() for p in players.slice(0,-1))
->   # Create each layer of factor nodes.  At the top we have priors
->   # initialized to the player's current skill estimate.
+271,274c274,277
+<   skill = [new PriorFactor(s, new Gaussian({"mu":pl.skill[0], "sigma":pl.skill[1] + GAMMA})) for (s, pl) in zip(ss, players)]
+<   skill_to_perf = [new LikelihoodFactor(s, p, pow(BETA,2)) for (s, p) in zip(ss, ps)]
+<   perf_to_team = [new SumFactor(t, [p], [1]) for (p, t) in zip(ps, ts)]
+<   team_diff = [new SumFactor(d, [t1, t2], [+1, -1]) for (d, t1, t2) in zip(ds, ts.slice(0,-1), ts.slice(1))]
+---
 >   skill = (new PriorFactor(i.s, new Gaussian({"mu":i.pl.skill[0], "sigma":i.pl.skill[1] + GAMMA})) for i in zip({s:ss, pl:players}))
 >   skill_to_perf = (new LikelihoodFactor(i.s, i.p, pow(BETA,2)) for i in zip({s:ss, p:ps}))
 >   perf_to_team = (new SumFactor(i.t, [i.p], [1]) for i in zip({p:ps, t:ts}))
 >   team_diff = (new SumFactor(i.d, [i.t1, i.t2], [+1, -1]) for i in zip({d:ds, t1:ts.slice(0,-1), t2:ts.slice(1)}))
->   # At the bottom we connect adjacent teams with a 'win' or 'draw'
->   # factor, as determined by the rank values.
+277c280
+<   trunc = [new TruncateFactor(d, Vdraw if pl1.rank == pl2.rank else Vwin, Wdraw if pl1.rank == pl2.rank else Wwin, EPSILON) for (d, pl1, pl2) in zip(ds, players.slice(0,-1), players.slice(1))]
+---
 >   trunc = (new TruncateFactor(i.d, (if i.pl1.rank == i.pl2.rank then Vdraw else Vwin), (if i.pl1.rank == i.pl2.rank then Wdraw else Wwin), EPSILON) for i in zip({d:ds, pl1:players.slice(0,-1), pl2:players.slice(1)}))
->   # Start evaluating the graph by pushing messages 'down' from the
->   # priors.
->   for f in skill
->     f.Start()
->   for f in skill_to_perf
->     f.UpdateValue()
->   for f in perf_to_team
->     f.UpdateSum()
->   # Because the truncation factors are approximate, we iterate,
->   # adjusting the team performance (t) and team difference (d)
->   # variables until they converge.  In practice this seems to happen
->   # very quickly, so I just do a fixed number of iterations.
->   #
->   # This order of evaluation is given by the numbered arrows in Figure
->   # 1 of the Herbrich paper.
->   for i in range(5)
->     for f in team_diff
->       f.UpdateSum()             # arrows (1) and (4)
->     for f in trunc
->       f.Update()                # arrows (2) and (5)
->     for f in team_diff
->       f.UpdateTerm(0)           # arrows (3) and (6)
->       f.UpdateTerm(1)
->   # Now we push messages back up the graph, from the teams back to the
->   # player skills.
->   for f in perf_to_team
->     f.UpdateTerm(0)
->   for f in skill_to_perf
->     f.UpdateMean()
->   # Finally, the players' new skills are the new values of the s
->   # variables.
+309,311c312,321
+<   for s, pl in zip(ss, players)
+<     pl.skill = s.value.MuSigma()
+< __all__ = ["AdjustPlayers", "SetParameters", "INITIAL_MU", "INITIAL_SIGMA"]
+\ No newline at end of file
+---
 >   #console.log require('util').inspect(ss, true, 10, true)
 >   for i in zip({s:ss, pl:players})
 >     i.pl.skill = i.s.value.MuSigma()

< # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
< # See the License for the specific language governing permissions and
< # limitations under the License.
< ###
< Implements the player skill estimation algorithm from Herbrich et al.,
< "TrueSkill(TM): A Bayesian Skill Rating System".
< ###
< from __future__ import print_function
< __author__ = "Doug Zongker <dougz@isotropic.org>"
< import sys

1,14d0
< #!/usr/bin/python
< # Copyright 2010 Doug Zongker
< #
