1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
# frozen_string_literal: true

# Port of https://github.com/philihp/openskill.js/blob/main/src/statistics.ts

module OpenSkill
  module Statistics
    # Standard normal PDF
    #: (Float x) -> Float
    def self.phi_minor(x)
      Math.exp(-0.5 * x * x) / Math.sqrt(2 * Math::PI)
    end

    # Standard normal CDF via Math.erf
    #: (Float x) -> Float
    def self.phi_major(x)
      0.5 * (1.0 + Math.erf(x / Math.sqrt(2)))
    end

    #: (Float x, Float t) -> Float
    def self.v(x, t)
      xt = x - t
      denom = phi_major(xt)
      return -xt if denom < Float::EPSILON

      phi_minor(xt) / denom
    end

    #: (Float x, Float t) -> Float
    def self.w(x, t)
      xt = x - t
      denom = phi_major(xt)
      if denom < Float::EPSILON
        return x.negative? ? 1.0 : 0.0
      end

      vxt = v(x, t)
      vxt * (vxt + xt)
    end

    #: (Float x, Float t) -> Float
    def self.vt(x, t)
      xx = x.abs
      b = phi_major(t - xx) - phi_major(-t - xx)
      if b < 1e-5
        return x.negative? ? -x - t : -x + t
      end

      a = phi_minor(-t - xx) - phi_minor(t - xx)
      (x.negative? ? -a : a) / b
    end

    #: (Float x, Float t) -> Float
    def self.wt(x, t)
      xx = x.abs
      b = phi_major(t - xx) - phi_major(-t - xx)
      return 1.0 if b < Float::EPSILON

      ((t - xx) * phi_minor(t - xx) + (t + xx) * phi_minor(-t - xx)) / b + vt(x, t)**2
    end
  end
end