add WeightedStdDevDirectionDeg

libwx.go

@@ -294,8 +294,17 @@ func WeightedAvgDirectionDeg(degrees []Degree, weights []float64) (Degree, error
 	return radToDeg(circularMean(degToRadSlice(clampedDegSlice(degrees)), weights)), nil
 }
 
-// StdDevDeg calculates the circular standard deviation of the given set of angles (in degrees).
+// StdDevDirectionDeg calculates the circular standard deviation of the given set of angles (in degrees).
 // This is useful to find e.g. the variability of wind direction.
-func StdDevDeg(degrees []Degree) Degree {
-	return radToDeg(circularStdDev(degToRadSlice(clampedDegSlice(degrees))))
+func StdDevDirectionDeg(degrees []Degree) Degree {
+	return radToDeg(circularStdDev(degToRadSlice(clampedDegSlice(degrees)), nil))
+}
+
+// WeightedStdDevDirectionDeg calculates the circular standard deviation of the given set of angles (in degrees).
+// This is useful to find e.g. the variability of wind direction, weighted by wind speed.
+func WeightedStdDevDirectionDeg(degrees []Degree, weights []float64) (Degree, error) {
+	if len(degrees) != len(weights) {
+		return 0.0, ErrMismatchedInputLength
+	}
+	return radToDeg(circularStdDev(degToRadSlice(clampedDegSlice(degrees)), weights)), nil
 }

stat.go

@@ -2,7 +2,6 @@ package libwx
 
 import "math"
 
-// borrowed from gonum:
 func circularMean(x, weights []float64) float64 {
 	if weights != nil && len(x) != len(weights) {
 		panic("stat: slice length mismatch")
@@ -24,12 +23,26 @@ func circularMean(x, weights []float64) float64 {
 	return math.Atan2(aY, aX)
 }
 
-// from my (open) gonum PR which doesn't (yet) support weights:
-func circularStdDev(x []float64) float64 {
+func circularStdDev(x []float64, weights []float64) float64 {
+	if weights != nil && len(x) != len(weights) {
+		panic("stat: slice length mismatch")
+	}
+
 	var aX, aY float64
-	for _, v := range x {
-		aX += math.Cos(v)
-		aY += math.Sin(v)
+	if weights != nil {
+		var sumW float64
+		for i, v := range x {
+			w := weights[i]
+			sumW += w
+			aX += w * math.Cos(v)
+			aY += w * math.Sin(v)
+		}
+		return math.Sqrt(-2 * math.Log(math.Hypot(aY, aX)/sumW))
+	} else {
+		for _, v := range x {
+			aX += math.Cos(v)
+			aY += math.Sin(v)
+		}
+		return math.Sqrt(-2 * math.Log(math.Hypot(aY, aX)/float64(len(x))))
 	}
-	return math.Sqrt(-2 * math.Log(math.Sqrt(aY*aY+aX*aX)/float64(len(x))))
 }