StudentT.h

Go to the documentation of this file.

// Honeycomb, Copyright (C) 2015 NewGamePlus Inc.  Distributed under the Boost Software License v1.0.
#pragma once

#include "Honey/Math/Random/Dist/ChiSqr.h"
#include "Honey/Math/Random/Dist/Beta.h"
#include "Honey/Math/Random/Random.h"

namespace honey
{


template<class Real>
class StudentT_ : public RandomDist<Real>
{
    typedef RandomDist<Real>    Super;
    RandomDist_imports();
    typedef Gaussian_<Double>   Gaussian;
    typedef Beta_<Double>       Beta;
    typedef BetaInc<Double>     BetaInc;
    typedef ChiSqr_<Double>     ChiSqr;
    typedef Random_<Real>       Random;
    
public:
    typedef typename Random::DistStats DistStats;
    typedef Vec<2,Real>     Vec2;

    StudentT_(optional<RandomGen&> gen, Real n)         : Super(gen), n(n) { assert(n > 0); }
    StudentT_(Real n)                                   : StudentT_(optnull, n) {}

    virtual Real next() const;
    virtual Real pdf(Real x) const;
    virtual Real cdf(Real x) const;
    virtual Real cdfComp(Real x) const;
    virtual Real cdfInv(Real P) const;
    virtual Real mean() const                           { return 0; }
    virtual Real variance() const                       { return n / (n - 2); }

    struct Stats
    {
        DistStats dist;         
        Vec2    meanCi;         
        Vec2    stdDevCi;       
        szt     df;             
        Real    alpha;          
        int     tail;           
        Real    t;              
        Real    p;              

        friend ostream& operator<<(ostream& os, const Stats& val)
        {
            int ci = 100*(1-val.alpha);
            return os           << stringstream::indentInc << "{" << endl
                                << "Dist:"              << endl << val.dist << endl
                                << "Mean CI "           << std::setw(2) << ci
                                <<           "%:    "   << val.meanCi << endl
                                << "Std Dev CI "        << std::setw(2) << ci
                                <<              "%: "   << val.stdDevCi << endl
                                << "DF:             "   << val.df << endl
                                << "t Value:        "   << val.t << endl
               << (val.tail == 0 ? "Pr > |t|:       " :
                   val.tail == 1 ? "Pr > t:         " :
                                   "Pr < t:         ")  << val.p 
                                << stringstream::indentDec << endl  << "}";
        }
    };


    template<class Range>
    static bool test(const Range& samples, optional<Stats&> stats = optnull, Real mu = 0, Real alpha = 0.05, int tail = 0)
    {
        assert(tail >= -1 && tail <= 1);
        DistStats d = Random::stats(samples);
        sdt df = d.n - 1;
        StudentT_ student(df);
        Real t = (d.mean - mu) / d.stdErr;
        Real p = student.cdfComp(tail == 0 ? Alge::abs(t) : tail == 1 ? t : -t);
        if (tail == 0) p *= 2;

        if (stats)
        {
            stats->dist = d;
            Real meanCi = student.cdfInv(1 - (tail == 0 ? alpha/2 : alpha)) * d.stdErr;
            stats->meanCi = Vec2(   tail == -1 ? -Real_::inf : d.mean - meanCi,
                                    tail == 1 ? Real_::inf : d.mean + meanCi);
            stats->stdDevCi = ChiSqr(df).stdDevCi(d.stdDev, alpha).template cast<Vec2>();
            stats->df = df;
            stats->alpha = alpha;
            stats->tail = tail;
            stats->t = t;
            stats->p = p;
        }

        return p <= alpha;
    }

    struct PooledStats
    {
        DistStats dist[2];      
        Real    mean;           
        Real    stdDev;         
        Real    stdErr;         
        Vec2    meanCi;         
        Vec2    stdDevCi;       
        szt     df;             
        Real    alpha;          
        int     tail;           
        Real    t;              
        Real    p;              

        friend ostream& operator<<(ostream& os, const PooledStats& val)
        {
            int ci = 100*(1-val.alpha);
            return os           << stringstream::indentInc << "{" << endl
                                << "Dist 1:"            << endl << val.dist[0] << endl
                                << "Dist 2:"            << endl << val.dist[1] << endl
                                << "Pooled Diff (1 - 2):" << endl
                                << stringstream::indentInc << "{" << endl
                                << "Mean:           "   << val.mean << endl
                                << "Std Dev:        "   << val.stdDev << endl
                                << "Std Err:        "   << val.stdErr << endl
                                << "Mean CI "           << std::setw(2) << ci
                                <<           "%:    "   << val.meanCi << endl
                                << "Std Dev CI "        << std::setw(2) << ci
                                <<              "%: "   << val.stdDevCi << endl
                                << "DF:             "   << val.df << endl
                                << "t Value:        "   << val.t << endl
               << (val.tail == 0 ? "Pr > |t|:       " :
                   val.tail == 1 ? "Pr > t:         " :
                                   "Pr < t:         ")  << val.p 
                                << stringstream::indentDec << endl  << "}"
                                << stringstream::indentDec << endl  << "}";
        }
    };


    template<class Range, class Range2>
    static typename std::enable_if<mt::isRange<Range2>::value, bool>::type
        test(const Range& samples1, const Range2& samples2, optional<PooledStats&> stats = optnull, Real mu = 0, Real alpha = 0.05, int tail = 0)
    {
        assert(tail >= -1 && tail <= 1);
        DistStats d1 = Random::stats(samples1);
        DistStats d2 = Random::stats(samples2);

        sdt df = d1.n + d2.n - 2;
        StudentT_ student(df);
        Real mean = d1.mean - d2.mean;
        Real stdDev = Alge::sqrt(((d1.n-1)*Alge::sqr(d1.stdDev) + (d2.n-1)*Alge::sqr(d2.stdDev)) / df);
        Real stdErr = stdDev * Alge::sqrt((1/Real(d1.n)) + (1/Real(d2.n)));
        Real t = (mean - mu) / stdErr;
        Real p = student.cdfComp(tail == 0 ? Alge::abs(t) : tail == 1 ? t : -t);
        if (tail == 0) p *= 2;

        if (stats)
        {
            stats->dist[0] = d1;
            stats->dist[1] = d2;
            stats->mean = mean;
            stats->stdDev = stdDev;
            stats->stdErr = stdErr;
            Real meanCi = student.cdfInv(1 - (tail == 0 ? alpha/2 : alpha)) * stdErr;
            stats->meanCi = Vec2(   tail == -1 ? -Real_::inf : mean - meanCi,
                                    tail == 1 ? Real_::inf : mean + meanCi);
            stats->stdDevCi = ChiSqr(df).stdDevCi(stdDev, alpha).template cast<Vec2>();
            stats->df = df;
            stats->alpha = alpha;
            stats->tail = tail;
            stats->t = t;
            stats->p = p;
        }

        return p <= alpha;
    }

    Real n;
};

typedef StudentT_<Real>     StudentT;
typedef StudentT_<Float>    StudentT_f;
typedef StudentT_<Double>   StudentT_d;

extern template class StudentT_<Float>;
extern template class StudentT_<Double>;

}