# 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.

• $\epsilon$ measures the error of approximation. It is meant to be neglected in calculations. The rest of your question asks about how to calculate formulas that you have very clearly set forth--so what exactly is the problem?
– whuber
Jul 8, 2016 at 18:50
• @whuber Thank you very much for your help. 1. as I wrote I can't figure out who's $x_p$. What it represents? 2. Is the $x_p$ for which I found the approximation I reported the same used in the formula for $t_p$ approximation?
– ela
Jul 8, 2016 at 19:20
• I should have asked this before answering, but is there a particular reason you're doing this? If you just need to analyze some data, there are easier ways--you could call R via Rcpp, call out to a library that implements the t-test (GSL, Cephes); I even think Boost now includes a t-test. On the other hand, if you're doing this to learn, it's not a bad exercise, but the Handbook is an awful source to learn from. It may be good as a desk reference (though in practice, no one writes their own inverse normal distribution these days), but it's way too pithy to learn from. Jul 11, 2016 at 8:45
• @MattKrause Actually yes, I'm trying to evaluate if a device is weak to Side Channel Attacks with the methods showed here: eprint.iacr.org/2015/207.pdf. I asked my teacher if I could use boost, but he said "no". He gave me the Handbook just to use formulae for approximation, the study from where I'm studyign statistics is the italian version of "Introduction to Probability and Statistics for Engineers and Scientists" written by Sheldon M. Ross P.S. I've answered to your comments below your answer asking for some more infos, did you take a look? =)
– ela
Jul 11, 2016 at 9:19
• @MattKrause by the way, for the first attempt I used boost for calculating the t-test; I followed this guide: boost.org/doc/libs/1_36_0/libs/math/doc/sf_and_dist/html/…. But comparing it with the methodology for the t-test used in the book where I'm studying from (Sheldon M. Ross), which is the one I described in my question, I noticed they're not the same. (follows)
– ela
Jul 11, 2016 at 9:28

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:

1. 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$.

2. 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.

3. 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 ). 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.

• Does the case $\nu=3$ have an "easy" closed form (i.e. not relying on inverse beta functions) for the inverse CDF? I was aware that it did for 1, 2 and 4. When $\nu=3$, the CDF contains both a rational function and an inverse tan, so I'd be surprised if there were a nice inverse. Eqns 24 and 27 of this piece by W. T. Shaw are relevant. The t CDF looks quite different in even and odd cases, because of the inverse tans in the odd cases, but when $\nu=1$ there's no rational bit so easy to invert. For 3, 5, 7...it looks tricky! Jul 9, 2016 at 0:42
• @Silverfish To my recollection you can just use integration by parts to get another such integral at lower df (IIRC it reduces by 2 each time) Jul 9, 2016 at 0:57
• @Glen_b The going-down-by-two makes sense (since the even and odd cases are so distinct) but I'm not sure how this helps when it's the CDF itself, rather than the PDF, which is so gruesome. I think I will ask this as a question in its own right, I hadn't intended to derail the comments thread (I was just taken by surprise by the presence of 3 in the list of easy cases!) Jul 9, 2016 at 1:18
• @Silverfish Uh, but the value cdf is an integral of a pdf, right? Jul 9, 2016 at 1:20
• @Silverfish, try the new link. Also, I think I was mistaken about the $\nu=3$ case: there's a reasonably short form for the density, but the distribution gets a little ugly and it's downhill from there. Sorry! Jul 9, 2016 at 1:51