Fixed a problem where stddev in C++ and Go was calculated from variance instead of sample variance. R uses sample variance so defaulting to that.

cpp/include/stats.hpp

@@ -30,9 +30,13 @@ struct Stats {
     return sum / n;
   }
 
-  inline double stddev() {
+  inline double stddev(bool use_sample=true) {
+    if(use_sample) {
+      return std::sqrt(sample_variance());
+    } else {
       return std::sqrt(variance());
     }
+  }
 
   inline double variance() {
     return (sumsq - (sum * sum / n)) / n;

cpp/include/welford.hpp

@@ -24,9 +24,13 @@ struct WelfordStats {
     return Mean;
   }
 
-  inline double stddev() {
+  inline double stddev(bool use_sample=true) {
+    if(use_sample) {
+      return std::sqrt(sample_variance());
+    } else {
       return std::sqrt(variance());
     }
+  }
 
   inline double variance() {
     return M2 / n;

cpp/tests/confirm.r

@@ -1,3 +1,9 @@
 samples <- read.table("samples.txt")
-print(samples)
+#print(samples)
+print("Mean V1")
+mean(samples$V1)
+print("Stddev V1")
+sd(samples$V1)
+print("Variance V1")
+var(samples$V1)
 t.test(samples$V1, samples$V2)

cpp/tests/stats_tests.cpp

@@ -68,7 +68,7 @@ namespace stats_tests {
     ALMOST_EQUAL(var, stats.variance(), 5);
 
     double stddev = std::sqrt(var);
-    ALMOST_EQUAL(stddev, stats.stddev(), 5);
+    ALMOST_EQUAL(stddev, stats.stddev(false), 5);
   }
 
   void test_welford() {
@@ -123,10 +123,10 @@ namespace stats_tests {
     ALMOST_EQUAL(stats_test.dof, welf_test.dof, 1);
     ALMOST_EQUAL(stats_test.p_val, welf_test.p_val, 5);
 
-    fmt::println("welford samples t-test: t-val: {}, dof: {}, p-val: {}",
+    fmt::println("mean: {}, sd: {} var: {}; welford samples t-test: t-val: {}, dof: {}, p-val: {}",
-        welf_test.t_stat, welf_test.dof, welf_test.p_val);
+        welf.mean(), welf.stddev(), welf.sample_variance(), welf_test.t_stat, welf_test.dof, welf_test.p_val);
-    fmt::println("naive algebra samples t-test: t-val: {}, dof: {}, p-val: {}",
+    fmt::println("mean: {}, sd: {}; var: {}; naive algebra samples t-test: t-val: {}, dof: {}, p-val: {}",
-        stats_test.t_stat, stats_test.dof, stats_test.p_val);
+        stats.mean(), stats.stddev(), stats.sample_variance(), stats_test.t_stat, stats_test.dof, stats_test.p_val);
   }
 
   void failing_bad_t_test() {

go/magic_stats/Makefile

@@ -1,4 +1,10 @@
-
+R_SCRIPT=C:/Users/lcthw/AppData/Local/Programs/R/R-4.6.0/bin/Rscript.exe
 
 all:
 	go build .
+
+test:
+	./magic_stats
+
+r_test: test
+	$(R_SCRIPT) ../../cpp/tests/confirm.r

go/magic_stats/go.mod

@@ -1,5 +1,3 @@
 module lcthw.dev/learn-code-the-hard-way/magic_stats/go/magic_stats
 
 go 1.25.3
-
-require gonum.org/v1/gonum v0.17.0

go/magic_stats/go.sum

@@ -1,2 +0,0 @@
-gonum.org/v1/gonum v0.17.0 h1:VbpOemQlsSMrYmn7T2OUvQ4dqxQXU+ouZFQsZOx50z4=
-gonum.org/v1/gonum v0.17.0/go.mod h1:El3tOrEuMpv2UdMrbNlKEh9vd86bmQ6vqIcDwxEOc1E=

go/magic_stats/main.go

@@ -3,6 +3,9 @@ package main
 import (
   "fmt"
   "math"
+  "math/rand/v2"
+  "os"
+  "log"
 )
 
 type Welford struct {
@@ -15,7 +18,7 @@ type Welford struct {
 }
 
 type TTest struct {
-  TState float64
+  TStat float64
   DOF float64
   PVal float64
 }
@@ -44,10 +47,26 @@ func (stats *Welford) SampleVariance() float64 {
 }
 
 func (stats *Welford) StdDev() float64 {
-  return math.Sqrt(stats.Variance())
+  return math.Sqrt(stats.SampleVariance())
 }
 
-func (stats *Welford) TTest(second *Welford) {
+func Simpson(a, b float64, n int, f func(float64)float64) float64 {
+  n_float := float64(n)
+  h := ( b - a ) / n_float
+  sum := 0.0
+  for i := 0 ; i < n; i++ {
+    x := a + float64(i) * h
+    sum += ( f(x) + 4.0 * f(x + h / 2.0) + f(x + h) ) / 6.0
+  }
+
+  return sum * h
+}
+
+func (stats *Welford) Dump() {
+  fmt.Println("Mean:", stats.Mean, "N:", stats.N, "M2:", stats.M2, "Min:", stats.Min, "Max:", stats.Max)
+}
+
+func (stats *Welford) TTest(second *Welford) TTest {
   n1 := stats.N
   n2 := second.N
   mean1 := stats.Mean
@@ -74,29 +93,38 @@ func (stats *Welford) TTest(second *Welford) {
     return math.Pow(r, dof / 2.0 - 1.0) / math.Sqrt(1.0 - r)
   }
 
-  p_val := Simpson(0.0, b, 10000, fn) / gamma
+  p_val := Simpson(0.0, b, 10000, fn) / gamm
 
-  fmt.Println("t_stats:", t_stat, "dof", dof, "gamm:", gamm, "b:", b)
+  return TTest{t_stat, dof, p_val}
 }
 
 func main() {
   var welf Welford
   var welf2 Welford
+  skew_factor := rand.Float64() * 10
 
-  welf.Sample(1.2)
+  samples_file, err := os.Create("samples.txt")
-  welf.Sample(2.2)
+  if err != nil { log.Fatal(err) }
-  welf.Sample(3.2)
+  defer samples_file.Close()
-  welf.Sample(4.2)
 
-  welf2.Sample(1.1)
+  for i := 0; i < 100; i++ {
-  welf2.Sample(2.1)
+    sample := rand.Float64() * 100
-  welf2.Sample(3.1)
+    welf.Sample(sample)
-  welf2.Sample(4.1)
+    welf2.Sample(sample + skew_factor)
+
+    sample_line := fmt.Sprintf("%f %f\n", sample, sample + skew_factor)
+    samples_file.Write([]byte(sample_line))
+  }
 
-  welf.TTest(&welf)
+  welf.Dump()
+  welf2.Dump()
+  result := welf.TTest(&welf2)
 
   fmt.Println(
-    "Mean:", welf.Mean,
+    "----\nTTest Results:\nMean:", welf.Mean,
-    "Variance: ", welf.Variance(),
+    "\nVariance: ", welf.Variance(),
-    "Stddev: ", welf.StdDev());
+    "\nStddev: ", welf.StdDev(),
+    "\nTStat: ", result.TStat,
+    "\nDOF: ", result.DOF,
+    "\nPVal: ", result.PVal);
 }