Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Ruby script to calculate a coefficent of determination for simple 2D regression.
#! /usr/local/bin/ruby
#*********************************************************
# Ruby script to calculate a coefficient of determination.
#*********************************************************
#
class CoefficientOfDetermination
# 説明変数と目的変数
X = [83, 71, 64, 69, 69, 64, 68, 59, 81, 91, 57, 65, 58, 62]
Y = [183, 168, 171, 178, 176, 172, 165, 158, 183, 182, 163, 175, 164, 175]
# Execution
def exec
puts "説明変数 X = {#{X.join(', ')}}"
puts "目的変数 Y = {#{Y.join(', ')}}"
puts "---"
# 単回帰曲線算出
reg_curve = X.reg_curve(Y)
puts " a = %20.16f" % reg_curve[:a]
puts " b = %20.16f" % reg_curve[:b]
puts " c = %20.16f" % reg_curve[:c]
# 推定値
y_e = calc_estimations(X, reg_curve[:a], reg_curve[:b], reg_curve[:c])
# 標本値 Y (目的変数)の平均
y_b = Y.inject(0) { |s, a| s += a } / Y.size.to_f
puts "決定係数"
# 解法-1. 決定係数 (= 推定値の変動 / 標本値の変動)
r_2 = calc_s_r(y_b, y_e) / calc_s_y2(y_b, Y)
puts " R2 (1) = %20.16f" % r_2
# 解法-2. 決定係数 (= 1 - 残差の変動 / 標本値の変動)
r_2 = 1.0 - calc_s_e(Y, y_e) / calc_s_y2(y_b, Y)
puts " R2 (2) = %20.16f" % r_2
rescue => e
$stderr.puts "[#{e.class}] #{e.message}"
e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
exit 1
end
private
# 推定値
#
# @param xs: 説明変数配列
# @param a: 回帰曲線の a
# @param b: 回帰曲線の b
# @param b: 回帰曲線の c
# @return y_e: 推定値配列
def calc_estimations(xs, a, b, c)
y_e = Array.new
begin
xs.each { |x| y_e << a + b * x + c * x * x }
return y_e
rescue => e
raise
end
end
# 推定値の変動
#
# @param y_b: 標本値(目的変数)の平均
# @param y_e: 推定値配列
# @return s_r: 推定値の変動
def calc_s_r(y_b, y_e)
s_r = 0.0
begin
y_e.each do |a|
v = a - y_b
s_r += v * v
end
return s_r
rescue => e
raise
end
end
# 標本値の変動
#
# @param y_b: 標本値(目的変数)の平均
# @param y_s: 標本値(目的変数)配列
# @return s_y2: 標本値の変動
def calc_s_y2(y_b, y_s)
s_y2 = 0.0
begin
y_s.each do |a|
v = a - y_b
s_y2 += v * v
end
return s_y2
rescue => e
raise
end
end
# 残差の変動
#
# @param y_s: 標本値(目的変数)配列
# @param y_e: 推定値配列
# @return s_e: 残差の変動
def calc_s_e(y_s, y_e)
s_e = 0.0
begin
y_s.zip(y_e).each do |a, b|
v = a - b
s_e += v * v
end
return s_e
rescue => e
raise
end
end
end
class Array
# 単回帰曲線(2次)
def reg_curve(y)
# 以下の場合は例外スロー
# - 引数の配列が Array クラスでない
# - 自身配列が空
# - 配列サイズが異なれば例外
raise "Argument is not a Array class!" unless y.class == Array
raise "Self array is nil!" if self.size == 0
raise "Argument array size is invalid!" unless self.size == y.size
n = self.size # number of items
m_x = self.sum / n.to_f # avg(X)
m_y = y.sum / n.to_f # avg(Y)
m_x2 = self.map { |x| x ** 2 }.sum / n.to_f # avg(X^2)
m_x3 = self.map { |x| x ** 3 }.sum / n.to_f # avg(X^3)
m_x4 = self.map { |x| x ** 4 }.sum / n.to_f # avg(X^4)
m_xy = self.zip(y).map { |a, b| a * b }.sum / n.to_f # avg(X * Y)
m_x2y = self.zip(y).map { |a, b| a * a * b }.sum / n.to_f # avg(X^2 * Y)
s_xx = m_x2 - m_x * m_x # Sxx
s_xy = m_xy - m_x * m_y # Sxy
s_xx2 = m_x3 - m_x * m_x2 # Sxx2
s_x2x2 = m_x4 - m_x2 * m_x2 # Sx2x2
s_x2y = m_x2y - m_x2 * m_y # Sx2y
b = s_xy * s_x2x2 - s_x2y * s_xx2
b /= s_xx * s_x2x2 - s_xx2 * s_xx2
c = s_x2y * s_xx - s_xy * s_xx2
c /= s_xx * s_x2x2 - s_xx2 * s_xx2
a = m_y - b * m_x - c * m_x2
{a: a, b: b, c: c}
end
end
CoefficientOfDetermination.new.exec if __FILE__ == $0
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.