tiny.simd.Stats

1pub const Stats = struct {2    count_value: i64 = 0,3    minimum: f32 = std.math.floatMax(f32),4    maximum: f32 = -std.math.floatMax(f32),5    sum_log: f64 = 0,6    moment_1: f64 = 0,7    moment_2: f64 = 0,8    moment_3: f64 = 0,9    moment_4: f64 = 0,1011    pub const max_count: i64 = 3_037_000_499;12    pub const no_count: u8 = 1;13    pub const no_mean_sd: u8 = 2;14    pub const no_min_max: u8 = 4;15    pub const no_skew_kurt: u8 = 8;16    pub const no_geometric_mean: u8 = 16;17    pub const all_exclusions: u8 = no_count |18        no_mean_sd |19        no_min_max |20        no_skew_kurt |21        no_geometric_mean;2223    pub fn notify(self: *@This(), value: f32) void {24        std.debug.assert(self.count_value >= 0);25        std.debug.assert(self.count_value < max_count);26        self.count_value += 1;27        self.minimum = @min(self.minimum, value);28        self.maximum = @max(self.maximum, value);29        self.sum_log += @log(@as(f64, value));3031        const count_f64: f64 = @floatFromInt(self.count_value);32        const delta = @as(f64, value) - self.moment_1;33        const delta_div_count = delta / count_f64;34        const delta2_count_minus_1_div_count = delta * (count_f64 - 1) * delta_div_count;35        const count_polynomial = count_f64 * count_f64 - 3 * count_f64 + 3;36        self.moment_1 += delta_div_count;37        self.moment_4 += delta_div_count *38            (delta_div_count *39                (delta2_count_minus_1_div_count * count_polynomial + 6 * self.moment_2) -40                4 * self.moment_3);41        self.moment_3 += delta_div_count *42            (delta2_count_minus_1_div_count * (count_f64 - 2) - 3 * self.moment_2);43        self.moment_2 += delta2_count_minus_1_div_count;44    }4546    pub fn assimilate(self: *@This(), other: *const @This()) void {47        std.debug.assert(self.count_value >= 0);48        std.debug.assert(other.count_value >= 0);49        std.debug.assert(self.count_value <= max_count - other.count_value);50        const total_count = self.count_value + other.count_value;51        if (total_count == 0) return;5253        self.minimum = @min(self.minimum, other.minimum);54        self.maximum = @max(self.maximum, other.maximum);55        self.sum_log += other.sum_log;5657        const own_count: f64 = @floatFromInt(self.count_value);58        const other_count: f64 = @floatFromInt(other.count_value);59        const total_count_f64: f64 = @floatFromInt(total_count);60        const product = own_count * other_count;61        const count_squared = own_count * own_count;62        const other_count_squared = other_count * other_count;63        const total_count_squared = total_count_f64 * total_count_f64;64        const total_count_cubed = total_count_squared * total_count_f64;65        const inverse_total_count = 1 / total_count_f64;66        const inverse_total_count_squared = 1 / total_count_squared;6768        const delta = other.moment_1 - self.moment_1;69        const delta_squared = delta * delta;70        const delta_cubed = delta * delta_squared;71        const delta_fourth = delta_squared * delta_squared;7273        self.moment_1 = (own_count * self.moment_1 + other_count * other.moment_1) * inverse_total_count;74        const new_moment_2 = self.moment_2 + other.moment_2 +75            delta_squared * product * inverse_total_count;76        const new_moment_3 = self.moment_3 + other.moment_3 +77            delta_cubed * product * (own_count - other_count) * inverse_total_count_squared +78            3 * delta * (own_count * other.moment_2 - other_count * self.moment_2) * inverse_total_count;79        self.moment_4 += other.moment_4 +80            delta_fourth * product * (count_squared - product + other_count_squared) /81                total_count_cubed +82            6 * delta_squared *83                (count_squared * other.moment_2 + other_count_squared * self.moment_2) *84                inverse_total_count_squared +85            4 * delta * (own_count * other.moment_3 - other_count * self.moment_3) *86                inverse_total_count;87        self.moment_2 = new_moment_2;88        self.moment_3 = new_moment_3;89        self.count_value = total_count;90    }9192    pub fn count(self: *const @This()) i64 {93        return self.count_value;94    }9596    pub fn min(self: *const @This()) f32 {97        return self.minimum;98    }99100    pub fn max(self: *const @This()) f32 {101        return self.maximum;102    }103104    pub fn geometricMean(self: *const @This()) f64 {105        if (self.count_value == 0) return 0;106        return @exp(self.sum_log / @as(f64, @floatFromInt(self.count_value)));107    }108109    pub fn mean(self: *const @This()) f64 {110        return self.moment_1;111    }112113    pub fn sampleVariance(self: *const @This()) f64 {114        if (self.count_value == 0) return 0;115        return self.moment_2 / @as(f64, @floatFromInt(self.count_value));116    }117118    pub fn variance(self: *const @This()) f64 {119        if (self.count_value == 0) return 0;120        if (self.count_value == 1) return self.moment_2;121        return self.moment_2 / @as(f64, @floatFromInt(self.count_value - 1));122    }123124    pub fn standardDeviation(self: *const @This()) f64 {125        return @sqrt(self.variance());126    }127128    pub fn sampleSkewness(self: *const @This()) f64 {129        if (@abs(self.moment_2) < 1e-7) return 0;130        const count_f64: f64 = @floatFromInt(self.count_value);131        return self.moment_3 * @sqrt(count_f64) /132            std.math.pow(f64, self.moment_2, 1.5);133    }134135    pub fn skewness(self: *const @This()) f64 {136        if (self.count_value == 0) return 0;137        const count_f64: f64 = @floatFromInt(self.count_value);138        const ratio = (count_f64 - 1) / count_f64;139        return self.sampleSkewness() * std.math.pow(f64, ratio, 1.5);140    }141142    pub fn sampleKurtosis(self: *const @This()) f64 {143        if (@abs(self.moment_2) < 1e-7) return 0;144        const count_f64: f64 = @floatFromInt(self.count_value);145        return self.moment_4 * count_f64 / (self.moment_2 * self.moment_2);146    }147148    pub fn kurtosis(self: *const @This()) f64 {149        if (self.count_value == 0) return 0;150        const count_f64: f64 = @floatFromInt(self.count_value);151        const ratio = (count_f64 - 1) / count_f64;152        return self.sampleKurtosis() * ratio * ratio;153    }154155    pub fn mu1(self: *const @This()) f64 {156        std.debug.assert(self.count_value != 0);157        return self.moment_1;158    }159160    pub fn mu2(self: *const @This()) f64 {161        return self.centralMoment(self.moment_2);162    }163164    pub fn mu3(self: *const @This()) f64 {165        return self.centralMoment(self.moment_3);166    }167168    pub fn mu4(self: *const @This()) f64 {169        return self.centralMoment(self.moment_4);170    }171172    pub fn write(173        self: *const @This(),174        writer: *std.Io.Writer,175        exclude: u8,176    ) std.Io.Writer.Error!void {177        if (self.count_value == 0) return writer.writeAll("(none)");178        if (exclude & no_count == 0) {179            try writer.print("Count={d:9} ", .{@as(u64, @intCast(self.count_value))});180        }181        if (exclude & no_mean_sd == 0) {182            try writer.writeAll("Mean=");183            try writeCScientific(writer, self.mean(), 10, 3);184            try writer.writeAll(" SD=");185            try writeCScientific(writer, self.standardDeviation(), 8, 2);186            try writer.writeByte(' ');187        }188        if (exclude & no_min_max == 0) {189            try writer.writeAll("Min=");190            try writeCScientific(writer, self.minimum, 10, 3);191            try writer.writeAll(" Max=");192            try writeCScientific(writer, self.maximum, 10, 3);193            try writer.writeByte(' ');194        }195        if (exclude & no_skew_kurt == 0) {196            try writer.print("Skew={d:5.2} Kurt={d:7.2} ", .{197                self.skewness(),198                self.kurtosis(),199            });200        }201        if (exclude & no_geometric_mean == 0) {202            try writer.print("GeoMean={d:9.6} ", .{self.geometricMean()});203        }204    }205206    pub fn reset(self: *@This()) void {207        self.* = .{};208    }209210    fn centralMoment(self: *const @This(), moment: f64) f64 {211        std.debug.assert(self.count_value != 0);212        return moment / @as(f64, @floatFromInt(self.count_value));213    }214};

1pub const Stats = stats.Stats;