Skip to content

Instantly share code, notes, and snippets.

@komasaru
Last active June 20, 2022 04:38
Show Gist options
  • Save komasaru/105de1f29556d4b998662a591835bd29 to your computer and use it in GitHub Desktop.
Save komasaru/105de1f29556d4b998662a591835bd29 to your computer and use it in GitHub Desktop.
Ruby script to calculate a simple regression curve.(5d)
#! /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 + f * x^5
# : 連立方程式を ガウスの消去法(ピボット選択)で解く方法
#*********************************************
#
class Array
def reg_curve_5d(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_x9 = self.inject(0) { |s, a| s += a * a * a * a * a * a * a * a * a }
sum_x10 = self.inject(0) { |s, a| s += a * a * 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]
}
sum_x5y = self.zip(y).inject(0) {
|s, a| s += a[0] * a[0] * a[0] * a[0] * a[0] * a[1]
}
mtx = [
[self.size, sum_x, sum_x2, sum_x3, sum_x4, sum_x5, sum_y],
[ sum_x, sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, sum_xy],
[ sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, sum_x2y],
[ sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, sum_x3y],
[ sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, sum_x4y],
[ sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, sum_x10, sum_x5y]
]
ans = solve_ge(mtx)
{
a: ans[0][-1], b: ans[1][-1], c: ans[2][-1],
d: ans[3][-1], e: ans[4][-1], f: ans[5][-1]
}
end
private
# ピボット選択
def pivot(a, k)
pv = k
v_max = a[k][k].abs
(k + 1).upto(a.size - 1) do |i|
if a[i][k].abs > v_max
pv = i
v_max = a[i][k].abs
end
end
if k != pv
dummy = a[k]
a[k] = a[pv]
a[pv] = dummy
end
return a
rescue => e
raise
end
# 連立方程式の解(ガウスの消去法(ピボット選択))
def solve_ge(a)
n = a.size
# 前進消去
(n - 1).times do |k|
a = pivot(a, 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_5d(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]}"
puts "f = #{reg_line[:f]}"
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment