From 72fdb25210c579beaabc35cae7ec803436887f20 Mon Sep 17 00:00:00 2001
From: wukong <wukong@longaeva>
Date: Sat, 20 Oct 2018 11:47:10 -0700
Subject: merged lin_reg python files, added comment on Ramanujan's version of
 gamma function to sterling_approx;

---
 lin_reg.py          | 50 --------------------------------------------------
 lin_reg_example.py  |  6 +++---
 sterling_approx.awk | 18 ++++++++++++------
 3 files changed, 15 insertions(+), 59 deletions(-)
 delete mode 100644 lin_reg.py

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)
 }
-- 
cgit v1.2.3