/*3:*/
#line 100 "statlib.w"


#ifndef STATLIB_HEADER_DEFINES
#define STATLIB_HEADER_DEFINES
#include <math.h>        
#include <assert.h> 

#include <iostream> 
#include <exception> 
#include <stdexcept> 
#include <string> 
#include <vector> 
#include <algorithm> 
using namespace std;
/*6:*/
#line 184 "statlib.w"


class probabilityDistribution{
private:
virtual string distributionName(void)= 0;

public:
static double Q_from_P(double x){
return 1-x;
}

static double P_from_Q(double x){
return 1-x;
}

virtual double mean(void)= 0;
virtual double stdev(void)= 0;
double variance(void){
return stdev()*stdev();
}
virtual double skewness(void)= 0;
virtual double kurtosisExcess(void)= 0;
double kurtosis(void){
return kurtosisExcess()+3;
}

virtual double CDF_P(double x)= 0;

double CDF_Q(double x){
return Q_from_P(CDF_P(x));
}


virtual void writeParameters(ostream&of){
of<<"Distribution: "<<distributionName()<<"\n";
of<<"    Mean = "<<mean()<<"  Stdev = "<<stdev()<<
"  Variance = "<<variance()<<"  Skewness = "<<
skewness()<<"\n";
of<<"    Kurtosis = "<<kurtosis()<<
"  KurtosisExcess = "<<kurtosisExcess()<<"\n";
}
};

/*:6*//*7:*/
#line 242 "statlib.w"


class normalDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"normal";
}

double mu;
double sigma;

public:


normalDistribution(){
set_mu_sigma();
}




normalDistribution(double c_mu,double c_sigma){
set_mu_sigma(c_mu,c_sigma);
}



void set_mu(double c_mu= 0){
mu= c_mu;
}

void set_sigma(double c_sigma= 1){
sigma= c_sigma;
}

void set_mu_sigma(double c_mu= 0,double c_sigma= 1){
set_mu(c_mu);
set_sigma(c_sigma);
}

double get_mu(void){
return mu;
}

double get_sigma(void){
return sigma;
}



double mean(void){
return get_mu();
}

double stdev(void){
return get_sigma();
}

double skewness(void){
return 0;
}

double kurtosisExcess(void){
return 0;
}



/*8:*/
#line 334 "statlib.w"

void set_variance(double var){
set_sigma(sqrt(var));
}

/*:8*/
#line 311 "statlib.w"
;
/*9:*/
#line 346 "statlib.w"

void approximateBinomial(unsigned int n)
{
set_mu_sigma(n/2,sqrt(n/4.0));
}

/*:9*/
#line 312 "statlib.w"
;



static double p_from_mu_sigma_x(double mu,double sigma,double x);
static double x_from_p_mu_sigma(double p,double mu,double sigma);
static double mu_from_p_x_sigma(double p,double x,double sigma);
static double sigma_from_p_x_mu(double p,double x,double mu);



/*10:*/
#line 358 "statlib.w"

double z_score(double m){
return(m-mean())/stdev();
}

/*:10*/
#line 323 "statlib.w"
;

double CDF_P(double x);
};

/*:7*//*13:*/
#line 467 "statlib.w"


class chiSquareDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"chi-square";
}

double k;

public:


chiSquareDistribution(){
set_k();
}

chiSquareDistribution(double c_k= 1){
set_k(c_k);
}



void set_k(double c_k= 1){
k= c_k;
}

double get_k(void){
return k;
}



double mean(void){
return get_k();
}

double stdev(void){
return sqrt(2*k);
}

double skewness(void){
return sqrt(8/k);
}

double kurtosisExcess(void){
return 12/k;
}



static double p_from_k_x(double k,double x);
static double x_from_p_k(double p,double k);
static double k_from_p_x(double p,double x);



double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    Degrees of freedom = "<<k<<"\n";
}
};

/*:13*//*16:*/
#line 613 "statlib.w"


class gammaDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"gamma";
}

double alpha;
double lambda;

public:


gammaDistribution(){
set_alpha_lambda();
}

gammaDistribution(double c_alpha= 10,double c_lambda= 10){
set_alpha_lambda(c_alpha,c_lambda);
}



void set_alpha(double c_alpha= 10){
alpha= c_alpha;
}

double get_alpha(void){
return alpha;
}

void set_lambda(double c_lambda= 10){
lambda= c_lambda;
}

double get_lambda(void){
return lambda;
}

void set_alpha_lambda(double c_alpha= 10,double c_lambda= 10){
set_alpha(c_alpha);
set_lambda(c_lambda);
}



double mean(void){
return alpha/lambda;
}

double stdev(void){
return sqrt(alpha)/lambda;
}

double skewness(void){
return 2/sqrt(alpha);
}

double kurtosisExcess(void){
return 6/alpha;
}



static double p_from_alpha_lambda_x(double alpha,double lambda,double x);
static double x_from_p_alpha_lambda(double p,double alpha,double lambda);
static double alpha_from_p_lambda_x(double p,double lambda,double x);
static double lambda_from_p_alpha_x(double p,double alpha,double x);

double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    Shape (alpha) = "<<alpha<<"  Scale (lambda)  = "<<lambda<<"\n";
}
};

/*:16*//*19:*/
#line 778 "statlib.w"


class betaDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"beta";
}

double a;
double b;

public:


betaDistribution(){
set_a_b();
}

betaDistribution(double c_a= 5,double c_b= 10){
set_a_b(c_a,c_b);
}



void set_a(double c_a= 5){
a= c_a;
}

double get_a(void){
return a;
}

void set_b(double c_b= 10){
b= c_b;
}

double get_b(void){
return b;
}

void set_a_b(double c_a= 5,double c_b= 10){
set_a(c_a);
set_b(c_b);
}



double mean(void){
return a/(a+b);
}

double stdev(void){
return sqrt(a*b)/((a+b)*sqrt(a+b+1));

}

double skewness(void){
return(2*(b-a)*sqrt(a+b+1))/(sqrt(a*b)*(a+b+2));

}

double kurtosisExcess(void){
return((3*(a+b+1)*((a*b*(a+b-6))+
(2*(a+b)*(a+b)))/
(a*b*(a+b+2)*(a+b+3))))-3;


}



static double p_from_a_b_x(double a,double b,double x);
static double x_from_p_a_b(double p,double a,double b);
static double a_from_p_b_x(double p,double b,double x);
static double b_from_p_a_x(double p,double a,double x);

double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    a = "<<a<<"  b = "<<b<<"\n";
}
};

/*:19*//*22:*/
#line 957 "statlib.w"


class tDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"t";
}

double k;

public:


tDistribution(){
set_k();
}

tDistribution(double c_k= 3){
set_k(c_k);
}



void set_k(double c_k= 3){
k= c_k;
}

double get_k(void){
return k;
}



double mean(void){
return 0;
}

double stdev(void){
return sqrt(k/(k-2));
}

double skewness(void){
return 0;
}

double kurtosisExcess(void){
return 6.0/(k-4);
}



static double p_from_k_x(double k,double x);
static double x_from_p_k(double p,double k);
static double k_from_p_x(double p,double x);



double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    Degrees of freedom = "<<k<<"\n";
}
};

/*:22*//*25:*/
#line 1101 "statlib.w"


class FDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"F";
}

double u;
double v;

public:


FDistribution(){
set_u_v();
}

FDistribution(double c_u= 8,double c_v= 8){
set_u_v(c_u,c_v);
}



void set_u(double c_u= 8){
u= c_u;
}

double get_u(void){
return u;
}

void set_v(double c_v= 8){
v= c_v;
}

double get_v(void){
return v;
}

void set_u_v(double c_u= 8,double c_v= 8){
set_u(c_u);
set_v(c_v);
}



double mean(void){
return v/(v-2);
}

double stdev(void){
return sqrt((2*(v*v)*(u+v-2))/(u*((v-2)*(v-2))*(v-4)));

}

double skewness(void){
return(2*sqrt(2.0)*sqrt(v-4)*(2*u+v-2))/
(sqrt(u)*(v-6)*sqrt(u+v-2));

}

double kurtosisExcess(void){
return(12*((v-4)*(v-2)*(v-2)+
u*(u+v-2)*(5*v-22)))/
(u*(v-8)*(v-6)*(u+v-2));


}



static double p_from_u_v_x(double u,double v,double x);
static double x_from_p_u_v(double p,double u,double v);
static double u_from_p_v_x(double p,double v,double x);
static double v_from_p_u_x(double p,double u,double x);

double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    Numerator DF = "<<u<<"  Denominator DF  = "<<v<<"\n";
}
};

/*:25*//*28:*/
#line 1269 "statlib.w"


class poissonDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"poisson";
}

double lambda;

public:


poissonDistribution(){
set_lambda();
}

poissonDistribution(double c_lambda= 0.5){
set_lambda(c_lambda);
}



void set_lambda(double c_lambda= 0.5){
lambda= c_lambda;
}

double get_lambda(void){
return lambda;
}



double mean(void){
return lambda;
}

double stdev(void){
return sqrt(lambda);
}

double skewness(void){
return 1/sqrt(lambda);
}

double kurtosisExcess(void){
return 1/lambda;
}



static double p_from_lambda_x(double lambda,double x);
static double x_from_p_lambda(double p,double lambda);
static double lambda_from_p_x(double p,double x);



double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    Arrival rate = "<<lambda<<"\n";
}
};

/*:28*//*31:*/
#line 1420 "statlib.w"


class binomialDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"binomial";
}

double n;
double r;

public:


binomialDistribution(){
set_n_r();
}

binomialDistribution(double c_n= 2,double c_r= 0.5){
set_n_r(c_n,c_r);
}



void set_n(double c_n= 2){
n= c_n;
}

double get_n(void){
return n;
}

void set_r(double c_r= 0.5){
r= c_r;
}

double get_r(void){
return r;
}

void set_n_r(double c_n= 2,double c_r= 0.5){
set_n(c_n);
set_r(c_r);
}



double mean(void){
return n*r;
}

double stdev(void){
return sqrt(n*r*(1-r));
}

double skewness(void){
return(1-(2*r))/stdev();
}

double kurtosisExcess(void){
return(1-(6*(1-r)*r))/(n*(1-r)*r);

}



static double p_from_n_r_s(double n,double r,double s);
static double s_from_p_r_n(double p,double r,double n);
static double n_from_p_r_s(double p,double r,double s);
static double r_from_p_n_s(double p,double n,double s);



/*10:*/
#line 358 "statlib.w"

double z_score(double m){
return(m-mean())/stdev();
}

/*:10*/
#line 1494 "statlib.w"
;

double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    Trials = "<<n<<"  Probability  = "<<r<<"\n";
}
};

/*:31*//*34:*/
#line 1600 "statlib.w"


class negativeBinomialDistribution:public probabilityDistribution{

private:
virtual string distributionName(void){
return"negative binomial";
}

double n;
double r;

public:


negativeBinomialDistribution(){
set_n_r();
}

negativeBinomialDistribution(double c_n= 2,double c_r= 0.5){
set_n_r(c_n,c_r);
}



void set_n(double c_n= 2){
n= c_n;
}

double get_n(void){
return n;
}

void set_r(double c_r= 0.5){
r= c_r;
}

double get_r(void){
return r;
}

void set_n_r(double c_n= 2,double c_r= 0.5){
set_n(c_n);
set_r(c_r);
}



double mean(void){
return(n*(1-r))/r;
}

double stdev(void){
return sqrt(n*(1-r))/r;
}

double skewness(void){
return(2-r)/sqrt(n*(1-r));
}

double kurtosisExcess(void){
return(6*(1-r)+(r*r))/(n*(1-r));

}



static double p_from_n_r_s(double n,double r,double s);
static double s_from_p_r_n(double p,double r,double n);
static double n_from_p_r_s(double p,double r,double s);
static double r_from_p_n_s(double p,double n,double s);

double CDF_P(double x);


void writeParameters(ostream&of){
probabilityDistribution::writeParameters(of);
of<<"    Successes = "<<n<<"  Probability  = "<<r<<"\n";
}
};

/*:34*//*37:*/
#line 1769 "statlib.w"

template<class T> class dataTable:public vector<T> 
{
public:
double mean(void);
double geometricMean(void);
double harmonicMean(void);
double median(void);
double RMS(void);
T mode(void);
double percentile(double k);
double quartile(int q)
{
assert(q>=1&&q<=3);
return percentile(0.25*q);
}

double variance(void);
double varianceMLE(void)
{
return(variance()*(this->size()-1))/this->size();
}
double stdev(void)
{
return sqrt(variance());
}
double stdevMLE(void)
{
return sqrt(varianceMLE());
}

double centralMoment(int k);
double skewness(void);
double kurtosis(void);
};

/*:37*//*38:*/
#line 1813 "statlib.w"

template<class T> double dataTable<T> ::mean(void)
{
typename dataTable<T> ::iterator p;
double m= 0;

for(p= this->begin();p!=this->end();p++){
m+= *p;
}
return m/this->size();
}

/*:38*//*39:*/
#line 1832 "statlib.w"

template<class T> double dataTable<T> ::geometricMean(void)
{
typename dataTable<T> ::iterator p;
double g= 1,ni= 1.0/this->size();

for(p= this->begin();p!=this->end();p++){
g*= pow(*p,ni);
}
return g;
}

/*:39*//*40:*/
#line 1851 "statlib.w"

template<class T> double dataTable<T> ::harmonicMean(void)
{
typename dataTable<T> ::iterator p;
double d= 0;

for(p= this->begin();p!=this->end();p++){
d+= 1.0/(*p);
}
return this->size()/d;
}

/*:40*//*41:*/
#line 1870 "statlib.w"

template<class T> double dataTable<T> ::RMS(void)
{
typename dataTable<T> ::iterator p;
double sum= 0;

for(p= this->begin();p!=this->end();p++){
sum+= (*p)*(*p);
}
return sqrt(sum/this->size());
}

/*:41*//*42:*/
#line 1894 "statlib.w"

template<class T> double dataTable<T> ::median(void)
{
dataTable<T> v(*this);

sort(v.begin(),v.end());
if(v.size()&1){
return(double)v[(v.size()+1)/2];
}else{
return(double)((v[(v.size()/2)-1]+v[v.size()/2])/2.0);
}
}

/*:42*//*43:*/
#line 1923 "statlib.w"

template<class T> T dataTable<T> ::mode(void)
{
dataTable<T> v(*this);
T cmode,bmode= 0;
int n= 0,most= 0;
typename dataTable<T> ::iterator p;

sort(v.begin(),v.end());
p= v.begin();
cmode= *p++;
while(p!=v.end()){
T cval= *p++;
if(cval==cmode){
n++;
}else{
if(n> most){
most= n;
bmode= cmode;
}
n= 1;
cmode= cval;
}
}
if(n> most){
most= n;
bmode= cmode;
}
return bmode;
}

/*:43*//*44:*/
#line 1963 "statlib.w"

template<class T> double dataTable<T> ::percentile(double k)
{
dataTable<T> v(*this);
double index,result;
int i;

assert(k>=0&&k<=1);

sort(v.begin(),v.end());

index= v.size()*k;
if(index!=floor(index)){
i= ((int)index);
result= v[i];
}else{
i= ((int)index)-1;
result= (v[i]+v[i+1])/2.0;
}
return result;
}

/*:44*//*45:*/
#line 1994 "statlib.w"

template<class T> double dataTable<T> ::variance(void)
{
typename dataTable<T> ::iterator p;
double mu= mean(),sum= 0;

for(p= this->begin();p!=this->end();p++){
double t= ((*p)-mu);
sum+= t*t;
}
return sum/(this->size()-1);
}

/*:45*//*46:*/
#line 2016 "statlib.w"

template<class T> double dataTable<T> ::centralMoment(int k)
{
typename dataTable<T> ::iterator p;
double mu= mean(),sum= 0;

for(p= this->begin();p!=this->end();p++){
sum+= pow((*p)-mu,(double)k);
}
return sum/this->size();
}

/*:46*//*47:*/
#line 2037 "statlib.w"

template<class T> double dataTable<T> ::skewness(void)
{
double sigma= stdevMLE();

return centralMoment(3)/(sigma*sigma*sigma);
}

/*:47*//*48:*/
#line 2055 "statlib.w"

template<class T> double dataTable<T> ::kurtosis(void)
{
double v= varianceMLE();

return centralMoment(4)/(v*v);
}

/*:48*/
#line 114 "statlib.w"

#endif

/*:3*/