Ruby - Array クラス拡張で単回帰曲線（4次回帰モデル）計算！
Updated:
Ruby で Array クラスを拡張して回帰式が4次の単回帰曲線を計算してみました。（連立方程式の解法にはガウスの消去法を使用）
前回は3次回帰モデルについて行なっています。
0. 前提条件Permalink
- LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
- Ruby 2.6.3 での作業を想定。
1. 単回帰曲線（4次回帰モデル）の求め方Permalink
求める曲線を y=a+bx+cx2+dx3+ex4 とすると、残差の二乗和 S は
S=N∑i=1(yi−a−bxi−cx2i−dx3−ex4)2となる。 a,b,c,d,e それぞれで偏微分したものを 0 とする。
∂S∂a=2N∑i=1(a+bxi+cx2i+dx3i+ex4i−yi)=0∂S∂b=2N∑i=1(axi+bx2i+cx3i+dx4i+ex5i−xiyi)=0∂S∂c=2N∑i=1(ax2i+bx3i+cx4i+dx5i+ex6i−x2iyi)=0∂S∂d=2N∑i=1(ax3i+bx4i+cx5i+dx6i+ex7i−x3iyi)=0∂S∂e=2N∑i=1(ax4i+bx5i+cx6i+dx7i+ex8i−x3iyi)=0これらを変形すると、
aN+bN∑i=1xi+cN∑i=1x2i+dN∑i=1x3i+eN∑i=1x4i=N∑i=1yiaN∑i=1xi+bN∑i=1x2i+cN∑i=1x3i+dN∑i=1x4i+eN∑i=1x5i=N∑i=1xiyiaN∑i=1x2i+bN∑i=1x3i+cN∑i=1x4i+dN∑i=1x5i+eN∑i=1x6i=N∑i=1x2iyiaN∑i=1x3i+bN∑i=1x4i+cN∑i=1x5i+dN∑i=1x6i+eN∑i=1x7i=N∑i=1x3iyiaN∑i=1x4i+bN∑i=1x5i+cN∑i=1x6i+dN∑i=1x7i+eN∑i=1x8i=N∑i=1x4iyiとなる。これらの連立方程式を解けばよい。
2. ガウスの消去法による連立方程式の解法についてPermalink
当ブログ過去記事を参照。
3. Ruby スクリプトの作成Permalink
- Shebang ストリング（1行目）では、フルパスでコマンド指定している。（当方の慣習）
File: regression_curve_4d.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
#! /usr/local/bin/ruby
#*********************************************
# Ruby script to calculate a simple regression curve.
# : y = a + b * x + c * x^2 + d * x^3 + e * x^4
# : 連立方程式を ガウスの消去法で解く方法
#*********************************************
#
class Array
def reg_curve_4d(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
sum_x = self.inject(0) { |s, a| s += a }
sum_x2 = self.inject(0) { |s, a| s += a * a }
sum_x3 = self.inject(0) { |s, a| s += a * a * a }
sum_x4 = self.inject(0) { |s, a| s += a * a * a * a }
sum_x5 = self.inject(0) { |s, a| s += a * a * a * a * a }
sum_x6 = self.inject(0) { |s, a| s += a * a * a * a * a * a }
sum_x7 = self.inject(0) { |s, a| s += a * a * a * a * a * a * a }
sum_x8 = self.inject(0) { |s, a| s += a * a * a * a * a * a * a * a }
sum_y = y.inject(0) { |s, a| s += a }
sum_xy = self.zip(y).inject(0) { |s, a| s += a[0] * a[1] }
sum_x2y = self.zip(y).inject(0) { |s, a| s += a[0] * a[0] * a[1] }
sum_x3y = self.zip(y).inject(0) { |s, a| s += a[0] * a[0] * a[0] * a[1] }
sum_x4y = self.zip(y).inject(0) { |s, a| s += a[0] * a[0] * a[0] * a[0] * a[1] }
mtx = [
[self.size, sum_x, sum_x2, sum_x3, sum_x4, sum_y],
[ sum_x, sum_x2, sum_x3, sum_x4, sum_x5, sum_xy],
[ sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, sum_x2y],
[ sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, sum_x3y],
[ sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, sum_x4y]
]
ans = solve_ge(mtx)
{a: ans[0][-1], b: ans[1][-1], c: ans[2][-1], d: ans[3][-1], e: ans[4][-1]}
end
private
# 連立方程式の解（ガウスの消去法）
def solve_ge(a)
n = a.size
# 前進消去
(n - 1).times do |k|
(k + 1).upto(n - 1) do |i|
d = a[i][k] / a[k][k].to_f
(k + 1).upto(n) do |j|
a[i][j] -= a[k][j] * d
end
end
end
# 後退代入
(n - 1).downto(0) do |i|
d = a[i][n]
(i + 1).upto(n - 1) do |j|
d -= a[i][j] * a[j][n]
end
a[i][n] = d / a[i][i].to_f
end
return a
end
end
# 説明変数と目的変数
#ary_x = [107, 336, 233, 82, 61, 378, 129, 313, 142, 428]
#ary_y = [286, 851, 589, 389, 158, 1037, 463, 563, 372, 1020]
ary_x = [83, 71, 64, 69, 69, 64, 68, 59, 81, 91, 57, 65, 58, 62]
ary_y = [183, 168, 171, 178, 176, 172, 165, 158, 183, 182, 163, 175, 164, 175]
puts "説明変数 X = {#{ary_x.join(', ')}}"
puts "目的変数 Y = {#{ary_y.join(', ')}}"
puts "---"
# 単回帰曲線算出
reg_line = ary_x.reg_curve_4d(ary_y)
puts "a = #{reg_line[:a]}"
puts "b = #{reg_line[:b]}"
puts "c = #{reg_line[:c]}"
puts "d = #{reg_line[:d]}"
puts "e = #{reg_line[:e]}"
4. Ruby スクリプトの実行Permalink
$ ./regression_curve_4d.rb
説明変数 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}
---
a = -8069.31709645464
b = 454.222121309131
c = -9.336623651523444
d = 0.08475031013977985
e = -0.0002862803682765933
5. 視覚的な確認Permalink
参考までに、上記スクリプトで使用した2変量の各点と作成された単回帰曲線を gnuplot で描画してみた。
6. その他Permalink
前回は3次回帰モデル、今回は4次回帰モデルについて行なったが、5次以降も同様に行える。（当ブログでは、5次以降については記事にしない）
以上。
Comments