4
\$\begingroup\$

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 ?

\$\endgroup\$
  • 1
    \$\begingroup\$ Unfortunately, there is no type safe way to handle it other than that, I think. Generics are bad for math. A way worse, than, for example templates in C++. \$\endgroup\$ – Dmitry Nogin Jun 21 '16 at 8:49
3
\$\begingroup\$

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 :-(

\$\endgroup\$
  • \$\begingroup\$ Please see my answer. It could be a way to handle it more efficiently. \$\endgroup\$ – Dmitry Nogin Jun 24 '16 at 5:18
2
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ I started with the real operators like +*/- etc but after I got the same chaos in one class as you have ;-P I discarded it. It's a pity that those operators are static and there is no way to override them in a derived class. Even if you create dummies in the abstract class c# will still pick them in spite of defining them for the derived object too. It's really a mess ;-) I'd go with the Add or AddOp or whatever instead of trying to reuse the normal operators. \$\endgroup\$ – t3chb0t Jun 24 '16 at 5:24
  • \$\begingroup\$ @t3chb0t Let's say that mess starts when things are getting tangled. Wide API is not always a bad thing by itself. I would personally prefer a big stupid library class combined with precise, concise, and nice looking business logic as it is usually more volatile. Also, value types can not be inherited, so it is the only way. As you see this solution can adapt multiple types, just add Value<Matrix>.AddOp = (x,y) => x + y, etc. somewhere. You could put setup code into different classes, if necessary - it is not inheritance, but works the same way. \$\endgroup\$ – Dmitry Nogin Jun 24 '16 at 5:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.