On-line evaluation of mean and variance in C#

Question

In various projects, I have to evaluate the mean and/or the variance of relatively large samples.

I wrote the following, to help me evaluate these quantities with a constant footprint. Basically, it adds elements one per one and updates the variables.

namespace DataStructures
{
    /// <summary>
    /// Holds a "ghost sample" into memory. It has a constant memory print
    /// and updates the size, the average and the variance of the sample.
    /// </summary>
    public interface IGhostSample<T>
    {
        void Add(T element);

        T Mean { get; }

        T Variance { get; }

        T StandardDev { get; }
    }
}

Which I specialized for doubles :

using System;

namespace DataStructures
{
    /// <summary>
    /// Holds a "ghost sample" into memory. It has a constant memory print
    /// and updates the size, the average and the variance of the sample.
    /// </summary>
    public class GhostSample : IGhostSample<double>
    {
        #region Private Attributes
        private double _mean = 0;
        private double _variance = 0;
        private int _size = 0;
        #endregion

        #region Accessors
        public int Size
        {
            get { return _size; }
        }

        public double Mean
        {
            get { return _mean; }
        }

        public double Variance
        {
            get { return _variance; }
        }

        public double StandardDev
        {
            get { return Math.Sqrt(_variance); }
        }
        #endregion

        #region Methods
        public void Add(double element)
        {
            double previousMean = _mean;
            _mean = (previousMean * _size + element) / (_size + 1);
            _variance = (_size * _variance + (element - previousMean) * (element - _mean)) / (_size + 1);
            _size++;
        }
        #endregion
    }
}

And Math.Net Vectors :

using MathNet.Numerics.LinearAlgebra.Double;
using System;

namespace DataStructures
{
    /// <summary>
    /// Holds a "ghost sample" into memory. It has a constant memory print
    /// and updates the size, the average and the variance of the sample.
    /// </summary>
    public class VectorGhostSample : IGhostSample<DenseVector>
    {
        #region Private attributes
        private int _size = 0;
        private int _length;
        private DenseVector _mean;
        private DenseVector _variance;
        #endregion

        #region Constructor
        /// <summary>
        /// Builds a new ghost sample.
        /// </summary>
        /// <param name="length">The length of the vectors.</param>
        public VectorGhostSample(int length)
        {
            _length = length;
            _mean = new DenseVector(_length);
            _variance = new DenseVector(_length);
        }
        #endregion

        #region Accessors
        /// <summary>
        /// The number of elements of the sample.
        /// </summary>
        public int Size
        {
            get { return _size; }
        }

        /// <summary>
        /// The (element-wise) mean of the sample
        /// </summary>
        public DenseVector Mean
        {
            get { return _mean; }
        }

        /// <summary>
        /// The (element-wise) variance of the sample
        /// </summary>
        public DenseVector Variance
        {
            get { return _variance; }
        }

        /// <summary>
        /// The (element-wise) standard deviation of the sample
        /// </summary>
        public DenseVector StandardDev
        {
            get
            {
                DenseVector std = new DenseVector(_variance.Count);
                for (int i = 0; i < _variance.Count; i++)
                    std[i] = Math.Sqrt(_variance[i]);
                return std;
            }
        }
        #endregion

        #region Methods
        /// <summary>
        /// Adds an element to the ghost sample (i.e. updates the mean and variance of the sample).
        /// </summary>
        /// <param name="element">The element to add to the sample</param>
        public void Add(DenseVector element)
        {
            DenseVector previousMean = _mean;
            _mean = (previousMean * _size + element) / (_size + 1);
            _variance = (DenseVector)(_size * _variance + (element - previousMean).PointwiseMultiply(element - _mean)).Divide(_size + 1);
            _size++;
        }
        #endregion
    }
}

I could not find a shorter way to do this. But vector and double behave in the same fashion, I only need addition and multiplication by a scalar to be able to perform these operations. Did I miss an obvious shortcut ?

Answer 1

Unfortunately there is no interface that defines the operators but you could define your own data type and with a litte over-engineering you get his:

Base class for the argument:

public abstract class Argument<T>
{
    protected Argument(T value)
    {
        Value = value;
    }

    public T Value { get; }

    public abstract Argument<T> Add(Argument<T> y);
    public abstract Argument<T> Subtract(Argument<T> y);
    public abstract Argument<T> Multiply(Argument<T> y);
    public abstract Argument<T> Divide(Argument<T> y);
    public abstract Argument<T> Increment();
    public abstract Argument<T> Sqrt();
}

Example double argument:

public class DoubleArgument : Argument<double>
{
    public static DoubleArgument Default = new DoubleArgument(0d);  
    public DoubleArgument(double value) : base(value) { }
    public override Argument<double> Add(Argument<double> y) => (DoubleArgument)(Value + y.Value);
    public override Argument<double> Subtract(Argument<double> y) => (DoubleArgument)(Value - y.Value);
    public override Argument<double> Multiply(Argument<double> y) => (DoubleArgument)(Value * y.Value);
    public override Argument<double> Divide(Argument<double> y) => (DoubleArgument)(Value / y.Value);
    public override Argument<double> Increment() => (DoubleArgument)(Value + 1.0);
    public override Argument<double> Sqrt() => (DoubleArgument)(Math.Sqrt(Value));
    public static implicit operator DoubleArgument(double value) => new DoubleArgument(value);
}

Calculation:

public class GhostSample<T>
{
    private Argument<T> _mean;
    private Argument<T> _variance;
    private Argument<T> _size;

    public GhostSample(Argument<T> size, Argument<T> mean, Argument<T> variance) 
    {
        _size = size;
        _mean = mean;
        _variance = variance;
    }

    public Argument<T> Size => _size;
    public Argument<T> Mean => _mean;
    public Argument<T> Variance => _variance;
    public Argument<T> StandardDev => _variance.Sqrt();    

    public void Add(Argument<T> element)
    {
        _size = _size.Increment();        
        var previousMean = _mean;
        _mean = previousMean.Multiply(_size).Add(element).Divide(_size); 
        _variance = _size.Multiply(_variance).Add(element.Subtract(previousMean).Multiply(element.Subtract(_mean)).Divide(_size));
    }
}

Usage:

var doubleGhostSample = new GhostSample<double>(DoubleArgument.Default, DoubleArgument.Default, DoubleArgument.Default);
doubleGhostSample.Add((DoubleArgument)2.0);
doubleGhostSample.Add((DoubleArgument)3.0);
doubleGhostSample.Add((DoubleArgument)8.0);

I didn't know how to get rid of the default-parameters :-(

Answer 2

Unfortunatly, Reference Types (classes – heap allocated) are very bad choice for the Math. Using them untolerabably slow downs calculations and significantly increases resource consumtions. We could define our own struct type, but it looks like a some kind of trade off – it is not 100% type safe. Let’s define an algorithm using a generic struct Value<T>:

public class GhostSample<T>
{
    public Value<T> Size { get; private set; }
    public Value<T> Mean { get; private set; }
    public Value<T> Variance { get; private set; }
    public Value<T> StandardDev => Variance.Sqrt;

    public void Add(Value<T> element)
    {
        var previousMean = Mean;
        Mean = (previousMean * Size + element) / (Size + 1);
        Variance = (Size * Variance + (element - previousMean) * (element - Mean)) / ++Size;
    }
}

We could use it this way (please note the initialization code – it should be repeated for every data type supported):

class Program
{
    static Program()
    {
        Value<double>.AddOp = (l, r) => l + r;
        Value<double>.AddScalarOp = (l, r) => l + r;
        Value<double>.DivideOp = (l, r) => l / r;
        Value<double>.DivideScalarOp = (l, r) => l / r;
        Value<double>.MultiplyOp = (l, r) => l * r;
        Value<double>.MultiplyScalarOp = (l, r) => l * r;
        Value<double>.SubtractOp = (l, r) => l - r;
        Value<double>.SubtractScalarOp = (l, r) => l - r;
        Value<double>.SqrtOp = v => Math.Sqrt(v);
    }

    static void Main(string[] args)
    {
        var doubleGhostSample = new GhostSample<double>();
        doubleGhostSample.Add(1);
        doubleGhostSample.Add(3);
        doubleGhostSample.Add(8);
        Console.WriteLine(doubleGhostSample.Mean); // 4
    }
}

Where library class would be:

public struct Value<T>
{
    public delegate T UnaryOp(T value);
    public delegate T BinaryOp(T left, T right);
    public delegate T ScalarOp(T left, double right);

    public static BinaryOp AddOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static ScalarOp AddScalarOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static BinaryOp SubtractOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static ScalarOp SubtractScalarOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static BinaryOp MultiplyOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static ScalarOp MultiplyScalarOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static BinaryOp DivideOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static ScalarOp DivideScalarOp { get; set; } = (l, r) => { throw new NotSupportedException(); };
    public static UnaryOp SqrtOp { get; set; } = (v) => { throw new NotSupportedException(); };

    public static Value<T> operator+(Value<T> left, Value<T> right) => AddOp(left, right);
    public static Value<T> operator+(Value<T> left, double right) => AddScalarOp(left, right);
    public static Value<T> operator -(Value<T> left, Value<T> right) => SubtractOp(left, right);
    public static Value<T> operator -(Value<T> left, double right) => SubtractScalarOp(left, right);
    public static Value<T> operator *(Value<T> left, Value<T> right) => MultiplyOp(left, right);
    public static Value<T> operator *(Value<T> left, double right) => MultiplyScalarOp(left, right);        
    public static Value<T> operator /(Value<T> left, Value<T> right) => DivideOp(left, right);
    public static Value<T> operator /(Value<T> left, double right) => DivideScalarOp(left, right);
    public static Value<T> operator ++(Value<T> value) => AddScalarOp(value, 1.0);
    public static Value<T> operator --(Value<T> value) => SubtractScalarOp(value, 1.0);        

    public static implicit operator T(Value<T> value) => value._value;
    public static implicit operator Value<T>(T value) => new Value<T>(value);

    Value(T value)
    {
        _value = value;
    }

    public Value<T> Sqrt => SqrtOp(this);

    readonly T _value;
}

Hope this helps.