Examples - Statistics

Standard deviation and Variance

x=(25, 24, 23, 22, 20, 18, 17, 16, 14, 10, 8);
print "mean=", ave x, "
(x-mean)=", listfor(i,0,count x-1, x[i]-Ans), "
Sum(x-mean)^2=", sumq Ans, "
Variance=", var x, "
Standard Deviation=", stdev x

Standard deviation with frequency list

Data=(14,13,12,11,10,9,8,7,6,5);
Freq=(1,2,3,3,5,3,2,1,2,1);
print "n=",n=sum Freq, "
sum x=", s=sum(Data*Freq), "
sum x^2=", s2=sumfor(i,0,count Data-1,Data[i]^2 * Freq[i]), "
mean=", s/n, "
variance=", (s2-s^2/n)/(n-1), "
standard deviation=", sqrt Ans

Two lists

X=(4,7,3); Y=(5,1,6);
ave X \   /* (4+7+3)/3 */
ave Y \   /* (5+1+6)/3 */
X*Y \     /* 4*5 + 7*1 + 3*6 */
sumq(X-Y) /* (4-5)^2 + (7-1)^2 + (3-6)^2 */

Two lists in a matrix

B=(4,5\7,1\3,6);
avex B \    /* (4+7+3)/3 */
avey B \    /* (5+1+6)/3 */
sumxy B \   /* 4*5 + 7*1 + 3*6 */
sumq(B[][0]-B[][1]) /* (4-5)^2 + (7-1)^2 + (3-6)^2 */

Student's t-test Wikipedia

x=(10, 9, 9, 8, 7, 7, 7, 6, 6, 5, 3, 2);
y=(14, 13, 13, 13, 12, 12, 10, 8, 8, 7, 7, 5, 5);
n1=count x; n2=count y;
print"Student t =",(ave(x)-ave(y))/
  sqrt(((n1-1)*var(x)+(n2-1)*var(y))/(n1+n2-2)*(1/n1+1/n2))

Student's t-test with frequency lists

Data1=(10, 9, 8, 7, 6, 5, 3, 2);
Freq1=(1, 2, 1, 3, 2, 1, 1, 1);
Data2=(14, 13, 12, 10, 8, 7, 5);
Freq2=(1, 3, 2, 1, 2, 2, 2);
print "n1=",n1=sum Freq1, "
mean1=",m1=sum(Data1*Freq1)/n1, "
variance1=", v1=(sumfor(i,0,count Data1-1,Data1[i]^2 * Freq1[i]) -
 m1^2*n1)/(n1-1), "
n2=",n2=sum Freq2, "
mean2=",m2=sum(Data2*Freq2)/n2, "
variance2=", v2=(sumfor(i,0,count Data2-1,Data2[i]^2 * Freq2[i]) -
 m2^2*n2)/(n2-1), "
t =",(m1-m2)/sqrt(((n1-1)*v1+(n2-1)*v2)/(n1+n2-2)*(1/n1+1/n2))

Pearson r Wikipedia

x=(100,40,95,90,92,85,55,60,98,20);
y=(0,95,5,20,30,40,50,70,0,100);
/* (n*(x*y)-sum(x)*sum(y))/sqrt((n*sumq(x)-sum(x)^2)*(n*sumq(y)-sum(y)^2) */
print "Pearson r=", lrr(x',y')

Spearman's rho Wikipedia

X=(86,97,99,100,101,103,106,110,112,113);
Y=(0,20,28,27,50,29,7,17,6,12);
k=0;
rank:
L=if(k==0,X,Y);
i=0;
loop:
b=0; e=0; m=L[i];
foreach(j,L, if(j>m, 0, if(j==m, e++, b++)));
j=b+(1+e)/2;
R= if(i==0, j, (R,j));
i++; if(i<count L, goto loop, 0);
if(k==0, Rx=R, Ry=R);
k++; if(k<2, goto rank, 0);
print "Spearman's rho=", lrr(Rx',Ry')

Rank sum test (Wilcoxon & Mann-Whitney) Wikipedia

L1=(12, 12, 11, 8, 7, 6, 6);
L2=(24, 24, 16, 15, 14, 12, 11, 10, 9);
S1=sort(L1);S2=sort(L2); n1=count L1;n2=count L2; 
R1=0;R2=0;i=0;j=0;o=1;t1=0;t2=0;last=0;ties=0;
merge:
if(i<n1, goto merge1, 0);
if(j<n2, goto take2, goto addrank);
merge1: if(j>=n2, goto take1, 0);
if(S1[i]<S2[j], goto take1, 0);
take2: b=0; item=S2[j]; j++; goto compare;
take1: b=1; item=S1[i]; i++;
compare: if(item==last, goto next, 0);
addrank:
a=t1+t2; if(a>0,ties=ties+a-1,0);
a=o-(a+1)/2; R1=R1+a*t1; R2=R2+a*t2;
if(o>n1+n2, goto end, 0);
t1=0;t2=0; last=item;
next: o++; if(b,t1++,t2++); goto merge;
end:
U1=R1-n1*(n1+1)/2; U2=R2-n2*(n2+1)/2;
if(U1+U2==n1*n2, gotor3,0); print "error"; return;
print "SUM OF RANKS 1=",R1, "
SUM OF RANKS 2=",R2, "
NUMBER OF TIES=",ties, "
U=",U=max(U1,U2), "
MEAN U=",Mu=n1*n2/2, "
STD DEV U=",Su=sqrt(n1*n2*(n1+n2+1)/12), "
Z=",(U-Mu)/Su

Wilcoxon match-pairs signed ranks test Wikipedia

L1=(51,49,46,45,46,39,38,40,41,40,49,55,45,37,44,52,37);
L2=(50,48,46,44,45,41,39,39,39,42,48,56,48,36,45,51,36);
D=sort(L1-L2); n=count D;
i=-1;
i++; if(i==n, gotor1, if(D[i]<0, gotor-1, 0));
if(i>0,gotor3,0);
print "Sum of ""-"" ranks= 0"; return;
S1=listforeach(x, reverse D[0,i-1], abs x);
i--;
i++; if(i==n, gotor1, if(D[i]==0, gotor-1, 0));
if(i<n,gotor3,0);
print "Sum of ""+"" ranks= 0"; return;
S2=D[i,n-1];
n1=count S1; n2=count S2;
R1=0;R2=0;i=0;j=0;o=1;t1=0;t2=0;last=0;
merge:
if(i<n1, goto merge1, 0);
if(j<n2, goto take2, goto addrank);
merge1: if(j>=n2, goto take1, 0);
if(S1[i]<S2[j], goto take1, 0);
take2: b=0; item=S2[j]; j++; goto compare;
take1: b=1; item=S1[i]; i++;
compare: if(item==last, goto next, 0);
addrank:
a=t1+t2;
a=o-(a+1)/2; R1=R1+a*t1; R2=R2+a*t2;
if(o>n1+n2, goto end, 0);
t1=0;t2=0; last=item;
next: o++; if(b,t1++,t2++); goto merge;
end:
print "Sum of ""-"" ranks=",R1, "
Sum of ""+"" ranks=",R2,"
N=",N=n1+n2,"
W=",abs(R1-R2),"
Normal distribution approximation:
Z=",z=(abs(R1-R2)-.5)/sqrt(N*(N+1)*(2*N+1)/6),"
p=",2*exp(-z^2/2)/sqrt(2*pi) * polynom(1/(1+0.2316419*z),
 (0,0.319381530,-0.356563782,1.781477937,-1.821255978,1.330274429))