1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
| | %prep
if ! zmodload zsh/mathfunc 2>/dev/null; then
ZTST_unimplemented="can't load the zsh/mathfunc module for testing"
fi
print -ru $ZTST_fd 'This test may take two seconds...'
function calc_chi2() {
# Calculate Chi Squared
integer n=$1
for ((k=1;k<=n;k++))
do
((delta=${observed[$k]}-expected))
#echo "$k=$delta \c"
((chi2+=delta**2/expected))
done
#echo
}
function ligf() {
float -E S=$1 Z=$2
float -gE RET
if (( Z < 0.0 )); then
RET=0.0
fi
float -E Sc=$((1.0/S))
((Sc *= Z**S ))
((Sc *= exp(-1*Z) ))
float -E Sum=1.0 Nom=1.0 Denom=1.0
integer I
# 200 iterations seems like enough for an approximation of an infinate
# series
for (( I=0;I<200;I++))
do
((Nom *= Z))
((S++))
((Denom *= S))
((Sum += (Nom/Denom) ))
done
RET=$((Sum * Sc))
}
#echo $observed
function calc_p() {
#echo "chi2=$chi2"
float chi2=0 cdf=0
float RET
calc_chi2 "$n"
ligf $(((n-1)*0.5)) $((chi2*0.5))
((cdf=$RET/gamma((n-1)*0.5) ))
((p=1-cdf))
}
%test
fail=0
for ((o=0;o<100;o++))
do
n=20
samples=1000
integer j=0
float -F expected=$((samples*1.0/n))
float -F p
typeset -a observed=([$n]=0)
for ((i=1;i<=samples;i++))
do
j=$(( zrand_int(n,1,1) ))
((observed[$j]++))
done
calc_p
if ((p<0.05)); then
((fail++))
fi
done
((fail <= 10))
0q:random integer $samples samples between 1-$n
fail=0
unset observed
for ((o=0;o<100;o++))
do
n=20
samples=1000
integer j=0
float -F expected=$((samples*1.0/n))
float -F p
typeset -a observed=([$n]=0)
for ((i=1;i<=samples;i++))
do
j=$(( ceil(zrand_float()*n) ))
((observed[$j]++))
done
calc_p
if ((p<0.05)); then
((fail++))
fi
done
((fail <= 10))
0q:random float $samples samples
fail=0
unset observed
for ((o=0;o<100;o++))
do
n=16
samples=1000
integer j=0
float -F expected=$((samples*1.0/n))
float -F p
typeset -a observed=([$n]=0)
for ((i=1;i<=samples;i++))
do
j=$(( (SRANDOM & 0x0F)+1 ))
((observed[$j]++))
done
calc_p
if ((p<0.05)); then
((fail++))
fi
done
((fail <= 10))
0q:SRANDOM mod 16 $samples samples
|