added fib.awk, a solution to fib eqn (example in Hamming text);

removed pi.sh (redundant);
added corrections to quad_reg

3 files changed, 26 insertions, 35 deletions

diff --git a/fib.awk b/fib.awk
new file mode 100644
index 0000000..6107c4c
--- /dev/null
+++ b/fib.awk
@@ -0,0 +1,10 @@
+#!/usr/bin/awk -f
+
+### fib.awk
+# fib sol'n from Hamming
+
+BEGIN {
+    n = ARGV[1]
+    printf(OFMT ORS,
+        (1/sqrt(5))*((1.0 + sqrt(5))/2.0)^n  - (1/sqrt(5))*((1.0 - sqrt(5))/2.0)^n)
+}
diff --git a/pi.sh b/pi.sh
deleted file mode 100644
index 83af08a..0000000
--- a/pi.sh
+++ /dev/null
@@ -1,24 +0,0 @@
-#!/bin/sh
-
-### pi.sh, https://en.wikipedia.org/wiki/Pi
-# In 1706 John Machin used the Gregory–Leibniz series to produce an algorithm
-# that converged much faster. Machin reached 100 digits of π with this
-# formula. Other mathematicians created variants, now known as Machin-like
-# formulae, that were used to set several successive records for calculating
-# digits of π. Machin-like formulae remained the best-known method for
-# calculating π well into the age of computers, and were used to set records
-# for 250 years, culminating in a 620-digit approximation in 1946 by Daniel
-# Ferguson–the best approximation achieved without the aid of a calculating
-# device.
-
-awk -v fig=${1} '
-    function pi() {
-        return 4*(4*atan2(1,5) - atan2(1,239))
-    }
-    BEGIN {
-        if (fig <= 0)
-            fig = 6
-        str = "%." fig "g\n"
-        printf(str, pi())
-    }
-'
diff --git a/quad_reg.awk b/quad_reg.awk
index 1b30afd..0620808 100644
--- a/quad_reg.awk
+++ b/quad_reg.awk
@@ -1,7 +1,7 @@
 #!/usr/bin/awk -f
 
 ### quad_reg.awk
-# simple linear regression between columns
+# quadratic regression along columns
 
 BEGIN {
     OFMT="%.9g"
@@ -36,8 +36,10 @@ NF > 0 {
 
             ### difference from the mean
             delta[y] = $y - mean[y]
+            delta2[y] = $y*$y - mean2[y]
             sum_delta[y] += delta[y]
-            sum_delta2[y] += delta[y]*delta[y]
+            sum2_delta[y] += delta[y]*delta[y]
+            sum_delta2[y] += delta2[y]
 
             ### sample variance
             (count[y] > 1) ? var[y] = sum_delta2[y]/(count[y] - 1) : var[y] = ""
@@ -48,14 +50,17 @@ NF > 0 {
                 sum_xy[x,y] += $x*$y
                 sum_x2y[x,y] += $x*$x*$y
                 sum_delta_xy[x,y] += delta[x]*delta[y]
+                sum_delta_xx2[x,y] += delta[x]*delta2[x]
+                sum_delta_x2y[x,y] += delta2[x]*delta[y]
+                sum_delta_x2x2[x,y] += delta2[x]*delta2[x]
 
                 # covariances
                 if (count[x,y] > 1) {
-                    s_xx[x,y] = sum2[x]/(count[x,y]) - mean[x]*mean[x]
-                    s_xy[x,y] = sum_xy[x,y]/(count[x,y]) - mean[x]*mean[y]
-                    s_xx2[x,y] = sum3[x]/(count[x,y]) - mean[x]*mean2[x]
-                    s_x2x2[x,y] = sum4[x]/(count[x,y]) - mean2[x]*mean2[x]
-                    s_x2y[x,y] = sum_x2y[x]/(count[x,y]) - mean2[x]*mean[y]
+                    s_xx[x,y] = sum2_delta[x]/(count[x,y] - 1)
+                    s_xy[x,y] = sum_delta_xy[x,y]/(count[x,y] - 1)
+                    s_xx2[x,y] = sum_delta_xx2[x,y]/(count[x,y] - 1)
+                    s_x2x2[x,y] = sum_delta_x2x2[x,y]/(count[x,y] - 1)
+                    s_x2y[x,y] = sum_delta_x2y[x,y]/(count[x,y] - 1)
                 }
 
                 bc_den[x,y] = (s_xx[x,y]*s_x2x2[x,y] - s_xx2[x,y]*s_xx2[x,y])
@@ -67,14 +72,14 @@ NF > 0 {
                     c[x,y] = 0
                     b[x,y] = 0
                 }
-                a[x,y] = mean[y] - b[x,y]*mean[x] - c[x,y]*mean[x]*mean[x]
+                a[x,y] = mean[y] - b[x,y]*mean[x] - c[x,y]*mean2[x]
 
                 # error estimate
                 err[x,y] = ($y - (a[x,y] + b[x,y]*$x + c[x,y]*$x*$x))
                 sum_err2[x,y] += err[x,y]*err[x,y]
 
                 # correlation
-                sum_delta2[y] ? r2[x,y] = sum_err2[x,y]/sum_delta2[y] : r2[x,y] = 1
+                sum_delta2[y] ? r[x,y] = sqrt(1 - sum_err2[x,y]/sum_delta2[y]) : r[x,y] = 0
             }
         }
         else
@@ -85,9 +90,9 @@ NF > 0 {
 END {
     for (x=1; x<=nf_max; x++) {
         for (y=1; y<=nf_max; y++) {
-            if (x != y && r2[x,y]) {
+            if (x != y && r[x,y]) {
                 printf(OFMT OFS "(%s)" OFS " = (" OFMT ")(%s)^2" OFS " + (" OFMT ")(%s)" OFS " + (" OFMT ")" ORS,
-                    10.0*log(r2[x,y])/log(10), header[y], c[x,y], header[x], b[x,y], header[x], a[x,y])
+                    10.0*log(r[x,y]*r[x,y])/log(10), header[y], c[x,y], header[x], b[x,y], header[x], a[x,y])
             }
         }
     }