mirror of
https://github.com/mihakralj/QuanTAlib.git
synced 2026-08-18 10:38:05 +00:00
Add Tukey's Biweight and WMAPE implementations with comprehensive tests and documentation
- Introduced Tukey's Biweight as a robust loss function, including mathematical foundation, usage patterns, and performance profile. - Added WMAPE (Weighted Mean Absolute Percentage Error) implementation, emphasizing its advantages for intermittent demand forecasting. - Created unit tests for WMAPE covering various scenarios including edge cases and batch calculations. - Documented both Tukey's Biweight and WMAPE with detailed explanations, properties, and common use cases.
This commit is contained in:
Miha Kralj
@@ -0,0 +1,476 @@
|
||||
using Xunit;
|
||||
|
||||
namespace QuanTAlib.Tests;
|
||||
|
||||
public class PseudoHuberTests
|
||||
{
|
||||
private const double Epsilon = 1e-10;
|
||||
private const int DefaultPeriod = 14;
|
||||
|
||||
#region Constructor Tests
|
||||
|
||||
[Fact]
|
||||
public void Constructor_ValidatesInput()
|
||||
{
|
||||
Assert.Throws<ArgumentException>(() => new PseudoHuber(0));
|
||||
Assert.Throws<ArgumentException>(() => new PseudoHuber(-1));
|
||||
Assert.Throws<ArgumentException>(() => new PseudoHuber(10, 0));
|
||||
Assert.Throws<ArgumentException>(() => new PseudoHuber(10, -1));
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Constructor_ValidPeriod_Succeeds()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(10);
|
||||
Assert.NotNull(pseudoHuber);
|
||||
Assert.Equal(10, pseudoHuber.WarmupPeriod);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Constructor_ValidDelta_Succeeds()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(10, 0.5);
|
||||
Assert.NotNull(pseudoHuber);
|
||||
Assert.Equal(0.5, pseudoHuber.Delta);
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
#region Property Tests
|
||||
|
||||
[Fact]
|
||||
public void Properties_Accessible()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(DefaultPeriod, 1.5);
|
||||
|
||||
Assert.Equal(0, pseudoHuber.Last.Value);
|
||||
Assert.False(pseudoHuber.IsHot);
|
||||
Assert.Contains("PseudoHuber", pseudoHuber.Name, StringComparison.Ordinal);
|
||||
Assert.Equal(1.5, pseudoHuber.Delta);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void IsHot_BecomesTrueWhenBufferFull()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(5);
|
||||
|
||||
Assert.False(pseudoHuber.IsHot);
|
||||
|
||||
for (int i = 1; i <= 4; i++)
|
||||
{
|
||||
pseudoHuber.Update(100 + i, 100.0);
|
||||
Assert.False(pseudoHuber.IsHot);
|
||||
}
|
||||
|
||||
pseudoHuber.Update(105, 100.0);
|
||||
Assert.True(pseudoHuber.IsHot);
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
#region Calculation Tests
|
||||
|
||||
[Fact]
|
||||
public void Calculate_PerfectPredictions_ReturnsZero()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(5);
|
||||
|
||||
for (int i = 0; i < 10; i++)
|
||||
{
|
||||
double value = 100 + i;
|
||||
pseudoHuber.Update(value, value);
|
||||
}
|
||||
|
||||
Assert.Equal(0.0, pseudoHuber.Last.Value, Epsilon);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_SmallErrors_ApproximatesL2()
|
||||
{
|
||||
// For small errors, Pseudo-Huber ≈ 0.5 * error²
|
||||
var pseudoHuber = new PseudoHuber(1, delta: 10.0);
|
||||
double error = 0.1; // Small relative to delta
|
||||
pseudoHuber.Update(100.0 + error, 100.0);
|
||||
|
||||
// Pseudo-Huber = δ² * (√(1 + (x/δ)²) - 1)
|
||||
// For small x/δ: √(1 + ε) ≈ 1 + ε/2, so loss ≈ δ² * (x/δ)²/2 = x²/2
|
||||
double expectedApprox = error * error / 2.0;
|
||||
double ratio = pseudoHuber.Last.Value / expectedApprox;
|
||||
|
||||
// Should be close to 1.0 for small errors
|
||||
Assert.InRange(ratio, 0.99, 1.01);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_LargeErrors_ApproximatesL1()
|
||||
{
|
||||
// For large errors, Pseudo-Huber ≈ δ * |error| - δ²/2
|
||||
var pseudoHuber = new PseudoHuber(1, delta: 1.0);
|
||||
double error = 100.0; // Large relative to delta
|
||||
pseudoHuber.Update(100.0 + error, 100.0);
|
||||
|
||||
// For large x: √(1 + (x/δ)²) ≈ |x/δ|
|
||||
// So loss ≈ δ² * (|x/δ| - 1) = δ|x| - δ²
|
||||
double expectedApprox = Math.Abs(error) - 1.0;
|
||||
double ratio = pseudoHuber.Last.Value / expectedApprox;
|
||||
|
||||
// Should be close to 1.0 for large errors
|
||||
Assert.InRange(ratio, 0.99, 1.01);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_SmoothTransition()
|
||||
{
|
||||
// Pseudo-Huber should be smooth across all error magnitudes
|
||||
var pseudoHuber = new PseudoHuber(1, delta: 1.0);
|
||||
double[] errors = { 0.01, 0.1, 0.5, 1.0, 2.0, 5.0, 10.0 };
|
||||
double[] losses = new double[errors.Length];
|
||||
|
||||
for (int i = 0; i < errors.Length; i++)
|
||||
{
|
||||
pseudoHuber.Reset();
|
||||
pseudoHuber.Update(100.0 + errors[i], 100.0);
|
||||
losses[i] = pseudoHuber.Last.Value;
|
||||
}
|
||||
|
||||
// Losses should be monotonically increasing
|
||||
for (int i = 1; i < losses.Length; i++)
|
||||
{
|
||||
Assert.True(losses[i] > losses[i - 1],
|
||||
$"Loss should increase: {losses[i - 1]} -> {losses[i]}");
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_Symmetry()
|
||||
{
|
||||
// Pseudo-Huber should be symmetric: loss(e) = loss(-e)
|
||||
var pseudoHuber1 = new PseudoHuber(1);
|
||||
var pseudoHuber2 = new PseudoHuber(1);
|
||||
|
||||
double error = 5.0;
|
||||
pseudoHuber1.Update(100.0 + error, 100.0); // Positive error
|
||||
pseudoHuber2.Update(100.0 - error, 100.0); // Negative error
|
||||
|
||||
Assert.Equal(pseudoHuber1.Last.Value, pseudoHuber2.Last.Value, Epsilon);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_DeltaEffectOnTransition()
|
||||
{
|
||||
// Larger delta means smoother transition, smaller delta means sharper
|
||||
var smallDelta = new PseudoHuber(1, delta: 0.5);
|
||||
var largeDelta = new PseudoHuber(1, delta: 2.0);
|
||||
|
||||
double error = 1.0; // Fixed error
|
||||
smallDelta.Update(100.0 + error, 100.0);
|
||||
largeDelta.Update(100.0 + error, 100.0);
|
||||
|
||||
// With large delta, the loss is more quadratic (smaller)
|
||||
// With small delta, the loss is more linear (larger relative to quadratic)
|
||||
// The raw loss values depend on the formula
|
||||
Assert.True(double.IsFinite(smallDelta.Last.Value));
|
||||
Assert.True(double.IsFinite(largeDelta.Last.Value));
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_ComparedToHuber()
|
||||
{
|
||||
// Pseudo-Huber should produce similar (but not identical) results to Huber
|
||||
var huber = new Huber(1, delta: 1.0);
|
||||
var pseudoHuber = new PseudoHuber(1, delta: 1.0);
|
||||
|
||||
// Test at various error magnitudes
|
||||
double[] errors = { 0.5, 1.0, 2.0 };
|
||||
|
||||
foreach (var error in errors)
|
||||
{
|
||||
huber.Reset();
|
||||
pseudoHuber.Reset();
|
||||
|
||||
huber.Update(100.0 + error, 100.0);
|
||||
pseudoHuber.Update(100.0 + error, 100.0);
|
||||
|
||||
// They should be in the same ballpark
|
||||
double ratio = pseudoHuber.Last.Value / huber.Last.Value;
|
||||
Assert.InRange(ratio, 0.5, 2.0); // Within factor of 2
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_AlwaysNonNegative()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(DefaultPeriod);
|
||||
var gbm = new GBM();
|
||||
|
||||
for (int i = 0; i < 100; i++)
|
||||
{
|
||||
var bar = gbm.Next();
|
||||
pseudoHuber.Update(bar.Close, bar.Close + (i % 2 == 0 ? 1 : -1) * (i + 1));
|
||||
Assert.True(pseudoHuber.Last.Value >= 0, "Pseudo-Huber loss should always be non-negative");
|
||||
}
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
#region State Management Tests
|
||||
|
||||
[Fact]
|
||||
public void Calculate_IsNew_False_UpdatesValue()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(5);
|
||||
|
||||
pseudoHuber.Update(100.0, 99.0);
|
||||
pseudoHuber.Update(101.0, 99.0, isNew: true);
|
||||
double beforeUpdate = pseudoHuber.Last.Value;
|
||||
|
||||
pseudoHuber.Update(105.0, 99.0, isNew: false);
|
||||
double afterUpdate = pseudoHuber.Last.Value;
|
||||
|
||||
Assert.NotEqual(beforeUpdate, afterUpdate);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void IterativeCorrections_RestoreToOriginalState()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(5);
|
||||
var gbm = new GBM();
|
||||
|
||||
// Feed 10 new values
|
||||
double tenthActual = 0, tenthPredicted = 0;
|
||||
for (int i = 0; i < 10; i++)
|
||||
{
|
||||
var bar = gbm.Next(isNew: true);
|
||||
tenthActual = bar.Close;
|
||||
tenthPredicted = bar.Close * 0.99;
|
||||
pseudoHuber.Update(tenthActual, tenthPredicted, isNew: true);
|
||||
}
|
||||
|
||||
// Remember state after 10 values
|
||||
double stateAfterTen = pseudoHuber.Last.Value;
|
||||
|
||||
// Generate 9 corrections with isNew=false (different values)
|
||||
for (int i = 0; i < 9; i++)
|
||||
{
|
||||
var bar = gbm.Next(isNew: false);
|
||||
pseudoHuber.Update(bar.Close, bar.Close * 1.01, isNew: false);
|
||||
}
|
||||
|
||||
// Feed the remembered 10th input again with isNew=false
|
||||
var finalResult = pseudoHuber.Update(tenthActual, tenthPredicted, isNew: false);
|
||||
|
||||
// State should match the original state after 10 values
|
||||
Assert.Equal(stateAfterTen, finalResult.Value, Epsilon);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Reset_ClearsState()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(DefaultPeriod);
|
||||
|
||||
pseudoHuber.Update(100.0, 99.0);
|
||||
pseudoHuber.Update(101.0, 99.0);
|
||||
double valueBefore = pseudoHuber.Last.Value;
|
||||
|
||||
pseudoHuber.Reset();
|
||||
|
||||
Assert.Equal(0, pseudoHuber.Last.Value);
|
||||
Assert.False(pseudoHuber.IsHot);
|
||||
|
||||
pseudoHuber.Update(50.0, 49.0);
|
||||
Assert.NotEqual(0, pseudoHuber.Last.Value);
|
||||
Assert.NotEqual(valueBefore, pseudoHuber.Last.Value);
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
#region Robustness Tests
|
||||
|
||||
[Fact]
|
||||
public void NaN_Input_UsesLastValidValue()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(5);
|
||||
|
||||
pseudoHuber.Update(100.0, 99.0);
|
||||
pseudoHuber.Update(101.0, 99.5);
|
||||
|
||||
var resultAfterNaN = pseudoHuber.Update(double.NaN, 100.0);
|
||||
Assert.True(double.IsFinite(resultAfterNaN.Value));
|
||||
|
||||
var resultAfterNaN2 = pseudoHuber.Update(102.0, double.NaN);
|
||||
Assert.True(double.IsFinite(resultAfterNaN2.Value));
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Infinity_Input_UsesLastValidValue()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(5);
|
||||
|
||||
pseudoHuber.Update(100.0, 99.0);
|
||||
pseudoHuber.Update(101.0, 99.5);
|
||||
|
||||
var resultAfterPosInf = pseudoHuber.Update(double.PositiveInfinity, 100.0);
|
||||
Assert.True(double.IsFinite(resultAfterPosInf.Value));
|
||||
|
||||
var resultAfterNegInf = pseudoHuber.Update(102.0, double.NegativeInfinity);
|
||||
Assert.True(double.IsFinite(resultAfterNegInf.Value));
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
#region Batch/Span Tests
|
||||
|
||||
[Fact]
|
||||
public void BatchCalc_MatchesIterativeCalc()
|
||||
{
|
||||
var gbm = new GBM();
|
||||
var actualSeries = new TSeries();
|
||||
var predictedSeries = new TSeries();
|
||||
|
||||
const int count = 100;
|
||||
for (int i = 0; i < count; i++)
|
||||
{
|
||||
var bar = gbm.Next();
|
||||
actualSeries.Add(bar.Time, bar.Close);
|
||||
predictedSeries.Add(bar.Time, bar.Close * (1.0 + (i % 2 == 0 ? 0.01 : -0.01)));
|
||||
}
|
||||
|
||||
// Calculate iteratively
|
||||
var iterative = new PseudoHuber(DefaultPeriod);
|
||||
var iterativeResults = new List<double>();
|
||||
for (int i = 0; i < count; i++)
|
||||
{
|
||||
iterativeResults.Add(iterative.Update(actualSeries[i].Value, predictedSeries[i].Value).Value);
|
||||
}
|
||||
|
||||
// Calculate batch
|
||||
var batchResults = PseudoHuber.Calculate(actualSeries, predictedSeries, DefaultPeriod);
|
||||
|
||||
// Compare
|
||||
Assert.Equal(iterativeResults.Count, batchResults.Count);
|
||||
for (int i = 0; i < batchResults.Count; i++)
|
||||
{
|
||||
Assert.Equal(batchResults[i].Value, iterativeResults[i], Epsilon);
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void SpanBatch_ValidatesInput()
|
||||
{
|
||||
double[] actual = [1, 2, 3, 4, 5];
|
||||
double[] predicted = [1.1, 2.1, 3.1, 4.1, 5.1];
|
||||
double[] output = new double[5];
|
||||
double[] wrongSizeOutput = new double[3];
|
||||
|
||||
Assert.Throws<ArgumentException>(() =>
|
||||
PseudoHuber.Batch(actual.AsSpan(), predicted.AsSpan(), wrongSizeOutput.AsSpan(), DefaultPeriod));
|
||||
|
||||
Assert.Throws<ArgumentException>(() =>
|
||||
PseudoHuber.Batch(actual.AsSpan(), predicted.AsSpan(), output.AsSpan(), 0));
|
||||
|
||||
Assert.Throws<ArgumentException>(() =>
|
||||
PseudoHuber.Batch(actual.AsSpan(), predicted.AsSpan(), output.AsSpan(), DefaultPeriod, 0));
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void SpanBatch_MatchesTSeriesBatch()
|
||||
{
|
||||
var gbm = new GBM();
|
||||
var actualSeries = new TSeries();
|
||||
var predictedSeries = new TSeries();
|
||||
double[] actualData = new double[100];
|
||||
double[] predictedData = new double[100];
|
||||
double[] output = new double[100];
|
||||
|
||||
for (int i = 0; i < 100; i++)
|
||||
{
|
||||
var bar = gbm.Next();
|
||||
actualData[i] = bar.Close;
|
||||
predictedData[i] = bar.Close * 0.99;
|
||||
actualSeries.Add(bar.Time, actualData[i]);
|
||||
predictedSeries.Add(bar.Time, predictedData[i]);
|
||||
}
|
||||
|
||||
var tseriesResult = PseudoHuber.Calculate(actualSeries, predictedSeries, DefaultPeriod);
|
||||
PseudoHuber.Batch(actualData.AsSpan(), predictedData.AsSpan(), output.AsSpan(), DefaultPeriod);
|
||||
|
||||
for (int i = 0; i < 100; i++)
|
||||
{
|
||||
Assert.Equal(tseriesResult[i].Value, output[i], Epsilon);
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void SpanBatch_HandlesNaN()
|
||||
{
|
||||
double[] actual = [100, 110, double.NaN, 120, 130];
|
||||
double[] predicted = [99, 109, 115, double.NaN, 129];
|
||||
double[] output = new double[5];
|
||||
|
||||
PseudoHuber.Batch(actual.AsSpan(), predicted.AsSpan(), output.AsSpan(), 3);
|
||||
|
||||
foreach (var val in output)
|
||||
{
|
||||
Assert.True(double.IsFinite(val), $"Expected finite value but got {val}");
|
||||
}
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
#region Error Handling Tests
|
||||
|
||||
[Fact]
|
||||
public void Update_ThrowsOnSingleInput()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(DefaultPeriod);
|
||||
var input = new TValue(DateTime.UtcNow, 100.0);
|
||||
|
||||
Assert.Throws<NotSupportedException>(() => pseudoHuber.Update(input));
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Prime_ThrowsNotSupported()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(DefaultPeriod);
|
||||
double[] data = [1, 2, 3, 4, 5];
|
||||
|
||||
Assert.Throws<NotSupportedException>(() => pseudoHuber.Prime(data.AsSpan()));
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Calculate_MismatchedSeriesLengths_Throws()
|
||||
{
|
||||
var actual = new TSeries();
|
||||
var predicted = new TSeries();
|
||||
|
||||
actual.Add(DateTime.UtcNow, 100);
|
||||
actual.Add(DateTime.UtcNow, 101);
|
||||
predicted.Add(DateTime.UtcNow, 99);
|
||||
|
||||
Assert.Throws<ArgumentException>(() => PseudoHuber.Calculate(actual, predicted, 5));
|
||||
}
|
||||
|
||||
#endregion
|
||||
|
||||
#region Resync Tests
|
||||
|
||||
[Fact]
|
||||
public void Resync_PreventsFloatingPointDrift()
|
||||
{
|
||||
var pseudoHuber = new PseudoHuber(10);
|
||||
var gbm = new GBM();
|
||||
|
||||
// Feed many values to trigger resync
|
||||
for (int i = 0; i < 2500; i++)
|
||||
{
|
||||
var bar = gbm.Next();
|
||||
pseudoHuber.Update(bar.Close, bar.Close * 0.99);
|
||||
}
|
||||
|
||||
// Should still produce valid results after many iterations
|
||||
Assert.True(double.IsFinite(pseudoHuber.Last.Value));
|
||||
Assert.True(pseudoHuber.Last.Value >= 0);
|
||||
}
|
||||
|
||||
#endregion
|
||||
}
|
||||
@@ -0,0 +1,258 @@
|
||||
using System.Runtime.CompilerServices;
|
||||
using System.Runtime.InteropServices;
|
||||
|
||||
namespace QuanTAlib;
|
||||
|
||||
/// <summary>
|
||||
/// PseudoHuber: Pseudo-Huber Loss (Charbonnier Loss)
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// The Pseudo-Huber loss is a smooth approximation to the Huber loss function.
|
||||
/// Unlike Huber loss which has a piecewise definition, Pseudo-Huber is smooth
|
||||
/// and differentiable everywhere, making it ideal for gradient-based optimization.
|
||||
///
|
||||
/// Formula:
|
||||
/// PseudoHuber = δ² * (√(1 + (error/δ)²) - 1)
|
||||
///
|
||||
/// Key properties:
|
||||
/// - Smooth and continuously differentiable everywhere
|
||||
/// - Approximates L2 (squared error) for small errors
|
||||
/// - Approximates L1 (absolute error) for large errors
|
||||
/// - δ (delta) controls the transition point
|
||||
/// - More computationally efficient than Huber's conditional logic
|
||||
/// - Also known as Charbonnier loss in image processing
|
||||
/// </remarks>
|
||||
[SkipLocalsInit]
|
||||
public sealed class PseudoHuber : AbstractBase
|
||||
{
|
||||
private readonly RingBuffer _lossBuffer;
|
||||
private readonly double _delta;
|
||||
private readonly double _deltaSquared;
|
||||
|
||||
[StructLayout(LayoutKind.Auto)]
|
||||
private record struct State(double LossSum, double LastValidActual, double LastValidPredicted, int TickCount);
|
||||
private State _state;
|
||||
private State _p_state;
|
||||
|
||||
private const int ResyncInterval = 1000;
|
||||
private const double DefaultDelta = 1.0;
|
||||
|
||||
public PseudoHuber(int period, double delta = DefaultDelta)
|
||||
{
|
||||
if (period <= 0)
|
||||
throw new ArgumentException("Period must be greater than 0", nameof(period));
|
||||
if (delta <= 0)
|
||||
throw new ArgumentException("Delta must be positive", nameof(delta));
|
||||
|
||||
_lossBuffer = new RingBuffer(period);
|
||||
_delta = delta;
|
||||
_deltaSquared = delta * delta;
|
||||
Name = $"PseudoHuber({period},{delta:F3})";
|
||||
WarmupPeriod = period;
|
||||
}
|
||||
|
||||
public double Delta => _delta;
|
||||
public override bool IsHot => _lossBuffer.IsFull;
|
||||
|
||||
/// <summary>
|
||||
/// Computes Pseudo-Huber loss: δ² * (√(1 + (x/δ)²) - 1)
|
||||
/// </summary>
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
private double PseudoHuberLoss(double x)
|
||||
{
|
||||
double ratio = x / _delta;
|
||||
double ratioSq = ratio * ratio;
|
||||
return _deltaSquared * (Math.Sqrt(1.0 + ratioSq) - 1.0);
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public TValue Update(TValue actual, TValue predicted, bool isNew = true)
|
||||
{
|
||||
double actualVal = actual.Value;
|
||||
double predictedVal = predicted.Value;
|
||||
|
||||
if (!double.IsFinite(actualVal))
|
||||
actualVal = double.IsFinite(_state.LastValidActual) ? _state.LastValidActual : 0.0;
|
||||
else
|
||||
_state.LastValidActual = actualVal;
|
||||
|
||||
if (!double.IsFinite(predictedVal))
|
||||
predictedVal = double.IsFinite(_state.LastValidPredicted) ? _state.LastValidPredicted : 0.0;
|
||||
else
|
||||
_state.LastValidPredicted = predictedVal;
|
||||
|
||||
double error = actualVal - predictedVal;
|
||||
double loss = PseudoHuberLoss(error);
|
||||
|
||||
if (isNew)
|
||||
{
|
||||
_p_state = _state;
|
||||
|
||||
double removedLoss = _lossBuffer.Count == _lossBuffer.Capacity ? _lossBuffer.Oldest : 0.0;
|
||||
_state.LossSum = _state.LossSum - removedLoss + loss;
|
||||
_lossBuffer.Add(loss);
|
||||
|
||||
_state.TickCount++;
|
||||
if (_lossBuffer.IsFull && _state.TickCount >= ResyncInterval)
|
||||
{
|
||||
_state.TickCount = 0;
|
||||
_state.LossSum = _lossBuffer.RecalculateSum();
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
_state = _p_state;
|
||||
|
||||
double removedLoss = _lossBuffer.Count == _lossBuffer.Capacity ? _lossBuffer.Oldest : 0.0;
|
||||
_state.LossSum = _state.LossSum - removedLoss + loss;
|
||||
_lossBuffer.UpdateNewest(loss);
|
||||
_state.LossSum = _lossBuffer.RecalculateSum();
|
||||
}
|
||||
|
||||
// Mean Pseudo-Huber Loss
|
||||
double result = _lossBuffer.Count > 0 ? _state.LossSum / _lossBuffer.Count : 0.0;
|
||||
|
||||
Last = new TValue(actual.Time, result);
|
||||
PubEvent(Last, isNew);
|
||||
return Last;
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public TValue Update(double actual, double predicted, bool isNew = true)
|
||||
{
|
||||
return Update(new TValue(DateTime.UtcNow, actual), new TValue(DateTime.UtcNow, predicted), isNew);
|
||||
}
|
||||
|
||||
public override TValue Update(TValue input, bool isNew = true)
|
||||
{
|
||||
throw new NotSupportedException("PseudoHuber requires two inputs. Use Update(actual, predicted).");
|
||||
}
|
||||
|
||||
public override TSeries Update(TSeries source)
|
||||
{
|
||||
throw new NotSupportedException("PseudoHuber requires two inputs. Use Calculate(actualSeries, predictedSeries, period, delta).");
|
||||
}
|
||||
|
||||
public override void Prime(ReadOnlySpan<double> source, TimeSpan? step = null)
|
||||
{
|
||||
throw new NotSupportedException("PseudoHuber requires two inputs.");
|
||||
}
|
||||
|
||||
public override void Reset()
|
||||
{
|
||||
_lossBuffer.Clear();
|
||||
_state = default;
|
||||
_p_state = default;
|
||||
Last = default;
|
||||
}
|
||||
|
||||
public static TSeries Calculate(TSeries actual, TSeries predicted, int period, double delta = DefaultDelta)
|
||||
{
|
||||
if (actual.Count != predicted.Count)
|
||||
throw new ArgumentException("Actual and predicted series must have the same length", nameof(predicted));
|
||||
|
||||
int len = actual.Count;
|
||||
var t = new List<long>(len);
|
||||
var v = new List<double>(len);
|
||||
CollectionsMarshal.SetCount(t, len);
|
||||
CollectionsMarshal.SetCount(v, len);
|
||||
|
||||
var tSpan = CollectionsMarshal.AsSpan(t);
|
||||
var vSpan = CollectionsMarshal.AsSpan(v);
|
||||
|
||||
Batch(actual.Values, predicted.Values, vSpan, period, delta);
|
||||
actual.Times.CopyTo(tSpan);
|
||||
|
||||
return new TSeries(t, v);
|
||||
}
|
||||
|
||||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||||
public static void Batch(ReadOnlySpan<double> actual, ReadOnlySpan<double> predicted, Span<double> output, int period, double delta = DefaultDelta)
|
||||
{
|
||||
if (actual.Length != predicted.Length || actual.Length != output.Length)
|
||||
throw new ArgumentException("All spans must have the same length", nameof(output));
|
||||
if (period <= 0)
|
||||
throw new ArgumentException("Period must be greater than 0", nameof(period));
|
||||
if (delta <= 0)
|
||||
throw new ArgumentException("Delta must be positive", nameof(delta));
|
||||
|
||||
int len = actual.Length;
|
||||
if (len == 0) return;
|
||||
|
||||
double deltaSquared = delta * delta;
|
||||
|
||||
const int StackAllocThreshold = 256;
|
||||
Span<double> lossBuffer = period <= StackAllocThreshold
|
||||
? stackalloc double[period]
|
||||
: new double[period];
|
||||
|
||||
double lossSum = 0;
|
||||
double lastValidActual = 0;
|
||||
double lastValidPredicted = 0;
|
||||
|
||||
for (int k = 0; k < len; k++)
|
||||
{
|
||||
if (double.IsFinite(actual[k])) { lastValidActual = actual[k]; break; }
|
||||
}
|
||||
for (int k = 0; k < len; k++)
|
||||
{
|
||||
if (double.IsFinite(predicted[k])) { lastValidPredicted = predicted[k]; break; }
|
||||
}
|
||||
|
||||
int bufferIndex = 0;
|
||||
int i = 0;
|
||||
|
||||
int warmupEnd = Math.Min(period, len);
|
||||
for (; i < warmupEnd; i++)
|
||||
{
|
||||
double act = actual[i];
|
||||
double pred = predicted[i];
|
||||
|
||||
if (double.IsFinite(act)) lastValidActual = act; else act = lastValidActual;
|
||||
if (double.IsFinite(pred)) lastValidPredicted = pred; else pred = lastValidPredicted;
|
||||
|
||||
double error = act - pred;
|
||||
double ratio = error / delta;
|
||||
double ratioSq = ratio * ratio;
|
||||
double loss = deltaSquared * (Math.Sqrt(1.0 + ratioSq) - 1.0);
|
||||
|
||||
lossSum += loss;
|
||||
lossBuffer[i] = loss;
|
||||
|
||||
output[i] = lossSum / (i + 1);
|
||||
}
|
||||
|
||||
int tickCount = 0;
|
||||
for (; i < len; i++)
|
||||
{
|
||||
double act = actual[i];
|
||||
double pred = predicted[i];
|
||||
|
||||
if (double.IsFinite(act)) lastValidActual = act; else act = lastValidActual;
|
||||
if (double.IsFinite(pred)) lastValidPredicted = pred; else pred = lastValidPredicted;
|
||||
|
||||
double error = act - pred;
|
||||
double ratio = error / delta;
|
||||
double ratioSq = ratio * ratio;
|
||||
double loss = deltaSquared * (Math.Sqrt(1.0 + ratioSq) - 1.0);
|
||||
|
||||
lossSum = lossSum - lossBuffer[bufferIndex] + loss;
|
||||
lossBuffer[bufferIndex] = loss;
|
||||
|
||||
bufferIndex++;
|
||||
if (bufferIndex >= period) bufferIndex = 0;
|
||||
|
||||
output[i] = lossSum / period;
|
||||
|
||||
tickCount++;
|
||||
if (tickCount >= ResyncInterval)
|
||||
{
|
||||
tickCount = 0;
|
||||
double recalcSum = 0;
|
||||
for (int k = 0; k < period; k++)
|
||||
recalcSum += lossBuffer[k];
|
||||
lossSum = recalcSum;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,156 @@
|
||||
# Pseudo-Huber: Smooth Huber Approximation
|
||||
|
||||
> "All the robustness of Huber, none of the discontinuities."
|
||||
|
||||
Pseudo-Huber Loss (also called Charbonnier Loss) is a smooth approximation to the Huber loss function. Unlike Huber which has a piecewise definition with a kink at δ, Pseudo-Huber is continuously differentiable everywhere, making it ideal for gradient-based optimization.
|
||||
|
||||
## Historical Context
|
||||
|
||||
The Pseudo-Huber function emerged from the optimization and machine learning communities as a way to get Huber-like robustness while maintaining smooth gradients. It's also known as Charbonnier loss in image processing, where it's used for edge-preserving smoothing and optical flow estimation.
|
||||
|
||||
## Architecture & Physics
|
||||
|
||||
Pseudo-Huber uses the formula δ²(√(1 + (x/δ)²) - 1), which smoothly interpolates between quadratic behavior for small errors and linear behavior for large errors. The transition is gradual rather than abrupt, with no discontinuity in derivatives.
|
||||
|
||||
### Properties
|
||||
|
||||
- **Smooth everywhere**: Infinitely differentiable (unlike Huber's kink)
|
||||
- **Non-negative**: Always ≥ 0, with 0 for perfect prediction
|
||||
- **Robust**: Large errors grow linearly, not quadratically
|
||||
- **Tunable**: δ (delta) controls the L2-to-L1 transition point
|
||||
|
||||
## Mathematical Foundation
|
||||
|
||||
### 1. Pseudo-Huber Function
|
||||
|
||||
For each error, compute:
|
||||
|
||||
$$L_\delta(e) = \delta^2 \left(\sqrt{1 + \left(\frac{e}{\delta}\right)^2} - 1\right)$$
|
||||
|
||||
Where:
|
||||
- $e = y - \hat{y}$ = prediction error
|
||||
- $\delta$ = tuning parameter (transition width)
|
||||
|
||||
### 2. Asymptotic Behavior
|
||||
|
||||
For small errors (|e| << δ):
|
||||
|
||||
$$L_\delta(e) \approx \frac{e^2}{2}$$
|
||||
|
||||
For large errors (|e| >> δ):
|
||||
|
||||
$$L_\delta(e) \approx \delta|e| - \delta^2$$
|
||||
|
||||
### 3. Gradient (Derivative)
|
||||
|
||||
$$\frac{dL}{de} = \frac{e}{\sqrt{1 + (e/\delta)^2}}$$
|
||||
|
||||
This approaches:
|
||||
- e for small errors (like L2)
|
||||
- δ·sign(e) for large errors (like L1)
|
||||
|
||||
### 4. Running Update (O(1))
|
||||
|
||||
QuanTAlib uses a ring buffer with running sum for O(1) updates:
|
||||
|
||||
$$S_{new} = S_{old} - L_{oldest} + L_{newest}$$
|
||||
|
||||
$$PseudoHuber = \frac{S_{new}}{n}$$
|
||||
|
||||
## Implementation Details
|
||||
|
||||
### Usage Patterns
|
||||
|
||||
```csharp
|
||||
// Streaming mode - with custom delta
|
||||
var pseudoHuber = new PseudoHuber(period: 20, delta: 1.0);
|
||||
var result = pseudoHuber.Update(actualValue, predictedValue);
|
||||
|
||||
// Batch mode - calculate for entire series
|
||||
var results = PseudoHuber.Calculate(actualSeries, predictedSeries, period: 20, delta: 1.0);
|
||||
|
||||
// Span mode - zero-allocation for high performance
|
||||
PseudoHuber.Batch(actualSpan, predictedSpan, outputSpan, period: 20, delta: 1.0);
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Default | Description |
|
||||
| :--- | :--- | :--- | :--- |
|
||||
| **period** | int | - | Lookback window for averaging (must be > 0) |
|
||||
| **delta** | double | 1.0 | Transition parameter (must be > 0) |
|
||||
|
||||
### Properties
|
||||
|
||||
| Property | Type | Description |
|
||||
| :--- | :--- | :--- |
|
||||
| **Last** | TValue | Most recent Pseudo-Huber value |
|
||||
| **IsHot** | bool | True when buffer is full |
|
||||
| **Delta** | double | Current delta parameter |
|
||||
| **Name** | string | Indicator name (e.g., "PseudoHuber(20,1.000)") |
|
||||
| **WarmupPeriod** | int | Number of periods before valid output |
|
||||
|
||||
## Performance Profile
|
||||
|
||||
| Metric | Score | Notes |
|
||||
| :--- | :--- | :--- |
|
||||
| **Throughput** | ~15 ns/bar | O(1) update, sqrt computation |
|
||||
| **Allocations** | 0 | Uses pre-allocated ring buffer |
|
||||
| **Complexity** | O(1) | Constant time per update |
|
||||
| **Accuracy** | 10/10 | Exact calculation |
|
||||
| **Timeliness** | 9/10 | No lag beyond the period |
|
||||
| **Smoothness** | 10/10 | Infinitely differentiable |
|
||||
|
||||
## Comparison with Huber
|
||||
|
||||
| Aspect | Huber | Pseudo-Huber |
|
||||
| :--- | :--- | :--- |
|
||||
| **Small errors** | e²/2 | ≈ e²/2 |
|
||||
| **Large errors** | δ\|e\| - δ²/2 | ≈ δ\|e\| - δ² |
|
||||
| **At e = δ** | Kink (C¹) | Smooth (C^∞) |
|
||||
| **Gradient** | Discontinuous 2nd derivative | Continuous all derivatives |
|
||||
| **Computation** | Conditional logic | Single formula |
|
||||
| **Optimization** | Can cause issues | Smooth convergence |
|
||||
|
||||
### Numerical Comparison
|
||||
|
||||
| Error (e) | Huber (δ=1) | Pseudo-Huber (δ=1) |
|
||||
| :--- | :--- | :--- |
|
||||
| **0.0** | 0.000 | 0.000 |
|
||||
| **0.5** | 0.125 | 0.118 |
|
||||
| **1.0** | 0.500 | 0.414 |
|
||||
| **2.0** | 1.500 | 1.236 |
|
||||
| **10.0** | 9.500 | 9.049 |
|
||||
|
||||
Pseudo-Huber produces slightly smaller values but follows the same qualitative behavior.
|
||||
|
||||
## Choosing δ
|
||||
|
||||
| δ Value | Behavior | Use Case |
|
||||
| :--- | :--- | :--- |
|
||||
| **0.1** | Quickly linear | Aggressive outlier handling |
|
||||
| **1.0** | Balanced | Standard choice |
|
||||
| **10.0** | Mostly quadratic | Near-MSE behavior |
|
||||
| **100.0** | Almost pure L2 | When outliers are rare |
|
||||
|
||||
## Common Use Cases
|
||||
|
||||
1. **Neural Network Training**: Smooth loss for gradient descent
|
||||
2. **Computer Vision**: Optical flow, stereo matching
|
||||
3. **Robust Regression**: When smoothness matters for optimization
|
||||
4. **Image Processing**: Edge-preserving filtering (Charbonnier)
|
||||
|
||||
## Edge Cases
|
||||
|
||||
- **Perfect Predictions**: Returns exactly 0
|
||||
- **NaN Handling**: Uses last valid value substitution
|
||||
- **Single Input**: Not supported (requires two series)
|
||||
- **δ = 0**: Invalid (division by zero)
|
||||
- **Large Errors**: Numerically stable (no overflow)
|
||||
|
||||
## Related Indicators
|
||||
|
||||
- [Huber](../huber/Huber.md) - Huber Loss (piecewise, with kink)
|
||||
- [LogCosh](../logcosh/LogCosh.md) - Log-Cosh Loss (different smooth approximation)
|
||||
- [MAE](../mae/Mae.md) - Mean Absolute Error (pure L1)
|
||||
- [MSE](../mse/Mse.md) - Mean Squared Error (pure L2)
|
||||
Reference in New Issue
Block a user