3
$\begingroup$

Trying to test code creating P-value manually against SciPy. The Scipy Documentation isn't the best, which makes it tought to know for sure what to do.

I am getting the correct t-stat and P-value with SciPy, but I'm not able to replicate the correct p-value manually - A friend steered me to scipy.stats.t.ppf - but I'm not getting a p-value from it.

What is the correct way to do scipy.stats.t.ppf()?

my version:

def t_test(sample, mu):
    mean = np.mean(sample)
    var = np.var(sample)
    sem = (var / len(sample)) ** .5
    t = abs(mu - mean)/sem
    df = len(sample) - 1
    p = scs.t.ppf(.95, df)
    return (t, p)

returns (0.081500599630942958, 1.7291328115213671)

scipy version:
scs.ttest_1samp(sample, 4.123)
returns (statistic=0.079436958358141435, pvalue=0.93751577779749051)

for testing, I'm using the following sample set and sample mu.

sample = [4.15848606,  3.86146363,  4.31545726,  3.3748772,
          4.67023082,  4.45950272,  3.85894915,  4.41089417,
          3.82360986,  3.79889443,  4.75884172,  3.27100914,
          4.08939402,  4.08904694,  5.62589842,  3.71445656,
          3.58463792,  4.42426443,  3.9671448 ,  4.39339124]

mu = 4.123 
$\endgroup$
2
  • $\begingroup$ I am unable to reproduce these results. For instance, the mean of the sample is 4.132523, which hardly differs from mu at all, so I cannot obtain anything like these statistics. What exactly is the distinction between mu and mean in your code? $\endgroup$
    – whuber
    Apr 1 '17 at 18:26
  • $\begingroup$ @whuber mu is a made up estimation of the population mean - while mean is the mean of the sample found here in sample I just copied and pasted the code into my terminal and it ran correctly - I edited the above to run my version within a function so you could run that easily $\endgroup$
    – be-ns
    Apr 1 '17 at 19:34
3
$\begingroup$

To get the same results, change two things:

  1. Change the estimation of the variance such that the divisor is N-1
  2. Calculate the p-value using the cdf, that is the probability of getting a more extreme value, here using that the t-distribution is symmetric around zero. Note that the function you're comparing with does a two-sided test, and therefore, so do I.

I've marked the relevant lines with ###. The result is now the same as from the ttest_1samp function.

def t_test(sample, mu):
mean = np.mean(sample)
var = np.var(sample, ddof = 1) ###
sem = (var / len(sample)) ** .5
t = abs(mu - mean)/sem
df = len(sample) - 1
p = 2*(1-scs.t.cdf(t, df)) ###
return (t, p)
$\endgroup$
4
  • $\begingroup$ worked like a charm. Just adding the definition of ddof for any future visitors. “Delta Degrees of Freedom”: the divisor used in the calculation is N - ddof, where N represents the number of elements. By default ddof is zero." $\endgroup$
    – be-ns
    Apr 1 '17 at 23:07
  • $\begingroup$ I think the p value should be 2*min(1-scs.t.cdf(t, df), scs.t.cdf(t, df)) - depending on which side of the t-distribution you are. $\endgroup$
    – reox
    Dec 2 '20 at 12:47
  • $\begingroup$ @reox: That would be correct if 't' were the signed t-statistic, but note that in swmo's code, it is the absolute value. This ensures that we are always on the right side of the distribution. $\endgroup$ Dec 2 '20 at 13:48
  • $\begingroup$ @RubenvanBergen ah yes you are right, I missed the abs. $\endgroup$
    – reox
    Dec 3 '20 at 14:18
0
$\begingroup$

Adding to the previous posts, p-value can be computed with either of the following two lines of code, for a 2-sided t-test, indicating whether the extreme value (of the test statistic) is falling in the left or the right critical region (with -ve and +ve t-statistic values, respectively), as shown in the next figure.

2*(t.cdf(-abs(t_stat), dof))
2*(1 - t.cdf(abs(t_stat), dof))

where t-stat is the value of the t-statistic and dof is the associated degrees of freedom.

from scipy.stats import t
t_stat = 2.25
dof = 15
# p-value for 2-sided test
2*(1 - t.cdf(abs(t_stat), dof))
# 0.03988800677091664
2*(t.cdf(-abs(t_stat), dof))
# 0.03988800677091648

The below figure shows how the critical region for 5% level of significance looks like for a 2-sided t-test. For the above example, we can see that the null hypothesis can be rejected.

enter image description here

$\endgroup$
4
  • $\begingroup$ Although implementation is often mixed with substantive content in questions, we are supposed to be a site for providing information about statistics, machine learning, etc., not code. It can be good to provide code as well, but please elaborate your substantive answer in text for people who don't read this language well enough to recognize & extract the answer from the code. $\endgroup$
    – Sycorax
    Sep 11 at 22:11
  • $\begingroup$ @Sycorax provided some information about hypothesis test $\endgroup$ Sep 11 at 22:16
  • $\begingroup$ That's all well and good, but it's not clear how this answers the question, which includes "how to compute the $t$ statistic" as one of its components. $\endgroup$
    – Sycorax
    Sep 11 at 22:24
  • $\begingroup$ @Sycorax I guess that part is already answered, just wanted to focus on different ways of computing the p-value that work correctly irrespective of the sign of the t-statistic (which is a part of the question too) $\endgroup$ Sep 11 at 22:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.