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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%A efficiency.tex DESIGN documentation Leonard Soicher
%
%
%
\def\DESIGN{\sf DESIGN}
\def\GRAPE{\sf GRAPE}
\def\nauty{\it nauty}
\def\Aut{{\rm Aut}\,}
\def\x{\times}
\Chapter{Matrices and efficiency measures for block designs}
In this chapter we describe functions to calculate certain matrices
associated with a block design, and the function `BlockDesignEfficiency'
which determines certain statistical efficiency measures of a 1-design.
We also document the utility function `DESIGN_IntervalForLeastRealZero',
which is used in the calculation of E-efficiency measures, but has much
wider application.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Matrices associated with a block design}
\>PointBlockIncidenceMatrix( <D> )
This function returns the point-block incidence matrix <N> of the
block design <D>. This matrix has rows indexed by the points of <D>
and columns by the blocks of <D>, with the $(i,j)$-entry of <N> being
the number of times point $i$ occurs in `<D>.blocks[<j>]'.
The returned matrix <N> is immutable.
\beginexample
gap> D:=DualBlockDesign(AGPointFlatBlockDesign(2,3,1));;
gap> BlockDesignBlocks(D);
[ [ 1, 2, 3, 4 ], [ 1, 5, 6, 7 ], [ 1, 8, 9, 10 ], [ 2, 5, 8, 11 ],
[ 2, 7, 9, 12 ], [ 3, 5, 10, 12 ], [ 3, 6, 9, 11 ], [ 4, 6, 8, 12 ],
[ 4, 7, 10, 11 ] ]
gap> PointBlockIncidenceMatrix(D);
[ [ 1, 1, 1, 0, 0, 0, 0, 0, 0 ], [ 1, 0, 0, 1, 1, 0, 0, 0, 0 ],
[ 1, 0, 0, 0, 0, 1, 1, 0, 0 ], [ 1, 0, 0, 0, 0, 0, 0, 1, 1 ],
[ 0, 1, 0, 1, 0, 1, 0, 0, 0 ], [ 0, 1, 0, 0, 0, 0, 1, 1, 0 ],
[ 0, 1, 0, 0, 1, 0, 0, 0, 1 ], [ 0, 0, 1, 1, 0, 0, 0, 1, 0 ],
[ 0, 0, 1, 0, 1, 0, 1, 0, 0 ], [ 0, 0, 1, 0, 0, 1, 0, 0, 1 ],
[ 0, 0, 0, 1, 0, 0, 1, 0, 1 ], [ 0, 0, 0, 0, 1, 1, 0, 1, 0 ] ]
\endexample
\>ConcurrenceMatrix( <D> )
This function returns the concurrence matrix <L> of the block design <D>.
This matrix is equal to $NN^{\rm T}$, where $N$ is the point-block
incidence matrix of <D> (see "PointBlockIncidenceMatrix") and
$N^{\rm T}$ is the transpose of $N$. If <D> is a binary block design
then the $(i,j)$-entry of its concurrence matrix is the number of blocks
containing points $i$ and $j$.
The returned matrix <L> is immutable.
\beginexample
gap> D:=DualBlockDesign(AGPointFlatBlockDesign(2,3,1));;
gap> BlockDesignBlocks(D);
[ [ 1, 2, 3, 4 ], [ 1, 5, 6, 7 ], [ 1, 8, 9, 10 ], [ 2, 5, 8, 11 ],
[ 2, 7, 9, 12 ], [ 3, 5, 10, 12 ], [ 3, 6, 9, 11 ], [ 4, 6, 8, 12 ],
[ 4, 7, 10, 11 ] ]
gap> ConcurrenceMatrix(D);
[ [ 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0 ],
[ 1, 3, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1 ],
[ 1, 1, 3, 1, 1, 1, 0, 0, 1, 1, 1, 1 ],
[ 1, 1, 1, 3, 0, 1, 1, 1, 0, 1, 1, 1 ],
[ 1, 1, 1, 0, 3, 1, 1, 1, 0, 1, 1, 1 ],
[ 1, 0, 1, 1, 1, 3, 1, 1, 1, 0, 1, 1 ],
[ 1, 1, 0, 1, 1, 1, 3, 0, 1, 1, 1, 1 ],
[ 1, 1, 0, 1, 1, 1, 0, 3, 1, 1, 1, 1 ],
[ 1, 1, 1, 0, 0, 1, 1, 1, 3, 1, 1, 1 ],
[ 1, 0, 1, 1, 1, 0, 1, 1, 1, 3, 1, 1 ],
[ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 0 ],
[ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 3 ] ]
\endexample
\>InformationMatrix( <D> )
This function returns the information matrix <C> of the block design <D>.
This matrix is defined as follows. Suppose <D> has $v$ points and $b$
blocks, let $R$ be the $v\times v$ diagonal matrix whose $(i,i)$-entry
is the replication number of the point $i$, let $N$ be the point-block
incidence matrix of <D> (see "PointBlockIncidenceMatrix"), and let $K$
be the $b\times b$ diagonal matrix whose $(j,j)$-entry is the length of
`<D>.blocks[<j>]'. Then the *information matrix* of <D> is
$C:=R-NK^{-1}N^{\rm T}$. If <D> is a 1-$(v,k,r)$ design then this expression
for <C> simplifies to $rI-k^{-1}L$, where $I$ is the $v\times v$ identity
matrix and $L$ is the concurrence matrix of <D> (see "ConcurrenceMatrix").
The returned matrix <C> is immutable.
\beginexample
gap> D:=DualBlockDesign(AGPointFlatBlockDesign(2,3,1));;
gap> BlockDesignBlocks(D);
[ [ 1, 2, 3, 4 ], [ 1, 5, 6, 7 ], [ 1, 8, 9, 10 ], [ 2, 5, 8, 11 ],
[ 2, 7, 9, 12 ], [ 3, 5, 10, 12 ], [ 3, 6, 9, 11 ], [ 4, 6, 8, 12 ],
[ 4, 7, 10, 11 ] ]
gap> InformationMatrix(D);
[ [ 9/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, 0, 0 ],
[ -1/4, 9/4, -1/4, -1/4, -1/4, 0, -1/4, -1/4, -1/4, 0, -1/4, -1/4 ],
[ -1/4, -1/4, 9/4, -1/4, -1/4, -1/4, 0, 0, -1/4, -1/4, -1/4, -1/4 ],
[ -1/4, -1/4, -1/4, 9/4, 0, -1/4, -1/4, -1/4, 0, -1/4, -1/4, -1/4 ],
[ -1/4, -1/4, -1/4, 0, 9/4, -1/4, -1/4, -1/4, 0, -1/4, -1/4, -1/4 ],
[ -1/4, 0, -1/4, -1/4, -1/4, 9/4, -1/4, -1/4, -1/4, 0, -1/4, -1/4 ],
[ -1/4, -1/4, 0, -1/4, -1/4, -1/4, 9/4, 0, -1/4, -1/4, -1/4, -1/4 ],
[ -1/4, -1/4, 0, -1/4, -1/4, -1/4, 0, 9/4, -1/4, -1/4, -1/4, -1/4 ],
[ -1/4, -1/4, -1/4, 0, 0, -1/4, -1/4, -1/4, 9/4, -1/4, -1/4, -1/4 ],
[ -1/4, 0, -1/4, -1/4, -1/4, 0, -1/4, -1/4, -1/4, 9/4, -1/4, -1/4 ],
[ 0, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, 9/4, 0 ],
[ 0, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, -1/4, 0, 9/4 ] ]
\endexample
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{The function BlockDesignEfficiency}
\>BlockDesignEfficiency( <D> )
\>BlockDesignEfficiency( <D>, <eps> )
\>BlockDesignEfficiency( <D>, <eps>, <includeMV> )
Let <D> be a 1-$(v,k,r)$ design with $v>1$, let <eps> be a positive
rational number (default: $10^{-6}$), and let <includeMV> be a boolean
(default: `false'). Then this function returns a record <eff> containing
information on statistical efficiency measures of <D>. These measures
are defined below. See \cite{Extrep}, \cite{BaCa} and \cite{BaRo}
for further details. All returned results are computed using exact
algebraic computation.
The component `<eff>.A' contains the A-efficiency measure for <D>,
`<eff>.Dpowered' contains the D-efficiency measure of <D> raised to the
power $v-1$, and `<eff>.Einterval' is a list $[a,b]$ of non-negative
rational numbers such that if $x$ is the E-efficiency measure of <D>
then $a\le x\le b$, $b-a\le$<eps>, and if $x$ is rational then $a=x=b$.
Moreover `<eff>.CEFpolynomial' contains the monic polynomial over the
rationals whose zeros (counting multiplicities) are the canonical
efficiency factors of the design <D>. If `<includeMV>=true' then
additional work is done to compute the MV- (also called E$'$-) efficiency
measure, and then `<eff>.MV' contains the value of this measure. (This
component may be set even if `<includeMV>=false', as a byproduct of
other computation.)
We now define the canonical efficiency factors and the A-, D-, E-,
and MV-efficiency measures of a 1-design.
Let $D$ be a 1-$(v,k,r)$ design with $v\ge 2$, let $C$ be the information
matrix of $D$ (see "InformationMatrix"), and let $F:=r^{-1}C$.
The eigenvalues of $F$ are all real and lie in the interval $[0,1]$.
At least one of these eigenvalues is zero: an associated eigenvector is
the all-$1$ vector. The remaining eigenvalues $\delta_1\le \delta_2\le
\cdots \le \delta_{v-1}$ of $F$ are called the *canonical efficiency
factors* of $D$. These are all non-zero if and only if $D$ is connected
(that is, the point-block incidence graph of $D$ is a connected graph).
If $D$ is not connected, then the A-, D-, E-, and MV-efficiency measures
of $D$ are all defined to be zero. Otherwise, the *A-efficiency
measure* is $(v-1)/\sum_{i=1}^{v-1}1/\delta_i$ (the harmonic mean
of the canonical efficiency factors), the *D-efficiency measure*
is $(\prod_{i=1}^{v-1}\delta_i)^{1/(v-1)}$ (the geometric mean of
the canonical efficiency factors), and the *E-efficiency measure* is
$\delta_1$ (the minimum of the canonical efficiency factors).
If $D$ is connected, and the MV-efficiency measure is required,
then it is computed as follows. Let $F:=r^{-1}C$ be as before,
and let $P:=v^{-1}J$, where $J$ is the $v\times v$ all-1 matrix. Set
$M:=(F+P)^{-1}-P$, making $M$ the ``Moore-Penrose inverse'' of $F$ (see
\cite{BaCa}). Then the *MV-efficiency measure* of $D$ is the minimum
value (over all $i,j\in \{1,\ldots,v\}$, $i\not=j$) of
$2/(M_{ii}+M_{jj}-M_{ij}-M_{ji})$.
\beginexample
gap> D:=DualBlockDesign(AGPointFlatBlockDesign(2,3,1));;
gap> BlockDesignBlocks(D);
[ [ 1, 2, 3, 4 ], [ 1, 5, 6, 7 ], [ 1, 8, 9, 10 ], [ 2, 5, 8, 11 ],
[ 2, 7, 9, 12 ], [ 3, 5, 10, 12 ], [ 3, 6, 9, 11 ], [ 4, 6, 8, 12 ],
[ 4, 7, 10, 11 ] ]
gap> BlockDesignEfficiency(D);
rec( A := 33/41,
CEFpolynomial := x_1^11-9*x_1^10+147/4*x_1^9-719/8*x_1^8+18723/128*x_1^7-106\
47/64*x_1^6+138159/1024*x_1^5-159813/2048*x_1^4+2067201/65536*x_1^3-556227/655\
36*x_1^2+89667/65536*x_1-6561/65536, Dpowered := 6561/65536,
Einterval := [ 3/4, 3/4 ] )
gap> BlockDesignEfficiency(D,10^(-4),true);
rec( A := 33/41,
CEFpolynomial := x_1^11-9*x_1^10+147/4*x_1^9-719/8*x_1^8+18723/128*x_1^7-106\
47/64*x_1^6+138159/1024*x_1^5-159813/2048*x_1^4+2067201/65536*x_1^3-556227/655\
36*x_1^2+89667/65536*x_1-6561/65536, Dpowered := 6561/65536,
Einterval := [ 3/4, 3/4 ], MV := 3/4 )
\endexample
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Section{Computing an interval for a certain real zero of a rational polynomial}
We document a {\DESIGN} package utility function used in the calculation
of the `Einterval' component above, but is more widely applicable.
\>DESIGN_IntervalForLeastRealZero( <f>, <a>, <b>, <eps> )
Suppose that <f> is a univariate polynomial over the rationals, <a>,
<b> are rational numbers with $<a>\le <b>$, and <eps> is a positive
rational number.
If <f> has no real zero in the closed interval $[<a>,<b>]$, then this
function returns the empty list. Otherwise, let $\alpha$ be the least
real zero of <f> such that $<a>\le \alpha\le <b>$. Then this function
returns a list $[c,d]$ of rational numbers, with $c\le \alpha\le d$
and $d-c\le <eps>$. Moreover, $c=d$ if and only if $\alpha$ is rational
(in which case $\alpha=c=d$).
\beginexample
gap> x:=Indeterminate(Rationals,1);
x_1
gap> f:=(x+3)*(x^2-3);
x_1^3+3*x_1^2-3*x_1-9
gap> L:=DESIGN_IntervalForLeastRealZero(f,-5,5,10^(-3));
[ -3, -3 ]
gap> L:=DESIGN_IntervalForLeastRealZero(f,-2,5,10^(-3));
[ -14193/8192, -7093/4096 ]
gap> List(L,Float);
[ -1.73254, -1.73169 ]
gap> L:=DESIGN_IntervalForLeastRealZero(f,0,5,10^(-3));
[ 14185/8192, 7095/4096 ]
gap> List(L,Float);
[ 1.73157, 1.73218 ]
gap> L:=DESIGN_IntervalForLeastRealZero(f,0,5,10^(-5));
[ 454045/262144, 908095/524288 ]
gap> List(L,Float);
[ 1.73204, 1.73205 ]
gap> L:=DESIGN_IntervalForLeastRealZero(f,2,5,10^(-5));
[ ]
\endexample
|
