Welch's t-test: Difference between revisions
m →{{header|Tcl}}: added zkl header |
→{{header|zkl}}: Added code |
||
Line 745: | Line 745: | ||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|C}} |
|||
⚫ | |||
<lang zkl> |
<lang zkl>fcn calculate_Pvalue(array1,array2){ |
||
if (array1.len()<=1 or array2.len()<=1) return(1.0); |
|||
mean1,mean2 := array1.sum(0.0),array2.sum(0.0); |
|||
if(mean1==mean2) return(1.0); |
|||
mean1/=array1.len(); |
|||
mean2/=array2.len(); |
|||
variance1:=array1.reduce('wrap(sum,x){ sum + (x-mean1).pow(2) },0.0); |
|||
variance2:=array2.reduce('wrap(sum,x){ sum + (x-mean2).pow(2) },0.0); |
|||
variance1/=(array1.len() - 1); |
|||
variance2/=(array2.len() - 1); |
|||
WELCH_T_STATISTIC:=(mean1-mean2)/ |
|||
(variance1/array1.len() + variance2/array2.len()).sqrt(); |
|||
DEGREES_OF_FREEDOM:= |
|||
( variance1/array1.len() + variance2/array2.len() ).pow(2) // numerator |
|||
/ ( |
|||
(variance1*variance1)/(array1.len().pow(2)*(array1.len() - 1)) + |
|||
(variance2*variance2)/(array2.len().pow(2)*(array2.len() - 1)) |
|||
); |
|||
a:=DEGREES_OF_FREEDOM/2; |
|||
x:=DEGREES_OF_FREEDOM/( WELCH_T_STATISTIC.pow(2) + DEGREES_OF_FREEDOM ); |
|||
N,h := 65535, x/N; |
|||
sum1,sum2 := 0.0, 0.0; |
|||
foreach i in (N){ |
|||
sum1+=((h*i + h/2.0).pow(a - 1))/(1.0 - (h*i + h/2.0)).sqrt(); |
|||
sum2+=((h*i).pow(a - 1))/(1.0 - h*i).sqrt(); |
|||
} |
|||
return_value:=((h/6.0)*( x.pow(a - 1)/(1.0 - x).sqrt() + |
|||
4.0*sum1 + 2.0*sum2) ) / |
|||
((0.0).e.pow(lngammal(a) + 0.57236494292470009 - lngammal(a + 0.5))); |
|||
if(return_value > 1.0) return(1.0); // or return_value is infinite, throws |
|||
return_value; |
|||
} |
|||
fcn lngammal(xx){ |
|||
var [const] cof=List( // static |
|||
76.18009172947146, -86.50532032941677, |
|||
24.01409824083091, -1.231739572450155, |
|||
0.1208650973866179e-2,-0.5395239384953e-5 |
|||
); |
|||
y:=x:=xx; |
|||
tmp:=x + 5.5 - (x + 0.5) * (x + 5.5).log(); |
|||
ser:=1.000000000190015; |
|||
foreach x in (cof){ ser+=(x/(y+=1)); } |
|||
return((2.5066282746310005 * ser / x).log() - tmp); |
|||
}</lang> |
|||
⚫ | |||
T(T(27.5,21.0,19.0,23.6,17.0,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19.0,21.7,21.4), |
|||
T(27.1,22.0,20.8,23.4,23.4,23.5,25.8,22.0,24.8,20.2,21.9,22.1,22.9,20.5,24.4)), |
|||
T(T(17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8), |
|||
T(21.5,22.8,21.0,23.0,21.6,23.6,22.5,20.7,23.4,21.8,20.7,21.7,21.5,22.5,23.6,21.5,22.5,23.5,21.5,21.8)), |
|||
T(T(19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0), |
|||
T(28.2,26.6,20.1,23.3,25.2,22.1,17.7,27.6,20.6,13.7,23.2,17.5,20.6,18.0,23.9,21.6,24.3,20.4,24.0,13.2)), |
|||
T(T(30.02,29.99,30.11,29.97,30.01,29.99), |
|||
T(29.89,29.93,29.72,29.98,30.02,29.98)), |
|||
T(T(3.0,4.0,1.0,2.1),T(490.2,340.0,433.9)) ); |
|||
foreach x,y in (testSets) |
|||
{ println("Test set 1 p-value = %f".fmt(calculate_Pvalue(x,y))); }</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Test set 1 p-value = 0.021378 |
|||
Test set 1 p-value = 0.148842 |
|||
Test set 1 p-value = 0.035972 |
|||
Test set 1 p-value = 0.090773 |
|||
Test set 1 p-value = 0.010752 |
|||
</pre> |
</pre> |
Revision as of 15:30, 18 September 2017
Given two lists of data, calculate the p-Value used for null hypothesis testing.
Task Description
Given two sets of data, calculate the p-value:
x = {3.0,4.0,1.0,2.1} y = {490.2,340.0,433.9}
Your task is to discern whether or not the difference in means between the two sets is statistically significant and worth further investigation. P-values are significance tests to gauge the probability that the difference in means between two data sets is significant, or due to chance. A threshold level, alpha, is usually chosen, 0.01 or 0.05, where p-values below alpha are worth further investigation and p-values above alpha are considered not significant. The p-value is not considered a final test of significance, only whether the given variable should be given further consideration.
There is more than on way of calculating the t-statistic, and you must choose which method is appropriate for you. Here we use Welch's t-test, which assumes that the variances between the two sets x
and y
are not equal. Welch's t-test statistic can be computed:
where
is the mean of set ,
and
is the number of variables in set ,
and
is the square root of the unbiased sample variance of set , i.e.
and the degrees of freedom, can be approximated:
The two-tailed p-value, , can be computed as a cumulative distribution function
where I is the incomplete beta function. This is the same as:
Keeping in mind that
and
can be calculated in terms of gamma functions and integrals more simply:
which simplifies to
The definite integral can be approximated with Simpson's Rule but other methods are also acceptable.
The , or lgammal(x)
function is necessary for the program to work with large a
values, as Gamma functions can often return values larger than can be handled by double
or long double
data types. The lgammal(x)
function is standard in math.h
with C99 and C11 standards.
C
Link with -lm
This program, for example, pvalue.c, can be compiled by
clang -o pvalue pvalue.c -Wall -pedantic -std=c11 -lm -O2
or
gcc -o pvalue pvalue.c -Wall -pedantic -std=c11 -lm -O2
.
This shows how pvalue can be calculated from any two arrays, using Welch's 2-sided t-test, which doesn't assume equal variance.
Smaller p-values converge more quickly than larger p-values.
const unsigned short int N = 65535
ensures integral convergence of about for p-values < 0.15, about for p-values approximately 0.5, but only for p-values approaching 1. <lang C>#include <stdio.h>
- include <math.h>
double calculate_Pvalue (const double *array1, const size_t array1_size, const double *array2, const size_t array2_size) { if (array1_size <= 1) { return 1.0; } if (array2_size <= 1) { return 1.0; } double mean1 = 0.0, mean2 = 0.0; for (size_t x = 0; x < array1_size; x++) { mean1 += array1[x]; } for (size_t x = 0; x < array2_size; x++) { mean2 += array2[x]; } if (mean1 == mean2) { return 1.0; } mean1 /= array1_size; mean2 /= array2_size; double variance1 = 0.0, variance2 = 0.0; for (size_t x = 0; x < array1_size; x++) { variance1 += (array1[x]-mean1)*(array1[x]-mean1); } for (size_t x = 0; x < array2_size; x++) { variance2 += (array2[x]-mean2)*(array2[x]-mean2); } if ((variance1 == 0.0) && (variance2 == 0.0)) { return 1.0; } variance1 = variance1/(array1_size-1); variance2 = variance2/(array2_size-1); const double WELCH_T_STATISTIC = (mean1-mean2)/sqrt(variance1/array1_size+variance2/array2_size); const double DEGREES_OF_FREEDOM = pow((variance1/array1_size+variance2/array2_size),2.0)//numerator / ( (variance1*variance1)/(array1_size*array1_size*(array1_size-1))+ (variance2*variance2)/(array2_size*array2_size*(array2_size-1)) ); const double a = DEGREES_OF_FREEDOM/2, x = DEGREES_OF_FREEDOM/(WELCH_T_STATISTIC*WELCH_T_STATISTIC+DEGREES_OF_FREEDOM); const unsigned short int N = 65535; const double h = x/N; double sum1 = 0.0, sum2 = 0.0; for(unsigned short int i = 0;i < N; i++) {
sum1 += (pow(h * i + h / 2.0,a-1))/(sqrt(1-(h * i + h / 2.0))); sum2 += (pow(h * i,a-1))/(sqrt(1-h * i));
} double return_value = ((h / 6.0) * ((pow(x,a-1))/(sqrt(1-x)) + 4.0 * sum1 + 2.0 * sum2))/(expl(lgammal(a)+0.57236494292470009-lgammal(a+0.5))); if ((isfinite(return_value) == 0) || (return_value > 1.0)) { return 1.0; } else { return return_value; } }
int main(void) { const double d1[] = {27.5,21.0,19.0,23.6,17.0,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19.0,21.7,21.4}; const double d2[] = {27.1,22.0,20.8,23.4,23.4,23.5,25.8,22.0,24.8,20.2,21.9,22.1,22.9,20.5,24.4}; const double d3[] = {17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8}; const double d4[] = {21.5,22.8,21.0,23.0,21.6,23.6,22.5,20.7,23.4,21.8,20.7,21.7,21.5,22.5,23.6,21.5,22.5,23.5,21.5,21.8}; const double d5[] = {19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0}; const double d6[] = {28.2,26.6,20.1,23.3,25.2,22.1,17.7,27.6,20.6,13.7,23.2,17.5,20.6,18.0,23.9,21.6,24.3,20.4,24.0,13.2}; const double d7[] = {30.02,29.99,30.11,29.97,30.01,29.99}; const double d8[] = {29.89,29.93,29.72,29.98,30.02,29.98}; const double x[] = {3.0,4.0,1.0,2.1}; const double y[] = {490.2,340.0,433.9};
printf("Test sets 1 p-value = %lf\n",calculate_Pvalue(d1,sizeof(d1)/sizeof(*d1),d2,sizeof(d2)/sizeof(*d2))); printf("Test sets 2 p-value = %lf\n",calculate_Pvalue(d3,sizeof(d3)/sizeof(*d3),d4,sizeof(d4)/sizeof(*d4))); printf("Test sets 3 p-value = %lf\n",calculate_Pvalue(d5,sizeof(d5)/sizeof(*d5),d6,sizeof(d6)/sizeof(*d6))); printf("Test sets 4 p-value = %lf\n",calculate_Pvalue(d7,sizeof(d7)/sizeof(*d7),d8,sizeof(d8)/sizeof(*d8))); printf("Test sets 5 p-value = %lf\n",calculate_Pvalue(x,sizeof(x)/sizeof(*x),y,sizeof(y)/sizeof(*y))); return 0; } </lang>
- Output:
Test sets 1 p-value = 0.021378 Test sets 2 p-value = 0.148842 Test sets 3 p-value = 0.035972 Test sets 4 p-value = 0.090773 Test sets 5 p-value = 0.010751
If your computer does not have lgammal
, add the following function before main
and replace lgammal
with lngammal
in the calculate_Pvalue
function:
<lang C>#include <stdio.h>
- include <math.h>
long double lngammal(const double xx) {
unsigned int j; double x,y,tmp,ser; const double cof[6] = { 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5 }; y = x = xx; tmp = x + 5.5 - (x + 0.5) * logl(x + 5.5); ser = 1.000000000190015; for (j=0;j<=5;j++) ser += (cof[j] / ++y); return(log(2.5066282746310005 * ser / x) - tmp);
}
</lang>
Go
<lang go>package main
import (
"fmt" "math"
)
var (
d1 = []float64{27.5, 21.0, 19.0, 23.6, 17.0, 17.9, 16.9, 20.1, 21.9, 22.6, 23.1, 19.6, 19.0, 21.7, 21.4} d2 = []float64{27.1, 22.0, 20.8, 23.4, 23.4, 23.5, 25.8, 22.0, 24.8, 20.2, 21.9, 22.1, 22.9, 20.5, 24.4} d3 = []float64{17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8} d4 = []float64{21.5, 22.8, 21.0, 23.0, 21.6, 23.6, 22.5, 20.7, 23.4, 21.8, 20.7, 21.7, 21.5, 22.5, 23.6, 21.5, 22.5, 23.5, 21.5, 21.8} d5 = []float64{19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0} d6 = []float64{28.2, 26.6, 20.1, 23.3, 25.2, 22.1, 17.7, 27.6, 20.6, 13.7, 23.2, 17.5, 20.6, 18.0, 23.9, 21.6, 24.3, 20.4, 24.0, 13.2} d7 = []float64{30.02, 29.99, 30.11, 29.97, 30.01, 29.99} d8 = []float64{29.89, 29.93, 29.72, 29.98, 30.02, 29.98} x = []float64{3.0, 4.0, 1.0, 2.1} y = []float64{490.2, 340.0, 433.9}
)
func main() {
fmt.Printf("%.6f\n", pValue(d1, d2)) fmt.Printf("%.6f\n", pValue(d3, d4)) fmt.Printf("%.6f\n", pValue(d5, d6)) fmt.Printf("%.6f\n", pValue(d7, d8)) fmt.Printf("%.6f\n", pValue(x, y))
}
func mean(a []float64) float64 {
sum := 0. for _, x := range a { sum += x } return sum / float64(len(a))
}
func sv(a []float64) float64 {
m := mean(a) sum := 0. for _, x := range a { d := x - m sum += d * d } return sum / float64(len(a)-1)
}
func welch(a, b []float64) float64 {
return (mean(a) - mean(b)) / math.Sqrt(sv(a)/float64(len(a))+sv(b)/float64(len(b)))
}
func dof(a, b []float64) float64 {
sva := sv(a) svb := sv(b) n := sva/float64(len(a)) + svb/float64(len(b)) return n * n / (sva*sva/float64(len(a)*len(a)*(len(a)-1)) + svb*svb/float64(len(b)*len(b)*(len(b)-1)))
}
func simpson0(n int, upper float64, f func(float64) float64) float64 {
sum := 0. nf := float64(n) dx0 := upper / nf sum += f(0) * dx0 sum += f(dx0*.5) * dx0 * 4 x0 := dx0 for i := 1; i < n; i++ { x1 := float64(i+1) * upper / nf xmid := (x0 + x1) * .5 dx := x1 - x0 sum += f(x0) * dx * 2 sum += f(xmid) * dx * 4 x0 = x1 } return (sum + f(upper)*dx0) / 6
}
func pValue(a, b []float64) float64 {
ν := dof(a, b) t := welch(a, b) g1, _ := math.Lgamma(ν / 2) g2, _ := math.Lgamma(.5) g3, _ := math.Lgamma(ν/2 + .5) return simpson0(2000, ν/(t*t+ν), func(r float64) float64 { return math.Pow(r, ν/2-1) / math.Sqrt(1-r) }) / math.Exp(g1+g2-g3)
}</lang>
- Output:
0.021378 0.148842 0.035972 0.090773 0.010751
J
Implementation:
<lang J>integrate=: adverb define
'a b steps'=. 3{.y,128 size=. (b - a)%steps size * +/ u |: 2 ]\ a + size * i.>:steps
) simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'
lngamma=: ^.@!@<:`(^.@!@(1 | ]) + +/@:^.@(1 + 1&| + i.@<.)@<:)@.(1&<:)"0 mean=: +/ % # nu=: # - 1: sampvar=: +/@((- mean) ^ 2:) % nu ssem=: sampvar % # welch_T=: -&mean % 2 %: +&ssem nu=: nu f. : ((+&ssem ^ 2:) % +&((ssem^2:)%nu)) B=: ^@(+&lngamma - lngamma@+)
p2_tail=:dyad define
t=. x welch_T y NB. need numbers for numerical integration v=. x nu y F=. ^&(_1+v%2) % 2 %: 1&- lo=. 0 hi=. v%(t^2)+v (F f. simpson integrate lo,hi) % 0.5 B v%2
)</lang>
integrate
and simpson
are from the Numerical integration task.
lngamma
is from http://www.jsoftware.com/pipermail/programming/2015-July/042174.html -- for values less than some convenient threshold (we use 1, but we could use a modestly higher threshold), we calculate it directly. For larger values we compute the fractional part directly and rebuild the log of the factorial using the sum of the logs.
mean
is classic J - most J tutorials will include this
The initial definition of nu
(degrees of freedom of a data set), as well as the combining form (approximating degrees of freedom for two sets of data) is from Welch's t test. (Verb definitions can be forward referenced, even in J's tacit definitions, but it seems clearer to specify these definitions so they only depend on previously declared definitions.)
sampvar
is sample variance (or: standard deviation squared)
ssem
is squared standard error of the mean
Also... please ignore the highlighting of v
in the definition of p2_tail. In this case, it's F that's the verb, v is just another number (the degrees of freedom for our two data sets. (But this is a hint that in explicit conjunction definitions, v would be the right verb argument. Unfortunately, the wiki's highlighting implementation is not capable of distinguishing that particular context from other contexts.)
Data for task examples: <lang J>d1=: 27.5 21 19 23.6 17 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19 21.7 21.4 d2=: 27.1 22 20.8 23.4 23.4 23.5 25.8 22 24.8 20.2 21.9 22.1 22.9 20.5 24.4 d3=: 17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8 d4=: 21.5 22.8 21 23 21.6 23.6 22.5 20.7 23.4 21.8 20.7 21.7 21.5 22.5 23.6 21.5 22.5 23.5 21.5 21.8 d5=: 19.8 20.4 19.6 17.8 18.5 18.9 18.3 18.9 19.5 22 d6=: 28.2 26.6 20.1 23.3 25.2 22.1 17.7 27.6 20.6 13.7 23.2 17.5 20.6 18 23.9 21.6 24.3 20.4 24 13.2 d7=: 30.02 29.99 30.11 29.97 30.01 29.99 d8=: 29.89 29.93 29.72 29.98 30.02 29.98 d9=: 3 4 1 2.1 da=: 490.2 340 433.9</lang>
Task examples: <lang J> d1 p2_tail d2 0.021378
d3 p2_tail d4
0.148842
d5 p2_tail d6
0.0359723
d7 p2_tail d8
0.0907733
d9 p2_tail da
0.0107377</lang>
Kotlin
This program brings in code from other tasks for gamma functions and integration by Simpson's rule as Kotlin doesn't have these built-in: <lang scala>// version 1.1.4-3
typealias Func = (Double) -> Double
fun square(d: Double) = d * d
fun sampleVar(da: DoubleArray): Double {
if (da.size < 2) throw IllegalArgumentException("Array must have at least 2 elements") val m = da.average() return da.map { square(it - m) }.sum() / (da.size - 1)
}
fun welch(da1: DoubleArray, da2: DoubleArray): Double {
val temp = sampleVar(da1) / da1.size + sampleVar(da2) / da2.size return (da1.average() - da2.average()) / Math.sqrt(temp)
}
fun degreesFreedom(da1: DoubleArray, da2: DoubleArray): Double {
val s1 = sampleVar(da1) val s2 = sampleVar(da2) val n1 = da1.size val n2 = da2.size val temp1 = square(s1 / n1 + s2 / n2) val temp2 = square(s1) / (n1 * n1 * (n1 - 1)) + square(s2) / (n2 * n2 * (n2 - 1)) return temp1 / temp2
}
fun gamma(d: Double): Double {
var dd = d val p = doubleArrayOf( 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ) val g = 7 if (dd < 0.5) return Math.PI / (Math.sin(Math.PI * dd) * gamma(1.0 - dd)) dd-- var a = p[0] val t = dd + g + 0.5 for (i in 1 until p.size) a += p[i] / (dd + i) return Math.sqrt(2.0 * Math.PI) * Math.pow(t, dd + 0.5) * Math.exp(-t) * a
}
fun lGamma(d: Double) = Math.log(gamma(d))
fun simpson(a: Double, b: Double, n: Int, f: Func): Double {
val h = (b - a) / n var sum = 0.0 for (i in 0 until n) { val x = a + i * h sum += (f(x) + 4.0 * f(x + h / 2.0) + f(x + h)) / 6.0 } return sum * h
}
fun p2Tail(da1: DoubleArray, da2: DoubleArray): Double {
val nu = degreesFreedom(da1, da2) val t = welch(da1, da2) val g = Math.exp(lGamma(nu / 2.0) + lGamma(0.5) - lGamma(nu / 2.0 + 0.5)) val b = nu / (t * t + nu) val f: Func = { r -> Math.pow(r, nu / 2.0 - 1.0) / Math.sqrt(1.0 - r) } return simpson(0.0, b, 10000, f) / g // n = 10000 seems more than enough here
}
fun main(args: Array<String>) {
val da1 = doubleArrayOf( 27.5, 21.0, 19.0, 23.6, 17.0, 17.9, 16.9, 20.1, 21.9, 22.6, 23.1, 19.6, 19.0, 21.7, 21.4 ) val da2 = doubleArrayOf( 27.1, 22.0, 20.8, 23.4, 23.4, 23.5, 25.8, 22.0, 24.8, 20.2, 21.9, 22.1, 22.9, 20.5, 24.4 ) val da3 = doubleArrayOf( 17.2, 20.9, 22.6, 18.1, 21.7, 21.4, 23.5, 24.2, 14.7, 21.8 ) val da4 = doubleArrayOf( 21.5, 22.8, 21.0, 23.0, 21.6, 23.6, 22.5, 20.7, 23.4, 21.8, 20.7, 21.7, 21.5, 22.5, 23.6, 21.5, 22.5, 23.5, 21.5, 21.8 ) val da5 = doubleArrayOf( 19.8, 20.4, 19.6, 17.8, 18.5, 18.9, 18.3, 18.9, 19.5, 22.0 ) val da6 = doubleArrayOf( 28.2, 26.6, 20.1, 23.3, 25.2, 22.1, 17.7, 27.6, 20.6, 13.7, 23.2, 17.5, 20.6, 18.0, 23.9, 21.6, 24.3, 20.4, 24.0, 13.2 ) val da7 = doubleArrayOf(30.02, 29.99, 30.11, 29.97, 30.01, 29.99) val da8 = doubleArrayOf(29.89, 29.93, 29.72, 29.98, 30.02, 29.98)
val x = doubleArrayOf(3.0, 4.0, 1.0, 2.1) val y = doubleArrayOf(490.2, 340.0, 433.9) val f = "%.6f" println(f.format(p2Tail(da1, da2))) println(f.format(p2Tail(da3, da4))) println(f.format(p2Tail(da5, da6))) println(f.format(p2Tail(da7, da8))) println(f.format(p2Tail(x, y)))
}</lang>
- Output:
0.021378 0.148842 0.035972 0.090773 0.010751
Perl 6
Perhaps "inspired by C example" may be more accurate. Gamma subroutine from Gamma function task.
<lang perl6>sub Γ(\z) {
constant g = 9; z < .5 ?? π / sin(π * z) / Γ(1 - z) !! τ.sqrt * (z + g - 1/2)**(z - 1/2) * exp(-(z + g - 1/2)) * [+] < 1.000000000000000174663 5716.400188274341379136 -14815.30426768413909044 14291.49277657478554025 -6348.160217641458813289 1301.608286058321874105 -108.1767053514369634679 2.605696505611755827729 -0.7423452510201416151527e-2 0.5384136432509564062961e-7 -0.4023533141268236372067e-8 > Z* 1, |map 1/(z + *), 0..*
}
sub p-value (@A, @B) {
return 1 if @A <= 1 or @B <= 1;
my $a-mean = @A.sum / @A; my $b-mean = @B.sum / @B; my $a-variance = @A.map( { ($a-mean - $_)² } ).sum / (@A - 1); my $b-variance = @B.map( { ($b-mean - $_)² } ).sum / (@B - 1); return 1 unless $a-variance && $b-variance;
my \Welsh-𝒕-statistic = ($a-mean - $b-mean)/($a-variance/@A + $b-variance/@B).sqrt;
my $DoF = ($a-variance / @A + $b-variance / @B)² / (($a-variance² / (@A³ - @A²)) + ($b-variance² / (@B³ - @B²)));
my $sa = $DoF / 2 - 1; my $x = $DoF / (Welsh-𝒕-statistic² + $DoF); my $N = 65355; my $h = $x / $N; my ( $sum1, $sum2 );
for ^$N »*» $h -> $i { $sum1 += (($i + $h / 2) ** $sa) / (1 - ($i + $h / 2)).sqrt; $sum2 += $i ** $sa / (1 - $i).sqrt; }
(($h / 6) * ( $x ** $sa / (1 - $x).sqrt + 4 * $sum1 + 2 * $sum2)) / ( Γ($sa + 1) * 1.77245385090551610 / Γ($sa + 1.5) );
}
- Testing
for (
[<27.5 21.0 19.0 23.6 17.0 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19.0 21.7 21.4>], [<27.1 22.0 20.8 23.4 23.4 23.5 25.8 22.0 24.8 20.2 21.9 22.1 22.9 20.5 24.4>],
[<17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8>], [<21.5 22.8 21.0 23.0 21.6 23.6 22.5 20.7 23.4 21.8 20.7 21.7 21.5 22.5 23.6 21.5 22.5 23.5 21.5 21.8>],
[<19.8 20.4 19.6 17.8 18.5 18.9 18.3 18.9 19.5 22.0>], [<28.2 26.6 20.1 23.3 25.2 22.1 17.7 27.6 20.6 13.7 23.2 17.5 20.6 18.0 23.9 21.6 24.3 20.4 24.0 13.2>],
[<30.02 29.99 30.11 29.97 30.01 29.99>], [<29.89 29.93 29.72 29.98 30.02 29.98>],
[<3.0 4.0 1.0 2.1>], [<490.2 340.0 433.9>]
) -> @left, @right { say p-value @left, @right }</lang>
- Output:
0.0213780014628669 0.148841696605328 0.0359722710297969 0.0907733242856673 0.010751534033393
R
<lang R>#!/usr/bin/R d1 <- c(27.5,21.0,19.0,23.6,17.0,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19.0,21.7,21.4) d2 <- c(27.1,22.0,20.8,23.4,23.4,23.5,25.8,22.0,24.8,20.2,21.9,22.1,22.9,20.5,24.4) d3 <- c(17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8) d4 <- c(21.5,22.8,21.0,23.0,21.6,23.6,22.5,20.7,23.4,21.8,20.7,21.7,21.5,22.5,23.6,21.5,22.5,23.5,21.5,21.8) d5 <- c(19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0) d6 <- c(28.2,26.6,20.1,23.3,25.2,22.1,17.7,27.6,20.6,13.7,23.2,17.5,20.6,18.0,23.9,21.6,24.3,20.4,24.0,13.2) d7 <- c(30.02,29.99,30.11,29.97,30.01,29.99) d8 <- c(29.89,29.93,29.72,29.98,30.02,29.98) x <- c(3.0,4.0,1.0,2.1) y <- c(490.2,340.0,433.9)
results <- t.test(d1,d2, alternative="two.sided", var.equal=FALSE) print(results$p.value) results <- t.test(d3,d4, alternative="two.sided", var.equal=FALSE) print(results$p.value) results <- t.test(d5,d6, alternative="two.sided", var.equal=FALSE) print(results$p.value) results <- t.test(d7,d8, alternative="two.sided", var.equal=FALSE) print(results$p.value) results <- t.test(x,y, alternative="two.sided", var.equal=FALSE) print(results$p.value) </lang>
- Output:
[1] 0.021378 [1] 0.1488417 [1] 0.03597227 [1] 0.09077332 [1] 0.01075156
Racket
, producing the same output.
<lang racket>#lang racket (require math/statistics math/special-functions)
(define (p-value S1 S2 #:n (n 11000))
(define σ²1 (variance S1 #:bias #t)) (define σ²2 (variance S2 #:bias #t)) (define N1 (sequence-length S1)) (define N2 (sequence-length S2)) (define σ²/sz1 (/ σ²1 N1)) (define σ²/sz2 (/ σ²2 N2)) (define degrees-of-freedom (/ (sqr (+ σ²/sz1 σ²/sz2)) (+ (/ (sqr σ²1) (* (sqr N1) (sub1 N1))) (/ (sqr σ²2) (* (sqr N2) (sub1 N2)))))) (define a (/ degrees-of-freedom 2)) (define a-1 (sub1 a)) (define x (let ((welch-t-statistic (/ (- (mean S1) (mean S2)) (sqrt (+ σ²/sz1 σ²/sz2))))) (/ degrees-of-freedom (+ (sqr welch-t-statistic) degrees-of-freedom)))) (define h (/ x n)) (/ (* (/ h 6) (+ (* (expt x a-1) (expt (- 1 x) -1/2)) (* 4 (for/sum ((i (in-range 0 n))) (* (expt (+ (* h i) (/ h 2)) a-1) (expt (- 1 (+ (* h i) (/ h 2))) -1/2)))) (* 2 (for/sum ((i (in-range 0 n))) (* (expt (* h i) a-1) (expt (- 1 (* h i)) -1/2)))))) (* (gamma a) 1.77245385090551610 (/ (gamma (+ a 1/2))))))
(module+ test
(list (p-value (list 27.5 21.0 19.0 23.6 17.0 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19.0 21.7 21.4) (list 27.1 22.0 20.8 23.4 23.4 23.5 25.8 22.0 24.8 20.2 21.9 22.1 22.9 20.5 24.4)) (p-value (list 17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8) (list 21.5 22.8 21.0 23.0 21.6 23.6 22.5 20.7 23.4 21.8 20.7 21.7 21.5 22.5 23.6 21.5 22.5 23.5 21.5 21.8)) (p-value (list 19.8 20.4 19.6 17.8 18.5 18.9 18.3 18.9 19.5 22.0) (list 28.2 26.6 20.1 23.3 25.2 22.1 17.7 27.6 20.6 13.7 23.2 17.5 20.6 18.0 23.9 21.6 24.3 20.4 24.0 13.2)) (p-value (list 30.02 29.99 30.11 29.97 30.01 29.99) (list 29.89 29.93 29.72 29.98 30.02 29.98)) (p-value (list 3.0 4.0 1.0 2.1) (list 490.2 340.0 433.9))))</lang>
- Output:
(0.021378001462867013 0.14884169660532798 0.035972271029796624 0.09077332428567102 0.01075139991904718)
Tcl
This is not particularly idiomatic Tcl, but perhaps illustrates some of the language's relationship with the Lisp family.
<lang Tcl>#!/usr/bin/tclsh
package require math::statistics package require math::special namespace path {::math::statistics ::math::special ::tcl::mathfunc ::tcl::mathop}
proc incf {_var {inc 1.0}} {
upvar 1 $_var var if {![info exists var]} { set var 0.0 } set var [expr {$inc + $var}]
}
proc sumfor {_var A B body} {
upvar 1 $_var var set var $A set res 0 while {$var < $B} { incf res [uplevel 1 $body] incr var } return $res
}
proc sqr {x} {expr {$x*$x}}
proc pValue {S1 S2 {n 11000}} {
set σ²1 [var $S1] set σ²2 [var $S2] set N1 [llength $S1] set N2 [llength $S2] set σ²/sz1 [/ ${σ²1} $N1] set σ²/sz2 [/ ${σ²2} $N2]
set d1 [/ [sqr ${σ²1}] [* [sqr $N1] [- $N1 1]]] set d2 [/ [sqr ${σ²2}] [* [sqr $N2] [- $N2 1]]] set DoF [/ [sqr [+ ${σ²/sz1} ${σ²/sz2}]] [+ $d1 $d2]]
set a [/ $DoF 2.0]
set welchTstat [/ [- [mean $S1] [mean $S2]] [sqrt [+ ${σ²/sz1} ${σ²/sz2}]]] set x [/ $DoF [+ [sqr $welchTstat] $DoF]] set h [/ $x $n]
/ [* [/ $h 6] \ [+ [* [** $x [- $a 1]] \ [** [- 1 $x] -0.5]] \ [* 4 [sumfor i 0 $n { * [** [+ [* $h $i] [/ $h 2]] [- $a 1]] \ [** [- 1 [* $h $i] [/ $h 2]] -0.5]}]] \ [* 2 [sumfor i 0 $n { * [** [* $h $i] [- $a 1]] [** [- 1 [* $h $i]] -0.5]}]]]] \ [* [Gamma $a] 1.77245385090551610 [/ 1.0 [Gamma [+ $a 0.5]]]]
}
foreach {left right} {
{ 27.5 21.0 19.0 23.6 17.0 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19.0 21.7 21.4 } { 27.1 22.0 20.8 23.4 23.4 23.5 25.8 22.0 24.8 20.2 21.9 22.1 22.9 20.5 24.4 }
{ 17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8 } { 21.5 22.8 21.0 23.0 21.6 23.6 22.5 20.7 23.4 21.8 20.7 21.7 21.5 22.5 23.6 21.5 22.5 23.5 21.5 21.8 } { 19.8 20.4 19.6 17.8 18.5 18.9 18.3 18.9 19.5 22.0 } { 28.2 26.6 20.1 23.3 25.2 22.1 17.7 27.6 20.6 13.7 23.2 17.5 20.6 18.0 23.9 21.6 24.3 20.4 24.0 13.2 } { 30.02 29.99 30.11 29.97 30.01 29.99 } { 29.89 29.93 29.72 29.98 30.02 29.98 } { 3.0 4.0 1.0 2.1 } { 490.2 340.0 433.9 }
} {
puts [pValue $left $right]
} </lang>
- Output:
0.021378001462853034 0.148841696604164 0.035972271029770915 0.09077332428458083 0.010751399918798182
zkl
<lang zkl>fcn calculate_Pvalue(array1,array2){
if (array1.len()<=1 or array2.len()<=1) return(1.0);
mean1,mean2 := array1.sum(0.0),array2.sum(0.0); if(mean1==mean2) return(1.0); mean1/=array1.len(); mean2/=array2.len();
variance1:=array1.reduce('wrap(sum,x){ sum + (x-mean1).pow(2) },0.0); variance2:=array2.reduce('wrap(sum,x){ sum + (x-mean2).pow(2) },0.0);
variance1/=(array1.len() - 1); variance2/=(array2.len() - 1);
WELCH_T_STATISTIC:=(mean1-mean2)/ (variance1/array1.len() + variance2/array2.len()).sqrt(); DEGREES_OF_FREEDOM:= ( variance1/array1.len() + variance2/array2.len() ).pow(2) // numerator / ( (variance1*variance1)/(array1.len().pow(2)*(array1.len() - 1)) +
(variance2*variance2)/(array2.len().pow(2)*(array2.len() - 1))
); a:=DEGREES_OF_FREEDOM/2; x:=DEGREES_OF_FREEDOM/( WELCH_T_STATISTIC.pow(2) + DEGREES_OF_FREEDOM ); N,h := 65535, x/N;
sum1,sum2 := 0.0, 0.0; foreach i in (N){ sum1+=((h*i + h/2.0).pow(a - 1))/(1.0 - (h*i + h/2.0)).sqrt(); sum2+=((h*i).pow(a - 1))/(1.0 - h*i).sqrt(); } return_value:=((h/6.0)*( x.pow(a - 1)/(1.0 - x).sqrt() + 4.0*sum1 + 2.0*sum2) ) / ((0.0).e.pow(lngammal(a) + 0.57236494292470009 - lngammal(a + 0.5)));
if(return_value > 1.0) return(1.0); // or return_value is infinite, throws return_value;
} fcn lngammal(xx){
var [const] cof=List( // static 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5 ); y:=x:=xx; tmp:=x + 5.5 - (x + 0.5) * (x + 5.5).log(); ser:=1.000000000190015; foreach x in (cof){ ser+=(x/(y+=1)); } return((2.5066282746310005 * ser / x).log() - tmp);
}</lang> <lang zkl>testSets:=T( T(T(27.5,21.0,19.0,23.6,17.0,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19.0,21.7,21.4),
T(27.1,22.0,20.8,23.4,23.4,23.5,25.8,22.0,24.8,20.2,21.9,22.1,22.9,20.5,24.4)),
T(T(17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8),
T(21.5,22.8,21.0,23.0,21.6,23.6,22.5,20.7,23.4,21.8,20.7,21.7,21.5,22.5,23.6,21.5,22.5,23.5,21.5,21.8)),
T(T(19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0),
T(28.2,26.6,20.1,23.3,25.2,22.1,17.7,27.6,20.6,13.7,23.2,17.5,20.6,18.0,23.9,21.6,24.3,20.4,24.0,13.2)),
T(T(30.02,29.99,30.11,29.97,30.01,29.99),
T(29.89,29.93,29.72,29.98,30.02,29.98)),
T(T(3.0,4.0,1.0,2.1),T(490.2,340.0,433.9)) );
foreach x,y in (testSets)
{ println("Test set 1 p-value = %f".fmt(calculate_Pvalue(x,y))); }</lang>
- Output:
Test set 1 p-value = 0.021378 Test set 1 p-value = 0.148842 Test set 1 p-value = 0.035972 Test set 1 p-value = 0.090773 Test set 1 p-value = 0.010752