diff options
| -rw-r--r-- | lin_reg.py | 50 | ||||
| -rw-r--r-- | lin_reg_example.py | 6 | ||||
| -rw-r--r-- | sterling_approx.awk | 18 |
3 files changed, 15 insertions, 59 deletions
diff --git a/lin_reg.py b/lin_reg.py deleted file mode 100644 index 8fb88b2..0000000 --- a/lin_reg.py +++ /dev/null @@ -1,50 +0,0 @@ -#!/usr/bin/python -#/usr/local/bin/python2.7 - -import numpy as np - - -def main(): - - dataz = np.genfromtxt('example.csv', delimiter=',', names=True) - - y = dataz['attitude'] - x = dataz['correct'] - print 'x = {} {}'.format(x, np.mean(x)) - print 'y = {} {}'.format(y, np.mean(y)) - - delta_x = x - np.mean(x) - delta_y = y - np.mean(y) - delta_xy = delta_x*delta_y - delta2_x = delta_x**2 - delta2_y = delta_y**2 - #print 'delta_x = {} {}'.format(delta_x, np.sum(delta_x)) - #print 'delta_y = {} {}'.format(delta_y, np.sum(delta_y)) - #print 'delta_xy = {} {}'.format(delta_xy, np.sum(delta_xy)) - #print 'delta2_x = {} {}'.format(delta2_x, np.sum(delta2_x)) - #print 'delta2_y = {} {}'.format(delta2_y, np.sum(delta2_y)) - - sx = np.sqrt(np.sum(delta2_x)/(delta2_x.size - 1)) - sy = np.sqrt(np.sum(delta2_y)/(delta2_y.size - 1)) - r = np.sum(delta_xy)/np.sqrt(np.sum(delta2_x)*np.sum(delta2_y)) - #print 'sx = {}'.format(sx) - #print 'sy = {}'.format(sy) - print 'r = {}'.format(r) - print 'r^2 = {}'.format(r**2) - - b = r*(sy/sx) - a = np.mean(y) - b*np.mean(x) - print 'b = {}'.format(b) - print 'a = {}'.format(a) - - err = y - (a + b*x) - b_err = np.sqrt(np.sum(err**2.0)/((x.size - 2)*np.sum(delta2_x))) - a_err = np.sqrt(np.sum(x**2.0)*np.sum(err**2.0)/(x.size*(x.size - 2)*np.sum(delta2_x))) - print 'b_err = {}'.format(b_err) - print 'a_err = {}'.format(a_err) - - print ' (attitude) \t= {}(correct) + {}'.format(b, a) - - -if __name__ == '__main__': - main() diff --git a/lin_reg_example.py b/lin_reg_example.py index 3a2043e..47f1a84 100644 --- a/lin_reg_example.py +++ b/lin_reg_example.py @@ -1,12 +1,11 @@ -#!/usr/bin/python -#/usr/local/bin/python2.7 +#!/usr/bin/env python import numpy as np def main(): - dataz = np.genfromtxt('example.csv', delimiter=',', names=True) + dataz = np.genfromtxt('data/example.csv', delimiter=',', names=True) y = dataz['attitude'] x = dataz['correct'] @@ -31,6 +30,7 @@ def main(): print 'Sy = {}'.format(Sy) print 'r = {}'.format(r) print 'r^2 = {}'.format(r**2) + print 'r^2 dB = {}'.format(10.0*np.log10(r**2)) b = r*(Sy/Sx) a = np.mean(y) - b*np.mean(x) diff --git a/sterling_approx.awk b/sterling_approx.awk index 5abe41d..ca119ac 100644 --- a/sterling_approx.awk +++ b/sterling_approx.awk @@ -1,17 +1,23 @@ #!/usr/bin/awk -f + +# https://en.wikipedia.org/wiki/Sterling_Approximation +# An alternative approximation for the Gamma function stated by Srinivasa +# Ramanujan (Ramanujan 1988) is +# Gamma(1+x) ~= sqrt(pi)((x/e)^x)(8x^3 + 4x^2 + x + 1/30)^(1/6) +# for x >= 0. The equivalent approximation for ln(n!) has an asymptotic error +# of 1/(1400*n^3) ... + + ### sterling_approx.awk # https://en.wikipedia.org/wiki/Stirling%27s_approximation BEGIN { ARGV[1] ? n = ARGV[1] : n = 0 pi = 4*atan2(1,1) - p = 0 + f = 0 if (n > 0) { - p = 1 - for (m=n; m>0; m--) - p *= n*exp(-1) - p = sqrt(2*pi*n)*p*(1 + 1/(12*n) + 1/(288*n*n) - 139/(51840*n*n*n) - 571/(2488320*n*n*n*n)) + f = sqrt(2*pi*n)*exp(n*log(n*exp(-1)))*(1 + 1/(12*n) + 1/(288*n*n) - 139/(51840*n*n*n) - 571/(2488320*n*n*n*n)) } - printf(OFMT ORS, p) + printf(OFMT ORS, f) } |