Ruby - Array クラスを拡張して重回帰分析（2次多項式モデル）！
Updated:
過去に、説明（独立）変数2個、目的（従属）変数１個の「重回帰式」の計算を Ruby の Array クラスを拡張する方法で実装しました。
今回は、重回帰式を2次多項式にしてみました。
0. 前提条件Permalink
- LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
- Ruby 2.6.4 での作業を想定。
1. 重回帰式（2次多項式モデル）の求め方Permalink
求める重回帰式を y=b0+b1x1+b2x2+b3x1x2+b4x12+b5x22 （説明変数が2個）とする場合、 x3=x1x2, x4=x12, x5=x22 と置くと、 y=b0+b1x1+b2x2+b3x3+b4x4+b5x5 （説明変数が5個）となるので、残差の二乗和 S は
S=∑(yi−b0−b1x1i−b2x2i−b3x3i−b4x4i−b5x5i)2となる。 (∑はN∑i=1)
b0,⋯,b5 それぞれで偏微分したものを 0 とする。
∂S∂b0=2∑(b0+b1x1i+b2x2i+b3x3i+b4x4i+b5x5i−yi)=0∂S∂b1=2∑(b0x1i+b1x1i2+b2x1ix2i+b3x1ix3i+b4x1ix4i+b5x1ix5i−x1iyi)=0∂S∂b2=2∑(b0x2i+b1x1ix2i+b2x2i2+b3x2ix3i+b4x2ix4i+b5x2ix5i−x2iyi)=0∂S∂b3=2∑(b0x3i+b1x1ix3i+b2x2ix3i+b3x3i2+b4x3ix4i+b5x3ix5i−x3iyi)=0∂S∂b4=2∑(b0x4i+b1x1ix4i+b2x2ix4i+b3x3ix4i+b4x4i2+b5x4ix5i−x4iyi)=0∂S∂b5=2∑(b0x5i+b1x1ix5i+b2x2ix5i+b3x3ix5i+b4x4ix5i+b5x5i2−x5iyi)=0これらを変形すると、
b0N+b1∑x1i+b2∑x2i+b3∑x3i+b4∑x4i+b5∑x5i=∑yib0∑x1i+b1∑x1i2+b2∑x1ix2i+b3∑x1ix3i+b4∑x1ix4i+b5∑x1ix5i=∑x1iyib0∑x2i+b1∑x1ix2i+b2∑x2i2+b3∑x2ix3i+b4∑x2ix4i+b5∑x2ix5i=∑x2iyib0∑x3i+b1∑x1ix3i+b2∑x2ix3i+b3∑x3i2+b4∑x3ix4i+b5∑x3ix5i=∑x3iyib0∑x4i+b1∑x1ix4i+b2∑x2ix4i+b3∑x3ix4i+b4∑x4i2+b5∑x4ix5i=∑x4iyib0∑x5i+b1∑x1ix5i+b2∑x2ix5i+b3∑x3ix5i+b4∑x4ix5i+b5∑x5i2=∑x5iyi(但し、x3=x1x2, x4=x12, x5=x22)となる。これらの b0,⋯,b5 の連立1次方程式を解けばよい。
※x3=x1x2, x4=x12, x5=x22 と置かず、直接計算してもよいが、偏微分や連立方程式が煩雑になり、分かりにくくなる。
2. ガウスの消去法による連立方程式の解法についてPermalink
当ブログ過去記事を参照。
3. Ruby スクリプトの作成Permalink
- Shebang ストリング（1行目）では、フルパスでコマンド指定している。（当方の慣習）
File: regression_multi_2d.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
#! /usr/local/bin/ruby
#*********************************************
# Ruby script to calculate a multiple regression function(2D).
# * y = b0 + b1x1 + b2x2 + b3x1x2 + b4x1^2 + b5x2^2
# これは、
# y = b0 + b1x1 + b2x2 + b3x3 + b4x4 + b5x5
# (但し、 x3 = x1x2, x4 = x1^2, x5 = x2^2)
# と同じ。
#*********************************************
#
class Array
def reg_multi_2d(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
x3 = x1.zip(x2).map { |a, b| a * b }
x4 = x1.map { |a| a * a }
x5 = x2.map { |a| a * a }
mtx = Array.new(6).map { Array.new(7, 0.0) }
# 左辺・対角成分
mtx[0][0] = s_size
mtx[1][1] = x1.map { |a| a * a }.sum
mtx[2][2] = x2.map { |a| a * a }.sum
mtx[3][3] = x3.map { |a| a * a }.sum
mtx[4][4] = x4.map { |a| a * a }.sum
mtx[5][5] = x5.map { |a| a * a }.sum
# 左辺・右上成分
mtx[0][1] = x1.sum
mtx[0][2] = x2.sum
mtx[0][3] = x3.sum
mtx[0][4] = x4.sum
mtx[0][5] = x5.sum
mtx[1][2] = x1.zip(x2).map { |a, b| a * b }.sum
mtx[1][3] = x1.zip(x3).map { |a, b| a * b }.sum
mtx[1][4] = x1.zip(x4).map { |a, b| a * b }.sum
mtx[1][5] = x1.zip(x5).map { |a, b| a * b }.sum
mtx[2][3] = x2.zip(x3).map { |a, b| a * b }.sum
mtx[2][4] = x2.zip(x4).map { |a, b| a * b }.sum
mtx[2][5] = x2.zip(x5).map { |a, b| a * b }.sum
mtx[3][4] = x3.zip(x4).map { |a, b| a * b }.sum
mtx[3][5] = x3.zip(x5).map { |a, b| a * b }.sum
mtx[4][5] = x4.zip(x5).map { |a, b| a * b }.sum
# 左辺・左下成分
mtx[1][0] = mtx[0][1]
mtx[2][0] = mtx[0][2]
mtx[2][1] = mtx[1][2]
mtx[3][0] = mtx[0][3]
mtx[3][1] = mtx[1][3]
mtx[3][2] = mtx[2][3]
mtx[4][0] = mtx[0][4]
mtx[4][1] = mtx[1][4]
mtx[4][2] = mtx[2][4]
mtx[4][3] = mtx[3][4]
mtx[5][0] = mtx[0][5]
mtx[5][1] = mtx[1][5]
mtx[5][2] = mtx[2][5]
mtx[5][3] = mtx[3][5]
mtx[5][4] = mtx[4][5]
# 右辺
mtx[0][6] = y.sum
mtx[1][6] = x1.zip(y).map { |a, b| a * b }.sum
mtx[2][6] = x2.zip(y).map { |a, b| a * b }.sum
mtx[3][6] = x3.zip(y).map { |a, b| a * b }.sum
mtx[4][6] = x4.zip(y).map { |a, b| a * b }.sum
mtx[5][6] = x5.zip(y).map { |a, b| a * b }.sum
# 連立方程式を解く (ガウスの消去法)
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
# 説明(独立)変数と目的(従属)変数
ary_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]
]
ary_y = [45, 38, 41, 34, 59, 47, 35, 43, 54, 52]
# 重回帰式算出(b0, b1, b2, ...)
reg_multi = ary_x.reg_multi_2d(ary_y)
# 結果出力
ary_x.each_with_index do |x, i|
puts "説明変数 X#{i + 1} = {#{ary_x[i].join(', ')}}"
end
puts "目的変数 Y = {#{ary_y.join(', ')}}"
puts "---"
p reg_multi
4. Ruby スクリプトの実行Permalink
$ ./regression_multi_2d.rb
説明変数 X1 = {17.5, 17.0, 18.5, 16.0, 19.0, 19.5, 16.0, 18.0, 19.0, 19.5}
説明変数 X2 = {30, 25, 20, 30, 45, 35, 25, 35, 35, 40}
目的変数 Y = {45, 38, 41, 34, 59, 47, 35, 43, 54, 52}
---
[-333.9032985981502, 45.8844170002146, -4.176065472613139, 0.21073632980137985, -1.384555099988024, 0.013710266583215294]
5. 視覚的な確認Permalink
参考までに、上記スクリプトで使用した2変量の各点と作成された単回帰曲線を gnuplot で描画してみた。
以上。
Comments