//
// kth central moments of a sequence in one pass 
//
// copyright (c) 2009 gordon anderson
// http://quantblog.wordpress.com
// released under BSD licence
//
#include <stdio.h>
#include <vector>
#include <boost/math/special_functions/binomial.hpp>        // our nCr impl

using namespace boost::math;
using namespace std;

#define _count_of(rg)  ((sizeof(rg))/(sizeof(rg[0])))

///

inline double nCr(int n, int r)
{
    assert(n<24);
    assert(r<=n);
    return binomial_coefficient<double>(n,r);               //todo - is this fast??
}

///

struct PrintVisitor
{
    void  operator()(double& v)
    {
        printf("val=%lf\n", v);
    };
    void  operator()(int& v)
    {
        printf("val=%d\n", v);
    };
};


// V value type of sequence elements
// R floating return type of stats functions
// C integral count type 
//
template<typename V=double, typename R=double, typename C=unsigned int>
struct StatsRecorder
{
    StatsRecorder(int ord=2)
    : order(std::max(ord,2)), moments(order+1, 0.)
    {
    }

    void  operator()(const V& x)
    {
        // accumulate simple kth powers

        V x_k=1.;
        for (int k=0;k<=order;k++)
        {
            moments[k]+=x_k;
            x_k*=x;
        }
    }

    R kth_moment(unsigned int k) const
    {
        assert(k<=order);

        int n=count();
        return (n>0) ? moments[k]/n : 0.;
    }

    R kth_central_moment(unsigned int k) const
    {
        // reconstruct central moment from k simple moments 
        assert(k<=order);

        R un=mean();
        R cmk=0.;
        for(int j=0;j<=k;j++)
        {
            cmk+=nCr(k,j)*pow(-un,(int)k-j)*moments[j];
        }
        int n=count();
        return (n>0) ? cmk/n : 0.;
    }

    R mean() const
    {
        int n=count();
        return (n>0) ? sum()/n : 0.;
    }
    R varn() const                             // population variance - denominator n
    {
        return kth_central_moment(2);
    }
    R var() const                              // sample variance - denominator n-1
    {
        int n=count();
        return (n>1) ? (n*varn()/(n-1)) : 0;
    }
    R sd() const { return sqrt(var()); }       // sample deviation
    R sdn() const { return sqrt(varn()); }     // population deviation

    C count() const { return moments[0]; }      // count is just the 0th simple moment
    R sum() const { return moments[1]; }        // the sum is just the 1st simple moment

    void print() const
    {
        printf("stats : sum=%lf mean=%lf var=%lf sd=%lf [varn=%lf sdn=%lf] n=%d\n", 
            (double)sum(), (double)mean(), 
            (double)var(), (double)sd(), 
            (double)varn(), (double)sdn(), count());

        for (int k=0;k<=order;k++)
        {
            printf("order=%d moment=%lf central_moment=%lf\n", k, kth_moment(k), kth_central_moment(k));
        }
    }

    int order;          // max k [we store this many moments]
    vector<R> moments;  // various kth order simple moments 
                        // nb. 0th is count, 1st is sum, 2nd is variance, (hence order kept >=2)
};


int main(void)
{
    const int order=6;

    double vals[]={0.5, 1.4, 3.2, 1.3, 1.0, 0.9, 2.7, 1.2};
    vector<double> v(&vals[0], &vals[0]+_count_of(vals));

    for_each(v.begin(), v.end(), PrintVisitor());    
    for_each(v.begin(), v.end(), StatsRecorder<>(order)).print();


    int ints[]={5, 4, 2, 3, 0, 9, 7, 2};
    vector<int> V(&ints[0], &ints[0]+_count_of(ints));

    for_each(V.begin(), V.end(), PrintVisitor());    
    for_each(V.begin(), V.end(), StatsRecorder<int>(order)).print();     

    return 0;
}