[SciPy-user] Welch's ttest
Angus McMorland
amcmorl at gmail.com
Mon Aug 13 16:24:20 EDT 2007
On 14/08/07, Peter PootyTang <peter.pootytang at gmail.com> wrote:
> Hello,
>
> I am converting some code from 'R' to using scipy, and the ttest
> results aren't matching. I figured out that the reason is because of
> the versions of the ttest that are being used. In 'R' I am doing a
> ttest without the assumtion that the deviations are the same for both
> samples. However, I can't find the same test in scipy. tt_ind might
> work, except my basline and sample sets have different sizes.
>
> Does anyone know how to implement to Welch's ttest using scipy?
>
> http://en.wikipedia.org/wiki/Welch%27s_t_test
You'll need to check my working, but I did have a go at implementing
this sometime back from Sokal and Rohlf's Biometry text. I'm sorry the
code is pretty ugly. Please let me know if you decide the
implementation is wrong, so I can patch it up. Ideally I (or someone
else) should create some unittests for this.
def welchs_approximate_ttest(n1, mean1, sem1, \
n2, mean2, sem2, alpha):
# calculate standard variances of the means
svm1 = sem1**2 * n1
svm2 = sem2**2 * n2
print "standard variance of the mean 1: %0.4f" % svm1
print "standard variance of the mean 2: %0.4f" % svm2
print ""
t_s_prime = (mean1 - mean2)/n.sqrt(svm1/n1+svm2/n2)
print "t'_s = %0.4f" % t_s_prime
print ""
t_alpha_df1 = scipy.stats.t.ppf(1-alpha/2, n1 - 1)
t_alpha_df2 = scipy.stats.t.ppf(1-alpha/2, n2 - 1)
print "t_alpha[%d] = %0.4f" % (n1-1, t_alpha_df1)
print "t_alpha[%d] = %0.4f" % (n2-1, t_alpha_df2)
print ""
t_alpha_prime = (t_alpha_df1 * sem1**2 + t_alpha_df2 * sem2**2) / \
(sem1**2 + sem2**2)
print "t'_alpha = %0.4f" % t_alpha_prime
print ""
if abs(t_s_prime) > t_alpha_prime:
print "Significantly different"
return True
else:
print "Not significantly different"
return False
Angus.
--
AJC McMorland, PhD Student
Physiology, University of Auckland
More information about the SciPy-User
mailing list