diff --git a/go/magic_stats/Makefile b/go/magic_stats/Makefile
new file mode 100644
index 0000000..907b0d2
--- /dev/null
+++ b/go/magic_stats/Makefile
@@ -0,0 +1,4 @@
+
+
+all:
+	go build .
diff --git a/go/magic_stats/go.mod b/go/magic_stats/go.mod
new file mode 100644
index 0000000..ddcb982
--- /dev/null
+++ b/go/magic_stats/go.mod
@@ -0,0 +1,5 @@
+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
diff --git a/go/magic_stats/go.sum b/go/magic_stats/go.sum
new file mode 100644
index 0000000..d45948d
--- /dev/null
+++ b/go/magic_stats/go.sum
@@ -0,0 +1,2 @@
+gonum.org/v1/gonum v0.17.0 h1:VbpOemQlsSMrYmn7T2OUvQ4dqxQXU+ouZFQsZOx50z4=
+gonum.org/v1/gonum v0.17.0/go.mod h1:El3tOrEuMpv2UdMrbNlKEh9vd86bmQ6vqIcDwxEOc1E=
diff --git a/go/magic_stats/main.go b/go/magic_stats/main.go
new file mode 100644
index 0000000..6bc5924
--- /dev/null
+++ b/go/magic_stats/main.go
@@ -0,0 +1,102 @@
+package main
+
+import (
+  "fmt"
+  "math"
+)
+
+type Welford struct {
+  Mean float64
+  N float64
+  M2 float64
+
+  Min float64
+  Max float64
+}
+
+type TTest struct {
+  TState float64
+  DOF float64
+  PVal float64
+}
+
+func (stats *Welford) Sample(s float64) {
+  stats.N += 1
+  old_mean := stats.Mean
+  stats.Mean += (s - stats.Mean) / stats.N
+  stats.M2 += (s - old_mean) * (s - stats.Mean)
+
+  if stats.N == 0 {
+    stats.Min = s
+    stats.Max = s
+  } else {
+    if stats.Min > s { stats.Min = s }
+    if stats.Max < s { stats.Max = s }
+  }
+}
+
+func (stats *Welford) Variance() float64 {
+  return stats.M2 / stats.N
+}
+
+func (stats *Welford) SampleVariance() float64 {
+  return stats.M2 / (stats.N - 1)
+}
+
+func (stats *Welford) StdDev() float64 {
+  return math.Sqrt(stats.Variance())
+}
+
+func (stats *Welford) TTest(second *Welford) {
+  n1 := stats.N
+  n2 := second.N
+  mean1 := stats.Mean
+  mean2 := second.Mean
+  var1 := stats.SampleVariance()
+  var2 := second.SampleVariance()
+
+  delta_mean := mean1 - mean2
+  pooled_se := math.Sqrt((var1 / n1) + (var2 / n2))
+  t_stat := delta_mean / pooled_se
+
+  num := math.Pow((var1 / n1) + (var2 / n2), 2)
+  den := (math.Pow(var1 / n2, 2) / (n1 - 1)) + (math.Pow(var2 / n2, 2) / (n2 - 1))
+  dof := num / den
+
+  lga_1, _ := math.Lgamma(dof / 2.0)
+  lga_2, _ := math.Lgamma(0.5)
+  lga_3, _ := math.Lgamma(dof / 2.0 + 0.5)
+  gamm := math.Exp(lga_1 + lga_2 - lga_3)
+
+  b := dof / (t_stat * t_stat + dof)
+
+  fn := func (r float64) float64 {
+    return math.Pow(r, dof / 2.0 - 1.0) / math.Sqrt(1.0 - r)
+  }
+
+  p_val := Simpson(0.0, b, 10000, fn) / gamma
+
+  fmt.Println("t_stats:", t_stat, "dof", dof, "gamm:", gamm, "b:", b)
+}
+
+func main() {
+  var welf Welford
+  var welf2 Welford
+
+  welf.Sample(1.2)
+  welf.Sample(2.2)
+  welf.Sample(3.2)
+  welf.Sample(4.2)
+
+  welf2.Sample(1.1)
+  welf2.Sample(2.1)
+  welf2.Sample(3.1)
+  welf2.Sample(4.1)
+
+  welf.TTest(&welf)
+
+  fmt.Println(
+    "Mean:", welf.Mean,
+    "Variance: ", welf.Variance(),
+    "Stddev: ", welf.StdDev());
+}