Last active
May 6, 2019 02:11
-
-
Save komasaru/032cc6392395fc9d94b0099f97a82b19 to your computer and use it in GitHub Desktop.
Ruby script to calculate a simple regression curve.(power)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#! /usr/local/bin/ruby | |
#********************************************* | |
# Ruby script to calculate a simple regression curve. | |
# : y = a * x**b | |
# : 連立方程式を ガウスの消去法で解く方法 | |
#********************************************* | |
# | |
class Array | |
def reg_curve_pow(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_lx = self.inject(0) { |s, a| s += Math.log(a) } | |
sum_lx2 = self.inject(0) { |s, a| s += Math.log(a) * Math.log(a) } | |
sum_ly = y.inject(0) { |s, a| s += Math.log(a) } | |
sum_lxly = self.zip(y).inject(0) { |s, a| s += Math.log(a[0]) * Math.log(a[1]) } | |
mtx = [ | |
[self.size, sum_lx, sum_ly], | |
[ sum_lx, sum_lx2, sum_lxly] | |
] | |
ans = solve_ge(mtx) | |
{a: Math.exp(ans[0][-1]), b: ans[1][-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_pow(ary_y) | |
puts "a = #{reg_line[:a]}" | |
puts "b = #{reg_line[:b]}" |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment