Ruby script to calculate an adjusted coefficient of determination.

adjusted_coefficient_of_determination.rb

#! /usr/local/bin/ruby
#*******************************************************************
# Ruby script to calculate an adjusted coefficient of determination.
#*******************************************************************
#
class AdjustedCoefficientOfDetermination
  # 説明変数と目的変数
  X = [
    [17.5, 17.0, 18.5, 16.0, 19.0, 19.5, 16.0, 18.0, 19.0, 19.5],
    [30  , 25  , 20  , 30  , 45  , 35  , 25  , 35  , 35  , 40  ]
  ]
  Y  = [45, 38, 41, 34, 59, 47, 35, 43, 54, 52]

  def exec
    puts "説明変数 X1 = [#{X[0].join(', ')}]"
    puts "説明変数 X2 = [#{X[1].join(', ')}]"
    puts "目的変数 Y  = [#{Y.join(', ')}]"
    puts "---"
    # 重回帰式算出
    reg = X.reg_multi(Y)
    puts "    a = %20.16f" % reg[0]
    puts "    b = %20.16f" % reg[1]
    puts "    c = %20.16f" % reg[2]
    # 推定値配列
    y_e = calc_estimations(X, Y, reg)
    # 標本値 Y （目的変数）の平均
    y_b = Y.inject(0) { |s, a| s += a } / Y.size.to_f
    # 自由度調整済み決定係数
    p, n = X.size, Y.size
    r_2  = calc_s_e(Y, y_e) / (n - p - 1)
    r_2 /= calc_s_r(Y, y_b) / (n - 1)
    r_2  = 1 - r_2
    puts "---"
    puts "自由度調整済み決定係数 R2f = %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: 説明変数配列(2個)
  # @param    y: 目的変数配列
  # @param  reg: 重回帰式の値 a, b, c
  # @return y_e: 推定値配列
  def calc_estimations(xs, y, reg)
    y_e = Array.new

    begin
      a, b, c = reg
      y.zip(xs[0], xs[1]).each do |y, x1, x2|
        y_e << a + b * x1 + c * x2
      end
      return y_e
    rescue => e
      raise
    end
  end

  # 標本値の変動
  #
  # @param  y_s: 標本値（目的変数）配列
  # @param  y_b: 標本値（目的変数）の平均
  # @return s_r: 標本値の変動
  def calc_s_r(y_s, y_b)
    s_r = 0.0

    begin
      y_s.each do |a|
        v = a - y_b
        s_r += v * v
      end
      return s_r
    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
  def reg_multi(y)
    # 元の数、変量内のサンプル数
    e_size, s_size = self.size, y.size
    # 以下の場合は例外スロー
    # - 引数の配列が Array クラスでない
    # - 自身配列が空
    # - 配列サイズが異なれば例外
    raise "Argument is not a Array class!"    unless y.class == Array
    raise "Self array is nil!"                if     e_size == 0
    raise "Argument array size is invalid!"   unless self[0].size == s_size
    1.upto(e_size - 1) do |i|
      raise "Argument array size is invalid!" unless self[0].size == self[i].size
    end
    x1, x2 = self
    sum_x1   = x1.sum
    sum_x2   = x2.sum
    sum_y    = y.sum
    sum_x1x1 = x1.map { |a| a * a }.sum
    sum_x1x2 = x1.zip(x2).map { |a, b| a * b }.sum
    sum_x1y  = x1.zip(y).map { |a, b| a * b }.sum
    sum_x2x1 = sum_x1x2
    sum_x2x2 = x2.map { |a| a * a }.sum
    sum_x2y  = x2.zip(y).map { |a, b| a * b }.sum
    mtx = [
      [s_size, sum_x1  , sum_x2  , sum_y  ],
      [sum_x1, sum_x1x1, sum_x1x2, sum_x1y],
      [sum_x2, sum_x2x1, sum_x2x2, sum_x2y]
    ]
    # 連立方程式を解く (ガウスの消去法)
    return gauss_e(mtx)
  end

private
  # ガウスの消去法
  def gauss_e(ary)
    # 行数
    n = ary.size
    # 前進消去
    0.upto(n - 2) do |k|
      (k + 1).upto(n - 1) do |i|
        if ary[k][k] == 0
          puts "解けない！"
          exit 1
        end
        d = ary[i][k] / ary[k][k].to_f
        (k + 1).upto(n) do |j|
          ary[i][j] -= ary[k][j] * d
        end
      end
    end
    # 後退代入
    (n - 1).downto(0) do |i|
      if ary[i][i] == 0
        puts "解けない！"
        exit 1
      end
      d = ary[i][n]
      (i + 1).upto(n - 1) do |j|
        d -= ary[i][j] * ary[j][n]
      end
      ary[i][n] = d / ary[i][i].to_f
    end
    return ary.map { |a| a[-1] }
  end
end

AdjustedCoefficientOfDetermination.new.exec if __FILE__ == $0