Quick implementation of the t-test in c++, compared with scipy

ttest_cpp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "require(['codemirror/mode/clike/clike'], function(Clike) { console.log('ROOTaaS - C++ CodeMirror module loaded'); });"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "IPython.CodeCell.config_defaults.highlight_modes['magic_text/x-c++src'] = {'reg':[/^%%cpp/]};"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to ROOTaaS 6.06/02\n"
     ]
    }
   ],
   "source": [
    "import ROOT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "    \n",
    "// c++ code, adapted from http://www.boost.org/doc/libs/1_36_0/libs/math/example/students_t_two_samples.cpp    \n",
    "    \n",
    "std::pair<double, double> two_samples_t_test_unequal_sd(\n",
    "        double Sm1,\n",
    "        double Sd1,\n",
    "        unsigned Sn1,\n",
    "        double Sm2,\n",
    "        double Sd2,\n",
    "        unsigned Sn2)\n",
    "{\n",
    "   //\n",
    "   // Sm1 = Sample Mean 1.\n",
    "   // Sd1 = Sample Standard Deviation 1.\n",
    "   // Sn1 = Sample Size 1.\n",
    "   // Sm2 = Sample Mean 2.\n",
    "   // Sd2 = Sample Standard Deviation 2.\n",
    "   // Sn2 = Sample Size 2.\n",
    "   //\n",
    "   // A Students t test applied to two sets of data.\n",
    "   // We are testing the null hypothesis that the two\n",
    "   // samples have the same mean and that any difference\n",
    "   // if due to chance.\n",
    "   // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda353.htm\n",
    "   //\n",
    "\n",
    "\n",
    "   // Print header:\n",
    "   cout <<\n",
    "      \"_________________________________________________\\n\"\n",
    "      \"Student t test for two samples (unequal variances)\\n\"\n",
    "      \"_________________________________________________\\n\\n\";\n",
    "   cout << setprecision(5);\n",
    "   cout << setw(55) << left << \"Number of Observations (Sample 1)\" << \"=  \" << Sn1 << \"\\n\";\n",
    "   cout << setw(55) << left << \"Sample 1 Mean\" << \"=  \" << Sm1 << \"\\n\";\n",
    "   cout << setw(55) << left << \"Sample 1 Standard Deviation\" << \"=  \" << Sd1 << \"\\n\";\n",
    "   cout << setw(55) << left << \"Number of Observations (Sample 2)\" << \"=  \" << Sn2 << \"\\n\";\n",
    "   cout << setw(55) << left << \"Sample 2 Mean\" << \"=  \" << Sm2 << \"\\n\";\n",
    "   cout << setw(55) << left << \"Sample 2 Standard Deviation\" << \"=  \" << Sd2 << \"\\n\";\n",
    "   //\n",
    "   // Now we can calculate and output some stats:\n",
    "   //\n",
    "   // Degrees of freedom:\n",
    "   double v = Sd1 * Sd1 / Sn1 + Sd2 * Sd2 / Sn2;\n",
    "   v *= v;\n",
    "   double t1 = Sd1 * Sd1 / Sn1;\n",
    "   t1 *= t1;\n",
    "   t1 /=  (Sn1 - 1);\n",
    "   double t2 = Sd2 * Sd2 / Sn2;\n",
    "   t2 *= t2;\n",
    "   t2 /= (Sn2 - 1);\n",
    "   v /= (t1 + t2);\n",
    "   cout << setw(55) << left << \"Degrees of Freedom\" << \"=  \" << v << \"\\n\";\n",
    "   // t-statistic:\n",
    "   double t_stat = (Sm1 - Sm2) / sqrt(Sd1 * Sd1 / Sn1 + Sd2 * Sd2 / Sn2);\n",
    "   cout << setw(55) << left << \"T Statistic\" << \"=  \" << t_stat << \"\\n\";\n",
    "   return make_pair(t_stat, v);\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ttest_indResult(statistic=1.0086894851058474, pvalue=0.34624784860876445)\n",
      "Ttest_indResult(statistic=1.0086894851058474, pvalue=0.3462478486087644)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "data1 = [1.1, 3, 4, 6, 7, 10., 0, 100.]\n",
    "data2 = [1.1, 2.9, 4.5, 5.9, 6.8]\n",
    "\n",
    "print stats.ttest_ind(data1, data2, equal_var=False)\n",
    "mean1, mean2 = np.mean(data1), np.mean(data2)\n",
    "std1, std2 = np.std(data1, ddof=1), np.std(data2, ddof=1)\n",
    "size1, size2 = len(data1), len(data2)\n",
    "print stats.stats.ttest_ind_from_stats(mean1, std1, size1, mean2, std2, size2, equal_var=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________\r\n",
      "Student t test for two samples (unequal variances)\r\n",
      "_________________________________________________\r\n",
      "\r\n",
      "Number of Observations (Sample 1)                      =  8\r\n",
      "Sample 1 Mean                                          =  16.387\r\n",
      "Sample 1 Standard Deviation                            =  33.939\r\n",
      "Number of Observations (Sample 2)                      =  5\r\n",
      "Sample 2 Mean                                          =  4.24\r\n",
      "Sample 2 Standard Deviation                            =  2.293\r\n",
      "Degrees of Freedom                                     =  7.102\r\n",
      "T Statistic                                            =  1.0087\r\n"
     ]
    }
   ],
   "source": [
    "# compute the ttest and the degree of freedom\n",
    "tresult = ROOT.two_samples_t_test_unequal_sd(mean1, std1, size1, mean2, std2, size2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.34624784860876456"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1 - ROOT.TMath.StudentI(np.abs(tresult.first), tresult.second)) * 2"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}