code Welch's t-test in C++ using approximations for student's t values I need to code a Welch's t-test between two populations in C++ without using external libraries like, for example, boost.
I know that given my two populations of size $N_1$ and $N_2$, I can calculate the T stat like this: 
$$ T = \frac {{\bar X_1} - {\bar X_2}}{\sqrt {\frac {s_1^2}{N_1} + \frac {s_2^2}{N_2} }}$$
and approximate the degrees of fredom $\nu$ like this:
$$ \nu \approx \frac {(\frac {s_1^2}{N_1} + \frac {s_2^2}{N_2})^2}{\frac {(\frac{s_1^2}{N1})^2}{N_1-1} + \frac {\frac{s_2^2}{N2})^2}{N_2-1}} $$
After calculating these two, I know I can do the $t$-$test$ like this: $$Refuse \ (\bar X_1 = \bar X_2) \ if: \  T > t_{\frac{\alpha}{2}, \nu-1}$$
$$Accept \ (\bar X_1 = \bar X_2) \ if: \  T \le t_{\frac{\alpha}{2}, \nu-1}$$
Where $t_{\frac{\alpha}{2}, \nu-1}$ is defined as the value for which we have: $$ P(T_{\nu-1} \ge t_{\frac{\alpha}{2}, \nu-1}) = \frac{\alpha}{2} $$
Where $T_{\nu-1}$ is a Student's t distribution with $\nu-1$ degrees of freedom.
My teacher gave me this link for the approximation formulae (it's a linkto the "Abramowitz and Stegun"'s "Handbook of Mathematical Functions" book): http://people.math.sfu.ca/~cbm/aands/page_949.htm
which, however, I think i didn't understand well.
The teacher told me to use this formula from the aforementioned link: $$ t_p \sim x_p + \frac {g_1(x_p)}{\nu} + \frac {g_2(x_p)}{\nu^2} +\frac {g_1(x_p)}{\nu^3} + \frac {g_1(x_p)}{\nu^4}$$ 
With: $A(t_p|\nu) = 1-2p$, $Q(x_p) = p$ and $g_1(x), g_2(x) ...$ defined as polynomial right below.
It also specifies (page 927 of the book: http://people.math.sfu.ca/~cbm/aands/page_927.htm) that in general it denotes:
with: 
$F(X) = P$ the c.d.f. of the random variable $X$
$Q(X) = 1-P$ the "upper tail area" of $X$
$A(X) = P-Q$ the c.d.f. of $|X|$
But still I can't understand what represent $t_p$, $x_p$ in the formula.
I thought that $t_p$ was the value such that 
$$ P(T_{\nu} \ge t_{p,\nu}) = p $$
and that $x_p$ was a value for which approximation is specified some pages backward (page 933 of the book) with the formula:
$$x_p = t - \frac{c_0+c_1t+c_2t^2}{1+d_1t+d_2t^2+d_3t^3} + \epsilon(p)$$
with:
$t=\sqrt\frac{1}{p^2}$
$c_i, d_i$ fixed constants
$|\epsilon(p)| < 4,5 \times 10^{-4}$
But I'm not sure about that, because this formula for $x_p$ is written in the part which is supposed to deal with gaussian distributions
Anyway, I can't understand what represents $\epsilon(p)$... And so how to calculate the values.
I hope that some of you can help because I'm really stucked in this.
 A: As you know, a $t$-test has two phases. First, you calculate the $t$ statistic from the your data. Then, you make an accept/reject decision by comparing your calculated $t$ value with a "critical" $t$-value that depends on your significance threshold $\alpha$ and the data's degrees of freedom $\nu$.  You seem to have the first part--calculating the $t$ statistic--well in hand, so all you need to do now is find the critical $t$-value and compare your calculated value with. 
How do you do that?
You could start with the $t$ distribution's cumulative distribution function (cdf), which gives you $p$ as a function of $T$ (and $\nu$), invert it, and then plug in the largest value of $p$ you're willing to tolerate. That value is either $\alpha$ or $\frac{\alpha}{2}$, depending on whether your test has one or two tails.
Unfortunately, the degrees of freedom parameter makes it hard to invert the $t$ distribution's CDF. When $\nu$ is 1, 2, or 4, there are simple formulas for the inverse, but in the general case you are stuck with approximating it using a power series to asymptotically approximate it. 
Several formula exist for this (see Hill, 1970). The one you are using works like this:


*

*Start by passing your $\frac{\alpha}{2}$ value though the inverse cdf of a standard normal distribution (i.e., N(0,1)) to get the quantity the Handbook calls $x_p$. There are several approximations, starting on Page 933, including the one you transcribed above. That formula comes with an error bound, but is only valid for $0 \lt p \le 0.5$; try plugging in $p=0.1$ and $p=0.9$; the answer should be be around $\pm 1.2816$, but 0.9 is considerably less accurate. Fortunately, the normal distribution is symmetric around zero, which lets you replace $p>0.5$ with $1-p$. 

*Next, adjust this factor by adding to terms from the expansion (i.e., $\frac{g_1(x_p)}{\nu} + \frac{g_2(x_p)}{\nu^2} + \cdots$). The approximation error will obviously decrease as you add more terms.

*If you somehow needed to indicate the error in this (or other approximations), it is customary to call that quantity $\epsilon$; it's used the same way in the inverse normal cdf. 
The inverse cdf is sometimes called the quantile function (which is why these functions all start with q in r). The notation in the Handbook is a little confusing. $Q$ could mean any quantile function, but here it specifically indicates the standard normal one.  The equations for $x_p$ starting on page 933 tell you how to approximate that (note that they're only valid for $0 \lt p \le 0.5$ but the normal distribution is symmetric!). From here, it should just be a matter of plugging and chugging. 
Incidentally, MATLAB falls back on this method when $\nu \ge 1000$, but uses a more direct method for smaller degrees of freedom. 

PS: There are several different forms of the t-test, but the general idea is the same. The numerator of the t-statistic tells you how far apart the means are (either from each other, in a two sample test, or from some pre-determined value [often, but not always zero] in a one-sample test). The denominator tells you how precisely the means are known (as in a standard error). 

The t-statistic is therefore small when the means are close together, relative to their uncertainty and becomes larger when means are farther apart and/or are more precisely known. 
