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
|
