Module: CompSci::Fit

Defined in:
lib/compsci/fit.rb

Class Method Summary collapse

Class Method Details

.best(xs, ys) ⇒ Object

Run logarithmic, linear, exponential, and power fits Return the stats for the best fit (highest r^2)

Takes x and y values and returns [a, b, r2, fn]



22
23
24
25
26
27
28
29
30
31
32
33
34
# File 'lib/compsci/fit.rb', line 22

def self.best xs, ys
  vals = []
  max_r2 = 0
  [:logarithmic, :linear, :exponential, :power].each { |fn|
    a, b, r2 = Fit.send(fn, xs, ys)
    # p [a, b, r2, fn]
    if r2 > max_r2
      vals = [a, b, r2, fn]
      max_r2 = r2
    end
  }
  vals
end

.constant(xs, ys) ⇒ Object

Fits the functional form: a (+ 0x)

Takes x and y values and returns [a, variance]



9
10
11
12
13
# File 'lib/compsci/fit.rb', line 9

def self.constant xs, ys
  y_bar = sigma(ys) / ys.size.to_f
  variance = sigma(ys) { |y| (y - y_bar) ** 2 }
  [y_bar, variance]
end

.error(xys, &blk) ⇒ Object

Takes an array of x/y pairs and calculates the general R^2 value to measure fit against a predictive function, which is the block supplied to error:

e.g. error(xys) { |x| 5 + 2 * x }

See: en.wikipedia.org/wiki/Coefficient_of_determination



62
63
64
65
66
67
68
# File 'lib/compsci/fit.rb', line 62

def self.error xys, &blk
  y_bar  = sigma(xys) { |_, y| y                   } / xys.size.to_f
  ss_tot = sigma(xys) { |_, y| (y - y_bar)    ** 2 }
  ss_res = sigma(xys) { |x, y| (yield(x) - y) ** 2 }

  1 - (ss_res / ss_tot)
end

.exponential(xs, ys) ⇒ Object

To fit a functional form: y = ae^(bx).

Takes x and y values and returns [a, b, r^2].

See: mathworld.wolfram.com/LeastSquaresFittingExponential.html



121
122
123
124
125
126
127
128
129
130
131
132
133
134
# File 'lib/compsci/fit.rb', line 121

def self.exponential xs, ys
  n     = xs.size
  xys   = xs.zip(ys)
  sxlny = sigma(xys) { |x, y| x * Math.log(y) }
  slny  = sigma(xys) { |_, y| Math.log(y)     }
  sx2   = sigma(xys) { |x, _| x * x           }
  sx    = sigma xs

  c = n * sx2 - sx ** 2
  a = (slny * sx2 - sx * sxlny) / c
  b = ( n * sxlny - sx * slny ) / c

  return Math.exp(a), b, self.error(xys) { |x| Math.exp(a + b * x) }
end

.linear(xs, ys) ⇒ Object

Fits the functional form: a + bx.

Takes x and y values and returns [a, b, r^2].

See: mathworld.wolfram.com/LeastSquaresFitting.html



99
100
101
102
103
104
105
106
107
108
109
110
111
112
# File 'lib/compsci/fit.rb', line 99

def self.linear xs, ys
  n   = xs.size
  xys = xs.zip(ys)
  sx  = sigma xs
  sy  = sigma ys
  sx2 = sigma(xs)  { |x|   x ** 2 }
  sxy = sigma(xys) { |x, y| x * y  }

  c = n * sx2 - sx**2
  a = (sy * sx2 - sx * sxy) / c
  b = ( n * sxy - sx * sy ) / c

  return a, b, self.error(xys) { |x| a + b * x }
end

.logarithmic(xs, ys) ⇒ Object

To fit a functional form: y = a + b*ln(x).

Takes x and y values and returns [a, b, r^2].

See: mathworld.wolfram.com/LeastSquaresFittingLogarithmic.html



77
78
79
80
81
82
83
84
85
86
87
88
89
90
# File 'lib/compsci/fit.rb', line 77

def self.logarithmic xs, ys
  n     = xs.size
  xys   = xs.zip(ys)
  slnx2 = sigma(xys) { |x, _| Math.log(x) ** 2 }
  slnx  = sigma(xys) { |x, _| Math.log(x)      }
  sylnx = sigma(xys) { |x, y| y * Math.log(x)  }
  sy    = sigma(xys) { |_, y| y                }

  c = n * slnx2 - slnx ** 2
  b = ( n * sylnx - sy * slnx ) / c
  a = (sy - b * slnx) / n

  return a, b, self.error(xys) { |x| a + b * Math.log(x) }
end

.power(xs, ys) ⇒ Object

To fit a functional form: y = ax^b.

Takes x and y values and returns [a, b, r^2].

See: mathworld.wolfram.com/LeastSquaresFittingPowerLaw.html



143
144
145
146
147
148
149
150
151
152
153
154
155
# File 'lib/compsci/fit.rb', line 143

def self.power xs, ys
  n       = xs.size
  xys     = xs.zip(ys)
  slnxlny = sigma(xys) { |x, y| Math.log(x) * Math.log(y) }
  slnx    = sigma(xs)  { |x   | Math.log(x)               }
  slny    = sigma(ys)  { |   y| Math.log(y)               }
  slnx2   = sigma(xs)  { |x   | Math.log(x) ** 2          }

  b = (n * slnxlny - slnx * slny) / (n * slnx2 - slnx ** 2)
  a = (slny - b * slnx) / n

  return Math.exp(a), b, self.error(xys) { |x| (Math.exp(a) * (x ** b)) }
end

.sigma(enum, &block) ⇒ Object

Enumerates over enum mapping block if given, returning the sum of the result. Eg:

sigma([1, 2, 3])                # => 1 + 2 + 3 => 6
sigma([1, 2, 3]) { |n| n ** 2 } # => 1 + 4 + 9 => 14


47
48
49
50
# File 'lib/compsci/fit.rb', line 47

def self.sigma enum, &block
  enum = enum.map(&block) if block
  enum.inject { |sum, n| sum + n }
end