Error estimate supports FLEETWIDE_MEANS

The error estimation tool now supports all 1.1 report types 🎉

Like *NumericStats reports, FleetwideMeans requires the user to pass in
an estimated lower bound on the denominator (number of contributing
devices).

Unit tests updated.

Change-Id: Ifb9e825f1a455c251290d04ab620f8bac733ff61
Reviewed-on: https://fuchsia-review.googlesource.com/c/cobalt/+/471660
Commit-Queue: Jared Weinstein <jaredweinstein@google.com>
Reviewed-by: Pasin Manurangsi <pasin@google.com>
Reviewed-by: Alexandre Zani <azani@google.com>
diff --git a/src/bin/config_parser/src/privacy/error_calculator.go b/src/bin/config_parser/src/privacy/error_calculator.go
index 34be217..af03066 100644
--- a/src/bin/config_parser/src/privacy/error_calculator.go
+++ b/src/bin/config_parser/src/privacy/error_calculator.go
@@ -37,8 +37,6 @@
 }
 
 // Given a |metric|, |report|, and |params|, estimates the report row error.
-//
-// TODO(jaredweinstein): implement estimates for other report types.
 func (e *ErrorCalculator) Estimate(metric *config.MetricDefinition, report *config.ReportDefinition, epsilon float64, population uint64, minDenominatorEstimate uint64) (estimate float64, err error) {
 	sparsity, err := getSparsityForReport(metric, report)
 	if err != nil {
@@ -51,9 +49,8 @@
 		return -1, err
 	}
 
-	reportType := report.GetReportType()
 	var errorEstimate float64
-	switch reportType {
+	switch report.GetReportType() {
 	case config.ReportDefinition_UNIQUE_DEVICE_HISTOGRAMS:
 		fallthrough // Calculate RMSE Error for each bucket.
 	case config.ReportDefinition_UNIQUE_DEVICE_COUNTS:
@@ -66,24 +63,23 @@
 	case config.ReportDefinition_FLEETWIDE_HISTOGRAMS:
 		errorEstimate = MultipleContributionRapporRMSE(population, privacyEncodingParams, 1)
 	case config.ReportDefinition_UNIQUE_DEVICE_NUMERIC_STATS:
-		if minDenominatorEstimate == 0 {
-			return -1, fmt.Errorf("User estimate for lower bound on unnoised denominator required for %s.", reportType)
-		}
-		if report.MaxValue == 0 && report.MinValue == 0 {
-			return -1, fmt.Errorf("MinValue and MaxValue required to estimate error for %s", reportType)
+		if err := meanReportConfigurationError(report, minDenominatorEstimate); err != nil {
+			return -1, err
 		}
 		errorEstimate = NumericStatsRapporRMSE(population, privacyEncodingParams, minDenominatorEstimate, report)
 	case config.ReportDefinition_HOURLY_VALUE_NUMERIC_STATS:
-		if minDenominatorEstimate == 0 {
-			return -1, fmt.Errorf("User estimate for lower bound on unnoised denominator required for %s.", reportType)
-		}
-		if report.MaxValue == 0 && report.MinValue == 0 {
-			return -1, fmt.Errorf("MinValue and MaxValue required to estimate error for %s", reportType)
+		if err := meanReportConfigurationError(report, minDenominatorEstimate); err != nil {
+			return -1, err
 		}
 		// (24 * population) pseudo-users to account for hourly values.
 		errorEstimate = NumericStatsRapporRMSE(24*population, privacyEncodingParams, minDenominatorEstimate, report)
+	case config.ReportDefinition_FLEETWIDE_MEANS:
+		if err := meanReportConfigurationError(report, minDenominatorEstimate); err != nil {
+			return -1, err
+		}
+		errorEstimate = FleetwideMeansRapporRMSE(24*population, privacyEncodingParams, minDenominatorEstimate, report)
 	default:
-		reportType := config.ReportDefinition_ReportType_name[int32(reportType)]
+		reportType := config.ReportDefinition_ReportType_name[int32(report.GetReportType())]
 		return -1, fmt.Errorf("Error estimation is not supported for reports of type %s.", reportType)
 	}
 
@@ -122,7 +118,25 @@
 	a := params.ProbBitFlip / (1 - 2*params.ProbBitFlip)
 	M := uint64(math.Max(math.Abs(float64(report.MinValue)), math.Abs(float64(report.MaxValue))))
 
-	// Numerator is aggregated like FleetwideOccurenceCounts.
+	// Numerator is aggregated like FleetwideOccurrenceCounts.
+	rmseNumerator := MultipleContributionRapporRMSE(population, params, M)
+	mseNumerator := math.Pow(rmseNumerator, 2)
+
+	// Denominator is aggregated like UniqueDeviceCounts.
+	rmseDenominator := SingleContributionRapporRMSE(population, params)
+	mseDenominator := math.Pow(rmseDenominator, 2)
+
+	return meanRapporRMSE(mseNumerator, mseDenominator, sigma, a, float64(M), float64(minDenominatorEstimate))
+}
+
+func FleetwideMeansRapporRMSE(population uint64, params PrivacyEncodingParams, minDenominatorEstimate uint64, report *config.ReportDefinition) float64 {
+	t := float64(report.MaxCount)
+	sigma := SingleContributionRapporRMSE(population, params) *
+		math.Sqrt(float64(t*(t+1)*(2*t+1)/6))
+	a := params.ProbBitFlip * t / (1 - 2*params.ProbBitFlip)
+	M := uint64(math.Max(math.Abs(float64(report.MinValue)), math.Abs(float64(report.MaxValue))))
+
+	// Numerator is aggregated like FleetwideOccurrenceCounts.
 	rmseNumerator := MultipleContributionRapporRMSE(population, params, M)
 	mseNumerator := math.Pow(rmseNumerator, 2)
 
@@ -165,3 +179,13 @@
 func h(u float64) float64 {
 	return (1+u)*math.Log(1+u) - u
 }
+
+func meanReportConfigurationError(report *config.ReportDefinition, minDenominatorEstimate uint64) error {
+	if minDenominatorEstimate == 0 {
+		return fmt.Errorf("user estimate for lower bound on unnoised denominator required for %s", report.GetReportType())
+	}
+	if report.MaxValue == 0 && report.MinValue == 0 {
+		return fmt.Errorf("MinValue and MaxValue required to estimate error for %s", report.GetReportType())
+	}
+	return nil
+}
diff --git a/src/bin/config_parser/src/privacy/error_calculator_test.go b/src/bin/config_parser/src/privacy/error_calculator_test.go
index eca0d0f..ebb4eab 100644
--- a/src/bin/config_parser/src/privacy/error_calculator_test.go
+++ b/src/bin/config_parser/src/privacy/error_calculator_test.go
@@ -75,9 +75,24 @@
 		ReportName: "UniqueDeviceNumericStatsMissingMaxValue",
 		ReportType: config.ReportDefinition_UNIQUE_DEVICE_NUMERIC_STATS,
 	}
-	unsupportedReport := config.ReportDefinition{
+	fleetwideMeans := config.ReportDefinition{
 		ReportName: "FleetwideMeans",
 		ReportType: config.ReportDefinition_FLEETWIDE_MEANS,
+		MaxValue:   1,
+		MaxCount:   1,
+	}
+	fleetwideMeansHighMaxCount := config.ReportDefinition{
+		ReportName: "FleetwideMeansHighMaxCount",
+		ReportType: config.ReportDefinition_FLEETWIDE_MEANS,
+		MaxValue:   1,
+		MaxCount:   4,
+	}
+	fleetwideMeansMissingMaxValue := config.ReportDefinition{
+		ReportName: "FleetwideMeansMissingMaxValue",
+		ReportType: config.ReportDefinition_FLEETWIDE_MEANS,
+	}
+	unsupportedReport := config.ReportDefinition{
+		ReportName: "UnsupportedReport",
 	}
 
 	type args struct {
@@ -121,12 +136,21 @@
 		{args{&testMetric, &hourlyValueNumericStats, 10, 10000, 500}, true, 0.09424308127635954},
 		{args{&testMetric, &hourlyValueNumericStats, 1, 20000, 500}, true, 0.20996453488325503},
 		{args{&testMetric, &hourlyValueNumericStats, 1, 10000, 2000}, true, 0.03425744055252169},
+		{args{&testMetric, &fleetwideMeans, 1, 10000, 500}, true, 0.14403990551821316},
+		{args{&testMetric, &fleetwideMeans, 10, 10000, 500}, true, 0.09424308127635954},
+		{args{&testMetric, &fleetwideMeans, 1, 20000, 500}, true, 0.20996453488325503},
+		{args{&testMetric, &fleetwideMeans, 1, 10000, 2000}, true, 0.03425744055252169},
+		{args{&testMetric, &fleetwideMeansHighMaxCount, 1, 10000, 500}, true, 14.67564563395404},
+		{args{&testMetric, &fleetwideMeansHighMaxCount, 10, 10000, 500}, true, 0.3612570115350953},
+		{args{&testMetric, &fleetwideMeansHighMaxCount, 1, 20000, 500}, true, 86.18119817858194},
+		{args{&testMetric, &fleetwideMeansHighMaxCount, 1, 10000, 2000}, true, 0.11189363618734707},
 
 		// Invalid input
 		// Mean Report missing minDenominatorEstimate
 		{args{&testMetric, &uniqueDeviceNumericStats, 1, 10000, 0}, false, 0},
 		{args{&testMetric, &hourlyValueNumericStats, 1, 10000, 0}, false, 0},
 		{args{&testMetric, &uniqueDeviceNumericStatsMissingMaxValue, 1, 10000, 500}, false, 0},
+		{args{&testMetric, &fleetwideMeansMissingMaxValue, 1, 10000, 500}, false, 0},
 
 		// This report type is not currently supported.
 		{args{&testMetric, &unsupportedReport, 1, 10000, 0}, false, 0},