mk-mode.com

mk-mode.com

Linux, Debian, IT, Server, PG, Ruby, Rails, Python, C++, Fortran, PC, MariaDB, math, GIS, etc...

Ruby - 単回帰分析（線形回帰）の決定係数計算！

Updated:


Ruby で2つの単回帰分析（線形回帰; 単回帰直線）の決定係数を計算してみました。

単回帰直線や相関係数の計算は Array クラスを拡張して行なっています。

0. 前提条件

  • LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
  • Ruby 2.6.3 での作業を想定。

1. 決定係数について

回帰分析において、目的変数の標本値（実測値）に対する目的変数の推測値（予測値）の説明力を表す指標（言い換えれば、説明変数（独立変数）が目的変数（従属変数）をどれくらい説明できているかを表す統計量）が

\[\begin{eqnarray*} 決定係数 R^2 \end{eqnarray*}\]

である。（「\(R\)の\(2\)乗」で表現するが、必ずしも何かの値の \(2\) 乗になるという意味ではない）

決定係数 \(R^2\) は次のように定義する。（定義の仕方は複数あるが、次の定義が最も一般的）

\[\begin{eqnarray*} 決定係数 R^2 = \frac{推定値の変動}{標本値の変動} = \frac{S_R}{S_{y^2}} \end{eqnarray*}\]

但し、

\[\begin{eqnarray*} 標本値の変動 &=& \sum_{i=1}^{N}(y_i - \bar{y})^2 = S_{y^2} \\ 推定値の変動 &=& \sum_{i=1}^{N}(Y_i - \bar{y})^2 = S_R \\ 残差の変動 &=& \sum_{i=1}^{N}(y_i - Y_i)^2 = S_E \end{eqnarray*}\]

これら3つの変動の間には次のような関係が成り立つ。

\[\begin{eqnarray*} 標本値の変動 = 推定値の変動 + 残差の変動 \end{eqnarray*}\]

これから、

\[\begin{eqnarray*} 1 &=& \frac{推定値の変動}{標本値の変動} + \frac{残差の変動}{標本値の変動} \\ \therefore \ \ 1 &=& 決定係数 R^2 + \frac{残差の変動}{標本値の変動} \\ \end{eqnarray*}\]

よって、

\[\begin{eqnarray*} 決定係数 R^2 &=& 1 - \frac{残差の変動}{標本値の変動} \end{eqnarray*}\]

となる。これは、「残差の変動が \(0\) に近ければ、決定係数が \(1\) が近くなる」ということで、「決定係数 \(R^2\) が \(1\) に近いほど、当てはまりが良い（説明変数が目的変数をより説明できている）」と表現できることになる。

ちなみに、単回帰分析（線形回帰）の場合、

\[\begin{eqnarray*} (xとyの相関係数)^2 = 決定係数 \end{eqnarray*}\]

となる。また、

\[\begin{eqnarray*} S_{y^2} &=& \sum_{i=1}^{N}y_i^2 - \frac{\left( \displaystyle{\sum_{i=1}^{N}y_i} \right)^2}{N} \\ S_{xy} &=& \sum_{i=1}^{N}x_{i}y_{i} - \frac{\displaystyle{\sum_{i=1}^{N}x_i} \sum_{i=1}^{N}y_i}{N} \\ S_R &=& b \cdot S_{xy} \\ \end{eqnarray*}\]

でるあること（導出方法は略）を利用して、 \(\begin{eqnarray*} 決定係数 R^2 &=& \frac{S_R}{S_{y^2}} \end{eqnarray*}\)

を求めることもできる。

2. Ruby スクリプトの作成

  • 以下のスクリプトでは4種の方法で決定係数を計算している。
  • Shebang ストリング（1行目）では、フルパスでコマンド指定している。（当方の慣習）

File: coefficient_of_determination_line.rb

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196

#! /usr/local/bin/ruby
#*********************************************************
# Ruby script to calculate a coefficient of determination.
#*********************************************************
#
class CoefficientOfDetermination
  # 説明変数と目的変数
  X = [107, 336, 233,  82,  61,  378, 129, 313, 142,  428]
  Y = [286, 851, 589, 389, 158, 1037, 463, 563, 372, 1020]

  # Execution
  def exec
    puts "説明変数 X = {#{X.join(', ')}}"
    puts "目的変数 Y = {#{Y.join(', ')}}"
    puts "---"
    # 単回帰直線算出（切片と傾き）
    reg_line = X.reg_line(Y)
    puts "    切片  a = %20.16f" % reg_line[:intercept]
    puts "    傾き  b = %20.16f" % reg_line[:slope]
    # 相関係数
    r = X.r(Y)
    puts "相関係数  r = %20.16f" % r
    # 推定値
    y_e = calc_estimations(X, Y, reg_line[:intercept], reg_line[:slope])
    # 標本値 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
    # 解法-3. 決定係数 (公式使用(解法-1,2の変形))
    r_2 = calc_r_2(X, Y, reg_line[:slope])
    puts "     R2 (3) = %20.16f" % r_2
    # 解法-4. 決定係数 (= 相関係数の2乗)
    r_2 = r * r
    puts "     R2 (4) = %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    y: 目的変数配列
  # @param    a: 回帰直線の切片
  # @param    b: 回帰直線の傾き
  # @return y_e: 推定値配列
  def calc_estimations(xs, y, a, b)
    y_e = Array.new

    begin
      xs.each { |x| y_e << a + b * 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

  # 決定係数 (公式使用)
  #
  # @param    x: 説明変数配列
  # @param    y: 目的変数配列
  # @param    b: 回帰直線の傾き
  # @return r_2: 残差の変動
  def calc_r_2(x, y, b)
    n = x.size
    sum_x  = x.inject(0.0) { |s, a| s += a }
    sum_y  = y.inject(0.0) { |s, a| s += a }
    sum_y2 = y.inject(0.0) { |s, a| s += a * a }
    sum_xy = x.zip(y).inject(0.0) { |s, a| s += a[0] * a[1] }
    s_y2 = sum_y2 - sum_y * sum_y / n.to_f
    s_xy = sum_xy - sum_x * sum_y / n.to_f
    s_r  = b * s_xy
    return s_r / s_y2
  rescue => e
    raise
  end
end

class Array
  # 単回帰直線
  def reg_line(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

    # x の総和
    sum_x = self.inject(0) { |s, a| s += a }
    # y の総和
    sum_y = y.inject(0) { |s, a| s += a }
    # x^2 の総和
    sum_xx = self.inject(0) { |s, a| s += a * a }
    # x * y の総和
    sum_xy = self.zip(y).inject(0) { |s, a| s += a[0] * a[1] }
    # 切片 a
    a  = sum_xx * sum_y - sum_xy * sum_x
    a /= (self.size * sum_xx - sum_x * sum_x).to_f
    # 傾き b
    b  = self.size * sum_xy - sum_x * sum_y
    b /= (self.size * sum_xx - sum_x * sum_x).to_f
    {intercept: a, slope: b}
  end

  # 相関係数
  def r(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

    # x の相加平均, y の相加平均 (arithmetic mean)
    mean_x = self.inject(0) { |s, a| s += a } / self.size.to_f
    mean_y = y.inject(0) { |s, a| s += a } / y.size.to_f
    # x と y の共分散の分子 (covariance)
    cov = self.zip(y).inject(0) { |s, a| s += (a[0] - mean_x) * (a[1] - mean_y) }
    # x の分散の分子, y の分散の分子 (variance)
    var_x = self.inject(0) { |s, a| s += (a - mean_x) ** 2 }
    var_y = y.inject(0) { |s, a| s += (a - mean_y) ** 2 }
    # 相関係数 (correlation coefficient)
    r = cov / Math.sqrt(var_x)
    r /= Math.sqrt(var_y)
  end
end

CoefficientOfDetermination.new.exec if __FILE__ == $0

4. Ruby スクリプトの実行

$ ./coefficient_of_determination_line.rb
説明変数 X = {107, 336, 233, 82, 61, 378, 129, 313, 142, 428}
目的変数 Y = {286, 851, 589, 389, 158, 1037, 463, 563, 372, 1020}
---
    切片  a =  99.0747587724579120
    傾き  b =   2.1445235003510281
相関係数  r =   0.9451950086576620
決定係数
     R2 (1) =   0.8933936043913584
     R2 (2) =   0.8933936043913584
     R2 (3) =   0.8933936043913585
     R2 (4) =   0.8933936043913578

決定係数が 約0.9 となっているので、単回帰直線の当てはまりは良いと言える。

5. 視覚的な確認

参考までに、上記スクリプトで使用した2変量の各点と作成された単回帰直線を gnuplot で描画してみた。

REGRESSION_LINE

以上。





 

Sponsored Link

 

Comments

You May Also Enjoy

Debian 12 (bookworm) - Ruby 3.3 インストール（ソースビルド）！  (2024-05-18)

Debian GNU/Linux 12 (bookworm) に Ruby 3.3 をソースをビルドしてインストール方法についての記録です。 以前古いバージョンでの作業時に残していた記録を参考に作業を行い、今回更新した作業記録を貼付する形式の内容となっています。 （当然ながら、興味がなければスルーしてください）