2021-04-22

直交表の数理～有限射影幾何を用いた直交表の構築編～

前回の記事で有限射影空間の理論を用いた $G$ 行列の構築について説明した。本稿では前回説明した手法で正しく直交表が構築できることを証明し、実際に直交表を構築するプログラム例を示したいと思う。

有限射影幾何を用いた直交表構築手法の証明

おさらい

$k$ を因子数、 $t$ を直交表の強さ、 $q$ を素数冪とする。 $k$ 行 $d$ 列の行列 $G$ があって、その任意の $t$ ( $t \le d$ ) 行は $\mathbb{F}_q$ 上の $d$ 次元ベクトル空間上で一次独立であるとする。このとき、 $d$ 行 $q^d$ 列の行列 $\Theta$ を用いて、直交表は $G\Theta$ という行列積により得られるとのことだった（詳細は前回の記事を参照）。以下ではこの事実の証明を考えていく。

事前に断らせて頂くと、証明はほぼ全て私が考えたものだが、有限体上のベクトル空間の扱いにあまり自信がない。有限体上のベクトル空間では、例えば内積が0でもベクトルが直交しないようなケースがあるなど、標数0の体で使える議論が使えないケースがある[1]。もしそういった点で誤りがあれば、Twitterやブログコメントにてご指摘頂けるとありがたい。

証明

${\bf x}_1, \cdots , {\bf x}_k$ を $\mathbb{F}_q$ 上の $d$ 次元ベクトル空間 $V$ 上の行ベクトルの集合とし、これらのうち任意の $t$ 個を選べばそれらは一次独立になっているとする。これを用いて $G$ は以下のように書ける。

$\displaystyle{ G = \left( \begin{array}{cc} {\bf x}_1 \\ \vdots \\ {\bf x}_k \end{array} \right) }$

同様に ${\bf y}_1, \cdots , {\bf y}_{q^d}$ を $V$ 上の列ベクトルの集合とし、これらは相異なるとする。ベクトルの数が $q^d$ であるため、これらのベクトルは $V$ のベクトル全体の集合となる。これを用いて $\Theta$ は以下のように書ける。

$\displaystyle{ \Theta = \left( \begin{array}{ccc} {\bf y}_1 & \cdots & {\bf y}_{q^d} \end{array} \right) }$

これより、行列積 $G \Theta$ は以下のようになる。

$\displaystyle{ G \Theta = \left( \begin{array}{ccc} {\bf x}_1 {\bf y}_1 & \cdots & {\bf x}_1 {\bf y}_{q^d} \\ \vdots & \ddots & \vdots \\ {\bf x}_k {\bf y}_1 & \cdots & {\bf x}_k {\bf y}_{q^d} \end{array} \right) }$

示したいことは、 $G \Theta$ の任意の $t$ 行を選んだとき、それらの各列の値から順序集合を形成したとき、同じ順序集合がちょうど $q^{d-t}$ 回ずつ現れることである。というのも、列は全部で $q^d$ 個あり、可能な順序集合のパターンは全部で $q^t$ 通りとなるため、 $t$ 個の値の組がそれぞれ $q^d/q^t = q^{d-t}$ 回ずつ出現すれば直交表としての性質を満たすためである。

以下、選ばれた $t$ 個の行の番号を $i_1, \cdots , i_t$ とする。

まず、全ての順序集合が $q^{d-t}$ 回未満しか出現しないと仮定する。このとき、可能な順序集合の組み合わせが $q^t$ 通りしかないため、最大でも $q^t (q^{d-t}-1) = q^d - q^t$ 個の列しか存在できなくなってしまう。そのため仮定は誤りで、少なくとも1つの順序集合は $q^{d-t}$ 回以上出現する。

一番出現回数が多い順序集合について、それが出現する列がちょうど $q^{d-t}$ 個であることを示す。そのために、このような列が $q^{d-t}+1$ 個あると仮定して矛盾を導く。今、 $l_1, \cdots , l_{q^{d-t}+1}$ という番号の列について、 $i_1, \cdots , i_t$ 行目の値から成る順序集合が一致したと仮定する。この時、 ${\bf x}_{i_1}$ について以下の式が成り立つ。

$\displaystyle{ \begin{eqnarray} {\bf x}_{i_1} {\bf y}_{l_\alpha} &=& {\bf x}_{i_1} {\bf y}_{l_\beta} \\ {\bf x}_{i_1} ({\bf y}_{l_\alpha} - {\bf y}_{l_\beta}) &=& 0 \end{eqnarray} }$

ただし、 $\alpha , \beta = 1, \cdots , q^{d-t}+1$ である。

${\bf x}_{i_2}, \cdots , {\bf x}_{i_t}$ についても同様に以下の式が得られる。

$\displaystyle{ \begin{eqnarray} {\bf x}_{i_2} ({\bf y}_{l_\alpha} - {\bf y}_{l_\beta}) &=& 0 \\ &\vdots& \\ {\bf x}_{i_t} ({\bf y}_{l_\alpha} - {\bf y}_{l_\beta}) &=& 0 \end{eqnarray} }$

この関係を線形写像の観点から整理してみる。そのために以下のような行列 $A$ を考える。

$\displaystyle{ A = \left( \begin{array}{c} {\bf x}_{i_1} \\ \vdots \\ {\bf x}_{i_t} \end{array} \right) }$

${\bf x}_{i_1}, \cdots , {\bf x}_{i_t}$ は一次独立なので、 $\mathrm{rank}A = t$ である。そのため、次元定理より $\mathrm{ker}A = d-t$ となる。 ${\bf y}_{l_\alpha} - {\bf y}_{l_\beta} \in \mathrm{ker}A$ なので、 ${\bf y}_{l_1}, \cdots , {\bf y}_{l_{q^{d-t}+1}}$ は $\mathbb{F}_q^d$ において全て同一の $d-t$ 次元空間に含まれる。そこで、 $\mathrm{ker}A$ のベクトル空間としての基底を ${\bf z}_1, \cdots , {\bf z}_{d-t}$ とすると、 ${\bf y}_{l_1}, \cdots , {\bf y}_{l_{q^{d-t}+1}}$ は全て以下の形式で表すことができる。

${\bf y}_{l_1} = c_1 {\bf z}_1 + \cdots + c_{d-t} {\bf z}_{d-t}$

ただし、 $c_1, \cdots , c_{d-t} \in \mathbb{F}_q$ である。これらの係数は $0, \cdots, q-1$ の $q$ 通りの値を取るため、この形式で表されるベクトルは全部で $q^{d-t}$ 本となる。これは ${\bf y}_{l_1}, \cdots , {\bf y}_{l_{q^{d-t}+1}}$ が相異なる $q^{d-t}+1$ 個のベクトルであることに矛盾する。

よって仮定は誤りで、一番出現回数が多い列は $q^{d-t}$ 個または $q^{d-t}+2$ 個以上ということになる。後者は $q^{d-t}+1$ 個の場合と同様の理由でありえないため、結局列数は $q^{d-t}$ 個ということになる。

今考えた一番出現回数が多い順序集合が出現する列を除いて、他の列について同様の議論を繰り返すことで、結局全ての順序集合の出現回数は $q^{d-t}$ 回と言える。よって $G \Theta$ は直交表になっている。

計算してみよう

証明ができて気分がすっきりしたところで、実際に計算してみよう。前回の記事で作成したgmat.pyをimportして直交表を構築するプログラム例を以下に示す。

#!/usr/bin/env sage

import gmat
import sys
from sage.all import *

def calcThetaMatrix(dim, order):
    K = GF(order, name="a", modulus="primitive")
    gamma = K.gen()
    Theta = calcThetaMatrixHelper(dim, order, gamma, 0, [0 for d in range(dim)])
    Theta = Matrix(Theta)
    return Theta.transpose()

def calcThetaMatrixHelper(dim, order, gamma, pos, currentVec):
    if pos == dim:
        return [currentVec.copy()]

    currentVec[pos] = 0
    Theta = calcThetaMatrixHelper(dim, order, gamma, pos+1, currentVec)
    for q in range(1, order):
        currentVec[pos] = gamma**q
        Theta += calcThetaMatrixHelper(dim, order, gamma, pos+1, currentVec)

    return Theta

def convertToInteger(oa, order):
    gf = GF(order, name="a", modulus="primitive")
    gamma = gf.gen()
    orderedGf = [0] + [gamma**(val) for val in range(order-1)]
    # Hasmap to convert elements of GF(q) to corresponding integers
    gfToInt = {}
    for i, val in enumerate(orderedGf):
        gfToInt[val] = i+1

    intOa = []
    for i in range(oa.nrows()):
        row = []
        for j in range(oa.ncols()):
            row.append(gfToInt[oa[i][j]])
        intOa.append(row)
    return Matrix(intOa)

def isCorrectOA(oa, strength):
    vectorIndexSet = gmat.genVectorIndex([i for i in range(len(oa[0]))], strength)
    for vis in vectorIndexSet:
        tupleCount = {}
        for row in range(oa.nrows()):
            rowTuple = tuple([oa[row][val] for val in vis])
            if rowTuple in tupleCount:
                tupleCount[rowTuple] += 1
            else:
                tupleCount[rowTuple] = 1

        count = 0
        for v in tupleCount.values():
            if count == 0:
                count = v
            elif count != v:
                return False
    return True

def main(argc, argv):
    if argc != 4:
        print("usage: " + argv[0] + " dim order strength")
        return

    dim = int(argv[1])
    order = int(argv[2])
    strength = int(argv[3])
    if order < 2 or dim < 1 or strength < 2:
        print("error: invalid parameter:")
        print("    dim={}, order={}, strength={}".format(dim, order, strength))
        return
    if dim + 1 < strength:
        print("error: dim+1 should not be less than strength.")
        return

    G = gmat.calcGMatrix(dim, order, strength)
    print("G matrix:")
    print(G)
    Theta = calcThetaMatrix(dim+1, order)
    print("Theta matrix:")
    print(Theta)
    oa = G * Theta
    oa = oa.transpose()
    intOa = convertToInteger(oa, order)
    intOa = Matrix(sorted(intOa))
    print("Resulting OA:")
    print(intOa)

    print("=== Sanity check ===")
    if isCorrectOA(intOa, strength):
        print("Correct OA.")
    else:
        print("Invalid OA.")

if __name__ == "__main__":
    main(len(sys.argv), sys.argv)

$\Theta$ を求めるところ（calcThetaMatrix関数）がややこしいが、単純なbacktrackアルゴリズムで生成している。

以下に実行例を示す。

shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py
usage: ./pg.py dim order strength
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 2 2 2
G matrix:
[0 0 1]
[1 0 1]
[0 1 1]
[1 1 1]
[0 1 0]
[1 1 0]
[1 0 0]
Theta matrix:
[0 0 0 0 1 1 1 1]
[0 0 1 1 0 0 1 1]
[0 1 0 1 0 1 0 1]
Resulting OA:
[1 1 1 1 1 1 1]
[1 1 2 2 2 2 1]
[1 2 1 2 1 2 2]
[1 2 2 1 2 1 2]
[2 1 1 2 2 1 2]
[2 1 2 1 1 2 2]
[2 2 1 1 2 2 1]
[2 2 2 2 1 1 1]
=== Sanity check ===
Correct OA.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 2 2 3
G matrix:
[0 0 1]
[1 0 1]
[0 1 1]
[1 1 1]
Theta matrix:
[0 0 0 0 1 1 1 1]
[0 0 1 1 0 0 1 1]
[0 1 0 1 0 1 0 1]
Resulting OA:
[1 1 1 1]
[1 1 2 2]
[1 2 1 2]
[1 2 2 1]
[2 1 1 2]
[2 1 2 1]
[2 2 1 1]
[2 2 2 2]
=== Sanity check ===
Correct OA.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 3 2 2
G matrix:
[0 0 0 1]
[1 0 0 1]
[0 1 0 1]
[1 1 0 1]
[0 0 1 1]
[1 0 1 1]
[0 1 1 1]
[1 1 1 1]
[0 0 1 0]
[1 0 1 0]
[0 1 1 0]
[1 1 1 0]
[0 1 0 0]
[1 1 0 0]
[1 0 0 0]
Theta matrix:
[0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1]
[0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1]
[0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1]
[0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]
Resulting OA:
[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
[1 1 1 1 2 2 2 2 2 2 2 2 1 1 1]
[1 1 2 2 1 1 2 2 1 1 2 2 2 2 1]
[1 1 2 2 2 2 1 1 2 2 1 1 2 2 1]
[1 2 1 2 1 2 1 2 1 2 1 2 1 2 2]
[1 2 1 2 2 1 2 1 2 1 2 1 1 2 2]
[1 2 2 1 1 2 2 1 1 2 2 1 2 1 2]
[1 2 2 1 2 1 1 2 2 1 1 2 2 1 2]
[2 1 1 2 1 2 2 1 2 1 1 2 2 1 2]
[2 1 1 2 2 1 1 2 1 2 2 1 2 1 2]
[2 1 2 1 1 2 1 2 2 1 2 1 1 2 2]
[2 1 2 1 2 1 2 1 1 2 1 2 1 2 2]
[2 2 1 1 1 1 2 2 2 2 1 1 2 2 1]
[2 2 1 1 2 2 1 1 1 1 2 2 2 2 1]
[2 2 2 2 1 1 1 1 2 2 2 2 1 1 1]
[2 2 2 2 2 2 2 2 1 1 1 1 1 1 1]
=== Sanity check ===
Correct OA.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./pg.py 2 4 3
G matrix:
[    0     0     1]
[    a     0     1]
[    0     a     1]
[    a     a     1]
[    a     1     0]
[a + 1     1     0]
Theta matrix:
[    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     a     a     a     a     a     a     a     a     a     a     a     a     a     a     a     a a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1 a + 1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1]
[    0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1     0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1     0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1     0     0     0     0     a     a     a     a a + 1 a + 1 a + 1 a + 1     1     1     1     1]
[    0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1     0     a a + 1     1]
Resulting OA:
[1 1 1 1 1 1]
[1 1 2 2 4 4]
[1 1 3 3 2 2]
[1 1 4 4 3 3]
[1 2 1 2 2 3]
[1 2 2 1 3 2]
[1 2 3 4 1 4]
[1 2 4 3 4 1]
[1 3 1 3 3 4]
[1 3 2 4 2 1]
[1 3 3 1 4 3]
[1 3 4 2 1 2]
[1 4 1 4 4 2]
[1 4 2 3 1 3]
[1 4 3 2 3 1]
[1 4 4 1 2 4]
[2 1 1 2 3 2]
[2 1 2 1 2 3]
[2 1 3 4 4 1]
[2 1 4 3 1 4]
[2 2 1 1 4 4]
[2 2 2 2 1 1]
[2 2 3 3 3 3]
[2 2 4 4 2 2]
[2 3 1 4 1 3]
[2 3 2 3 4 2]
[2 3 3 2 2 4]
[2 3 4 1 3 1]
[2 4 1 3 2 1]
[2 4 2 4 3 4]
[2 4 3 1 1 2]
[2 4 4 2 4 3]
[3 1 1 3 4 3]
[3 1 2 4 1 2]
[3 1 3 1 3 4]
[3 1 4 2 2 1]
[3 2 1 4 3 1]
[3 2 2 3 2 4]
[3 2 3 2 4 2]
[3 2 4 1 1 3]
[3 3 1 1 2 2]
[3 3 2 2 3 3]
[3 3 3 3 1 1]
[3 3 4 4 4 4]
[3 4 1 2 1 4]
[3 4 2 1 4 1]
[3 4 3 4 2 3]
[3 4 4 3 3 2]
[4 1 1 4 2 4]
[4 1 2 3 3 1]
[4 1 3 2 1 3]
[4 1 4 1 4 2]
[4 2 1 3 1 2]
[4 2 2 4 4 3]
[4 2 3 1 2 1]
[4 2 4 2 3 4]
[4 3 1 2 4 1]
[4 3 2 1 1 4]
[4 3 3 4 3 2]
[4 3 4 3 2 3]
[4 4 1 1 3 3]
[4 4 2 2 2 2]
[4 4 3 3 4 4]
[4 4 4 4 1 1]
=== Sanity check ===
Correct OA.

上記のプログラムでは $G \Theta$ を転置して $\Theta^t G^t$ の計算結果を出力している点、及び入力する次元の値は射影空間の次元の値としている (つまり $d-1$ の値を入力する仕様である) 点に注意されたい。

また、有限体の元をそのまま出力すると分かりづらいので、以前の記事と同じ方法で整数値に変換している。その後、行を昇順にsortしたものを最終出力としている。

今回も最後に正しく直交表が生成できているかチェックしているが、全てパスしている。

ちなみに

最近気づいたのだが、実はsageには直交表を直接的に求めるための関数が存在している[2]。「なんだ、こんな便利な関数があるのか」とガックリしたのだが、少し触ってみたところ、指数2以上の直交表は構築できないような挙動をしている。単なる私の理解不足かもしれないが、なんだかいまいちな仕様であり、本稿で作ったプログラムも無意味ではなかったかなと思っている。

以下に実行例を示す。第1引数が因子数、第2引数が水準数で、 $t$ というパラメータに強さの値を与えればよい。

sage: for row in designs.orthogonal_arrays.build(2, 2, t=2):
....:     print(row)
....:
[0, 0]
[0, 1]
[1, 0]
[1, 1]
sage: for row in designs.orthogonal_arrays.build(3, 2, t=2):
....:     print(row)
....:
[0, 0, 0]
[0, 1, 1]
[1, 0, 1]
[1, 1, 0]
sage: for row in designs.orthogonal_arrays.build(4, 2, t=2):
....:     print(row)
....:
---------------------------------------------------------------------------
EmptySetError                             Traceback (most recent call last)
<ipython-input-156-349912d12039> in <module>
----> 1 for row in designs.orthogonal_arrays.build(Integer(4), Integer(2), t=Integer(2)):
      2     print(row)
      3

/usr/lib/python3/dist-packages/sage/combinat/designs/orthogonal_arrays.py in build(k, n, t, resolvable)
   2142
   2143         """
-> 2144         return orthogonal_array(k,n,t,resolvable=resolvable)
   2145
   2146     @staticmethod

/usr/lib/python3/dist-packages/sage/combinat/designs/orthogonal_arrays.py in orthogonal_array(k, n, t, resolvable, check, existence, explain_construction)
    891         if explain_construction:
    892             return msg
--> 893         raise EmptySetError(msg)
    894
    895     elif k <= t:

EmptySetError: There exists no OA(4,2) as k(=4)>n+t-1=3

まとめ

本稿では有限射影幾何の理論を用いた直交表構築手法の証明を行った。その後、実際に直交表が構築できることをsageを使って確かめた。

もともとは単に直交表の生成アルゴリズムを調べていただけなのだが、有限幾何という非常に面白い分野に出くわしたのは幸運であった。せっかくなので、もう少し有限幾何で遊んでみたいと思う。

また、改めてsageの凄さを思い知らされた。計算機上で高度な数学的な計算が実行できるというのは、確実に数学の研究を加速させる有用なテクノロジーだと思う。こういった分野の歴史や手法について調べてみるのも面白いかもしれない。

参考

[1] tasusu.hatenablog.com
[2] doc.sagemath.org

2021-04-20

直交表の数理～有限射影幾何を用いたG行列の構築編～

前回の記事ではMOLSを用いた直交表の構築方法について説明した。この方法では構築できる直交表のパターンに限界があるため、他のより優れた方法について調べてみた。

直交表の構築方法にはいろいろあるようだが、本稿では私が特に数学的に面白そうだと感じた有限射影幾何を用いた方法について説明したいと思う。ただし、有限射影幾何とはどういったものかという説明だけでもなかなかの文章量になりそうなので、やはり記事をいくつかに分けて説明したいと思う。その走りとして、本稿では手法の概要と有限射影幾何の基礎について述べた後、直交表の構築に必要な $G$ 行列というものを構築する方法について説明したいと思う。

$G$ 行列を用いた直交表の構築手法

手法の概要

論文[1]によると、強さ $t\ (t=2, 3, 4)\$ の直交表を作るためには、ある性質を持つ行列 $G$ があればよいらしい。その性質とは「任意の $t$ 行を選んだとき、それらが一次独立になっている」ということである。もしそのような行列 $G$ が得られれば、以下のようにして直交表が生成できるとのことだ。

$\displaystyle{ G\Theta }$

ただし、 $k$ を因子数とすると、 $G$ は $k$ 行 $d$ 列の行列、 $\Theta$ は $d$ 行 $q^d$ 列の行列である。特に、 $\Theta$ は $\mathbb{F}_q$ 上の $d$ 次元ベクトル空間 $\mathbb{F}_q^d$ の相異なる元を列ベクトルとして $q^d$ 本並べることで得られる行列である。 $d$ は次元のパラメータであるが、詳細は後述する。

生成される行列 $G\Theta$ は $k$ 行 $q^d$ 列となる。この行列を転置すれば、前回の記事と同様のフォーマットの直交表が得られる。すなわち、 $(G\Theta)^t = \Theta^t G^t$ が求める直交表となる。

疑問点のオンパレードと解決の筋道

さて、これだけ見ると疑問点だらけでわけが分からない。少なくとも、以下のような疑問が湧いてくる。

$G$ 行列はどのようにして求めればよいのか？
$d$ というパラメータがいきなり出てきたが、これは何に基づいてどういう値に設定すべきものなのか？
なぜ $G \Theta$ が直交表になると言えるのか？
構築される直交表はどういうパターンになるか？

これらの疑問について、まずは現時点で私が理解していることの概要について述べる。

まず1点目について、これこそがまさに有限射影幾何が威力を発揮するところである。具体的には、有限体上の射影空間の点を求めて、それらから必要な部分集合を抜き出すことで解決できる。本稿の主題であり、難しいポイントなので詳細は後述する。

2点目について、dは因子数や水準数に依存して決まる。詳細は後述する。

3点目について、これは現時点で証明が見つけられていない。論文[1]で引用されている書籍に説明があるようだが、手元にないため参照できない。というわけで、今回も証明を自力で考えてみようと思うが、現時点での見通しはなく証明できるか分からない。少なくとも、これについては次回以降に回す。

4点目について、まず強さは2～4までいけるようである。5以上でも同様の理屈で行けそうだが、強さに応じて1点目で述べた「射影空間の点から必要な部分集合」なるものを抜き出すステップの計算量が爆発していくのと、実用上あまり必要になるとは思えないため、強さ5以上は考えないことにする。

有限体を利用している関係上、水準数は相変わらず素数冪に限られる。しかし、因子数は自由に選べるというのがこの手法の利点となる。

結局、構築できる直交表の型は $t\text{-}(q, k, \lambda)$ ( $2 \le t \le 4$ , $q$ は素数冪, $1 \le k, 1 \le \lambda)$ となる。 $\lambda$ は他のパラメータに応じて自動的に決まるため、無駄に大きくすることはできるが入力としての自由度は基本的にないものと考えてよい。

有限体上の射影空間

では、 $G$ 行列を求めるための有限射影幾何について考えていこう。具体的には有限体上の射影空間を考えるのだが、先に射影空間の一般論から説明する。

射影空間

射影空間の定義はいろいろなバリエーションがあるが、今回使用するものをWikipedia[2]から引用する。

Projective space

Given a vector space $V$ over a field $K$ , the projective space $P(V)$ is the set of equivalence classes of $V \setminus \{0\}$ under the equivalence relation $\sim$ defined by $x \sim y$ if there is a nonzero element $\lambda$ of $K$ such that $x = \lambda y$ .

つまり射影空間とは、互いにスカラー倍の関係にあるベクトル同士を同一視して、それらを改めて点と捉えた空間のことである。元のベクトル空間 $V$ の次元が $n+1$ なのに対して、そこから得られる射影空間の次元は $n$ であることに注意されたい。

また、射影空間における直線というものも考えることができる。まず、元の空間において原点を通る2つの異なる直線 $l_1, l_2$ があったとして、それらを両方とも含む平面 $S$ を考える。射影空間上では $l_1, l_2$ は点になるが、この平面 $S$ が射影空間上では点 $l_1, l_2$ を結ぶ直線となる。

もとのベクトル空間上での直線を点、平面を直線と捉えるため、少々ややこしい空間であるが、こういうのは慣れるしかない。

有限体の場合

さて、射影空間を得るためにはまずベクトル空間を考える必要があるが、体 $K$ としては何を選んでもよい。つまり、有限体上のベクトル空間であっても良いわけである。有限体上の射影空間のことを有限射影空間と呼ぶ。

以下では $\mathbb{F}_q$ ( $q$ は素数冪) 上の $n+1$ 次元ベクトル空間から得られる $n$ 次元射影空間を $\mathrm{PG}(n, q)$ と表記する。

有限射影空間の完全代表系

有限射影空間上の全ての点から完全代表系を求め、それを再び元のベクトル空間 $V$ の元と考えると、任意の2点から得られる $V$ 上のベクトルは一次独立となる。これは有限射影空間の定義を考えれば当然である。

有限射影空間の存在条件

有限射影空間は $q$ が素数冪の場合には存在することが知られているが、そうでない場合に存在するかどうかは未解決問題である[4]。そのため、以降の議論では $q$ が素数冪の場合のみを扱う。

有限射影空間上の点の数

$\mathbb{F}_q$ 上の $n$ 次元ベクトル空間 $V$ における点 (ベクトル) の数は $q^n$ 個に限られる。さらに、 $\mathrm{PG}(n, q)$ の場合は互いにスカラー倍の関係にある点は同一視されるので、さらに点の数が減ることになる。これは具体的にはいくつになるのだろうか？

次数が2の場合、つまり $\mathrm{PG}(2, q)$ 上の点の数は $q^2+q+1$ 個になるようである[3]。ちなみに、2次元の射影空間のことを特別に射影平面と呼ぶ。

次数が3以上の場合はどうなるのだろうか？これはGaussian binomial coefficientと呼ばれるものを使って以下のように算出できる[4][5]。

$\displaystyle{ \left( \begin{array}{c} n+1 \\ 1 \\ \end{array} \right)_q }$

私自身、Gaussian binomial coefficientなどというものは初めて見たのだが、計算の詳細は本稿の主眼ではないため、[5]などを参照されたい。ここではいくつかの具体的な $n, q$ について、 $\mathrm{PG}(n, q)$ 上の点の数がどうなるかをsageで計算してみよう。

sage: for dim in range(2,8):
....:     for order in range(2,8):
....:         if(is_prime_power(order)):
....:             print("The number of points in PG({}, {}) = {}".format(dim, order, q_binomial(dim+1, 1, q=order)))
....:
The number of points in PG(2, 2) = 7
The number of points in PG(2, 3) = 13
The number of points in PG(2, 4) = 21
The number of points in PG(2, 5) = 31
The number of points in PG(2, 7) = 57
The number of points in PG(3, 2) = 15
The number of points in PG(3, 3) = 40
The number of points in PG(3, 4) = 85
The number of points in PG(3, 5) = 156
The number of points in PG(3, 7) = 400
The number of points in PG(4, 2) = 31
The number of points in PG(4, 3) = 121
The number of points in PG(4, 4) = 341
The number of points in PG(4, 5) = 781
The number of points in PG(4, 7) = 2801
The number of points in PG(5, 2) = 63
The number of points in PG(5, 3) = 364
The number of points in PG(5, 4) = 1365
The number of points in PG(5, 5) = 3906
The number of points in PG(5, 7) = 19608
The number of points in PG(6, 2) = 127
The number of points in PG(6, 3) = 1093
The number of points in PG(6, 4) = 5461
The number of points in PG(6, 5) = 19531
The number of points in PG(6, 7) = 137257
The number of points in PG(7, 2) = 255
The number of points in PG(7, 3) = 3280
The number of points in PG(7, 4) = 21845
The number of points in PG(7, 5) = 97656
The number of points in PG(7, 7) = 960800

さらっと計算したが、sageは本当になんでも計算できてすごい。上記の結果から、 $n, k$ のいずれが増加しても点の数が爆発的に増加する様子が見て取れる。

有限射影空間と $G$ 行列

$G$ 行列の求め方

強さ2の場合

さて、話を直交表に戻そう。有限射影空間を用いることで、直ちに強さ2の直交表を求めることができる。というのも、前述の通り射影空間の完全代表系を元のベクトル空間に戻してやれば、そのうちの任意の2ベクトルは互いに一次独立となっており、これはまさに $G$ 行列の要件を満たすからである。

つまり、射影空間の完全代表系を元のベクトル空間に戻したものを行ベクトルとして縦に並べれば、それが $G$ 行列となる。

強さ3以上の場合

強さ3以上の場合は、まず強さ2の $G$ 行列を構築する。その上で、そこから任意の $t$ 行を選んだ時に、それらが一次独立になるように一生懸命探索して部分集合を選別するしかない。これをまじめにやると計算量が指数オーダーになってしまうので、以下のようなgreedyなアルゴリズムを使うなど、計算の工夫が必要となる。

選抜したベクトルの集合を空集合で初期化する。
強さ2の $G$ 行列の行を順に選んで以下を実行する。
今注目しているベクトルと現時点で選別されたベクトルの集合を合わせた集合について、そこから任意の $t$ 個のベクトルを選ぶと一次独立になるかを調べる。YESであれば今のベクトルを選別したベクトルの集合に加える。

ただし、これだと $G$ 行列の行数を最大化できないということが普通に起こりえる。また、それでなくとも完全代表系の選び方によっては $G$ 行列の行数を最大化できないケースがあるのではないのか？という疑問もある。しかし、実用上は $G$ の行数は多ければ多いほどよいが、必ずしも最大の行数となる必要はない。そのため、手法や計算の容易性のために $G$ 行列が多少小さくなることは許容できるものと考えている。

計算してみよう

射影空間上の点、及びその完全代表系はsageを使って計算が可能である。また、ベクトルの集合が一次独立であるかどうかもsageを使えば簡単に確かめられる。以下に $G$ 行列を求めるプログラムを示す。

#!/usr/bin/env sage

import sys
from sage.all import *

def getProjectiveSpacePoints(dim, order):
    PG = ProjectiveSpace(dim, GF(order, name="a", modulus="primitive"))
    return PG.rational_points()

def calcGMatrix(dim, order, strength):
    points = getProjectiveSpacePoints(dim, order)
    # Convert points in PG(dim, order) to vectors in vector space over F_{order}
    vectors = [vector(p) for p in points]
    if strength == 2:
        return Matrix(vectors)
    else:
        return Matrix(getLinearIndependentVectors(vectors, order, strength))

def getLinearIndependentVectors(vectors, order, strength):
    V = VectorSpace(GF(order, name="a", modulus="primitive"), len(vectors[0]))
    currentSet = []
    for v in vectors:
        if isRowIndependent(Matrix([v] + currentSet), order, strength):
            currentSet += [v]
    return currentSet

def genVectorIndex(index, numVecs):
    vectorsSet = []
    if len(index) == numVecs:
        vectorsSet.append([ind for ind in index])
    else:
        for i in range(len(index)-numVecs+1):
            if numVecs == 1:
                vectorsSet += [[index[i]]]
            else:
                subVectorsSet = genVectorIndex(index[i+1:], numVecs-1)
                for j in range(len(subVectorsSet)):
                    subVectorsSet[j] = [index[i]] + subVectorsSet[j]
                vectorsSet += subVectorsSet
    return vectorsSet

def isRowIndependent(G, order, strength):
    V = VectorSpace(GF(order, name="a", modulus="primitive"), len(G[0]))
    # If the number of rows of G matrix is relatively small,
    # all we have to do is just to check the linear independency of all rows of G.
    if G.nrows() <= strength:
        if V.linear_dependence([v for v in G]) != []:
            return False
    else:
        # If the dimension of the vector space is too small,
        # it is impossible to choose linear independent vector set.
        if G.ncols() < strength:
            return False
        vectorIndexSet = genVectorIndex([i for i in range(G.nrows())], strength)
        for vis in vectorIndexSet:
            vs = [G[i] for i in vis]
            if V.linear_dependence(vs) != []:
                return False
    return True

def main(argc, argv):
    if argc != 4:
        print("usage: " + argv[0] + " dim order strength")
        return

    dim = int(argv[1])
    order = int(argv[2])
    strength = int(argv[3])
    if order < 2 or dim < 1 or strength < 2:
        print("error: invalid parameter:")
        print("    dim={}, order={}, strength={}".format(dim, order, strength))
        return
    if dim + 1 < strength:
        print("error: dim+1 should not be less than strength.")
        return

    G = calcGMatrix(dim, order, strength)
    print("G matrix:")
    print(G)

    print("=== Sanity check ===")
    if isRowIndependent(G, order, strength):
        print("Correct G matrix.")
    else:
        print("Invalid G matrix.")

if __name__ == "__main__":
    main(len(sys.argv), sys.argv)

以下に実行例を示す。

shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py
usage: ./gmat.py dim order strength
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py 2 2 2
G matrix:
[0 0 1]
[1 0 1]
[0 1 1]
[1 1 1]
[0 1 0]
[1 1 0]
[1 0 0]
=== Sanity check ===
Correct G matrix.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py 2 2 3
G matrix:
[0 0 1]
[1 0 1]
[0 1 1]
[1 1 1]
=== Sanity check ===
Correct G matrix.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py 2 2 4
error: dim+1 should not be less than strength.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./gmat.py 3 3 4
G matrix:
[0 0 0 1]
[1 0 0 1]
[0 1 0 1]
[0 0 1 1]
[1 1 1 1]
=== Sanity check ===
Correct G matrix.

実行時引数として受け取っている次元の値は射影空間の次元としている。これは前述のパラメータ $d$ よりも1小さい値となっているので注意されたい。

得られた $G$ 行列は最後に任意の $t$ 行が一次独立になるかのチェックを行っているが、それは全てパスしている。

次元 $d$ の決め方

前述した疑問点の2つ目についてここで詳細に考えてみよう。直交表を作るにあたり、各因子は $G$ 行列の行に対応することになる。行を削る分には構わないが、行が足りないとなった場合には、次元を上げて行数を増やす必要がある。そのため、 $d, q$ の値をいろいろと変えたときに $G$ 行列の行数がどうなるかをあらかじめ調べておき、テスト要件として求められる因子数を満たすような値を指定すればよい。

例として、因子数6、水準数2、強さ3の直交表を構築したいとする。先ほど示したプログラムで計算すると、 $d=2, q=2, t=3$ の $G$ 行列の行数は4、 $d=3, q=2, t=3$ の $G$ 行列の行数は8となる。そのため、次数として2では足りず、3であればOKということになる。

まとめ

本稿では有限体上の射影空間の理論を用いた直交表の構築手法について説明した。特に、直交表を構築する上で最も重要となる $G$ 行列の構築方法にフォーカスした。

本稿で述べた方法によりなぜ直交表が構築できるのかという部分について疑問が残っているため、それについてはできる限り次回に解決するようにしたい。また、最終的に直交表を構築するところまでプログラムを紹介できなかったため、それも次回に説明したいと思う。

参考

[1] 須田健二, 五味弘, 有限射影幾何を用いたソフトウェアテスト向けの直交表自動生成プログラムの開発とその応用, 2016, https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiAtemPvonwAhVKxYsBHa7IDUwQFjABegQIAxAD&url=https%3A%2F%2Fipsj.ixsq.nii.ac.jp%2Fej%2F%3Faction%3Drepository_uri%26item_id%3D174252%26file_id%3D1%26file_no%3D1&usg=AOvVaw2fjazMjJ9ADLdgt_Y1nTAs
[2] Projective space - Wikipedia
[3] Projective plane - Wikipedia
[4] Finite geometry - Wikipedia
[5] Gaussian binomial coefficient - Wikipedia

2021-04-16

直交表の数理～ MOLSを用いた直交表の構築編～

実験計画法

前回の記事ではmutually orthogonal Latin squares (MOLS) というテーブル群を構築する方法について説明した。本稿ではMOLSを用いて直交表を構築する方法について説明する。

直交表の構築

基本の構築方法

$p$ を素数、 $k$ を1以上の整数、 $q = p^k$ とする。このとき、MOLS( $q$ )からは $2\text{-}(q, m, 1)$ ( $m = 2, \cdots, q+1$ ) というタイプの直交表が得られる。以下ではMOLS(3)から直交表 $2\text{-}(3, 4, 1)$ を得るケースを例にして、直交表の構築方法を説明する。全体的に英語版Wikipedia [1]を参考にした。

まず、直交表の最初の2列を順序対が辞書順になるように埋める。

#	因子1	因子2
1	1	1
2	1	2
3	1	3
4	2	1
5	2	2
6	2	3
7	3	1
8	3	2
9	3	3

次に、今埋めた2列の値をそれぞれLatin Squareの行・列と見なして、MOLSに属する各表から値を取り出す。今は例としてMOLS(3)を考えているので、これに属するLatin squareは $L_1, L_2$ の2つである。そのため、3列目と4列目をそれぞれ $L_1(i, j), L_2(i, j)$ の値で埋めればよい。

分かりやすさのために $L_1, L_2$ を具体的に示しておく。

$L_1$

1	2	3
2	3	1
3	1	2

$L_2$

1	2	3
3	1	2
2	3	1

これより、以下のように直交表が得られる。

#	因子1	因子2	因子3	因子4
1	1	1	1	1
2	1	2	2	2
3	1	3	3	3
4	2	1	2	3
5	2	2	3	1
6	2	3	1	2
7	3	1	3	2
8	3	2	1	3
9	3	3	2	1

直交表が得られていることの証明

なんだかよくわからないうちに直交表が構築できてしまったように見えるのだが、これは本当に直交表になっているのだろうか？相変わらずWikipedia[1]には証明がないので、自分で考えてみた。なお、記述をシンプルにするため、これ以降はLatin Squareや直交表の行・列の番号は0始まりとする。

上記の方法によって得られるテーブルの行数は $q^2$ となる。そのため、（強さ2で考えているので）任意の2列を選んで各行から順序対を形成したときに、それらがすべて異なることを示せばよい。そうすれば、ちょうど異なる $q^2$ 個の順序対が出尽くすことになるため、指数1の直交表が得られたことになる。

0列目、1列目、及び2列目以降で性質が異なるので、以下の組み合わせのパターンで場合分けをする。

0列目と1列目
0列目と2列目以降
1列目と2列目以降
2列目以降同士

準備

各ケースでの証明を始める前に、それぞれの列の $i$ 行目の値がどのように書き表せるかを最初にまとめておこう。有限体 $\mathbb{F}_q$ の元は前回の記事に倣って $a_0, \cdots, a_{q-1}$ と書く。

0列目： $\left\lfloor \frac{i}{q} \right\rfloor$ +1
1列目： $(i \bmod q)$ +1
$r+1$ ( $1 \le r \le q-1$ ) 列目： $a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor}$

0列目と1列目

このケースは順序対が辞書順になるように値を決めているので自明である。

0列目と2列目以降

1列目と $r+1$ ( $1 \le r \le q-1$ ) 列目について順序対を形成したとき、 $i$ 行目と $i'$ 行目の順序対が一致したと仮定する。すると、以下の2つの等式が成立する。

$\displaystyle{ \begin{eqnarray} && \left\lfloor \frac{i}{q} \right\rfloor = \left\lfloor \frac{i'}{q} \right\rfloor \\ && a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} = a_{(i' \bmod q)} + a_r a_{\left\lfloor \frac{i'}{q} \right\rfloor} \end{eqnarray} }$

1つ目の式を2つ目の式に代入して整理すると以下のようになる。

$\displaystyle{ \begin{eqnarray} a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} &=& a_{(i' \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} \\ a_{(i \bmod q)} &=& a_{(i' \bmod q)} \\ \end{eqnarray} }$

先ほどの1つ目の式より $|i-i'| < q$ となるため、 $i \equiv i' \pmod q$ となるためには $i = i'$ となる他ない。

1列目と2列目以降

2列目と $r+1$ ( $1 \le r \le q-1$ ) 列目について順序対を形成したとき、 $i$ 行目と $i'$ 行目の順序対が一致したと仮定する。すると、以下の2つの等式が成立する。

$\displaystyle{ \begin{eqnarray} && (i \bmod q) = (i' \bmod q) \\ && a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} = a_{(i' \bmod q)} + a_r a_{\left\lfloor \frac{i'}{q} \right\rfloor} \end{eqnarray} }$

1つ目の式を2つ目の式に代入して整理すると以下のようになる。

$\displaystyle{ \begin{eqnarray} a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i}{q} \right\rfloor} &=& a_{(i \bmod q)} + a_r a_{\left\lfloor \frac{i'}{q} \right\rfloor} \\ a_{\left\lfloor \frac{i}{q} \right\rfloor} &=& a_{\left\lfloor \frac{i'}{q} \right\rfloor} \end{eqnarray} }$

2つ目の等号では $r\ne 0$ より $a_r \ne 0$ であることを用いた。

先ほどの1つ目の式より $|i-i'|$ は $q$ の倍数となるため、 $\left\lfloor \frac{i}{q} \right\rfloor = \left\lfloor \frac{i'}{q} \right\rfloor$ となるためには $i = i'$ となる他ない。

2列目以降同士

2列目以降同士の組み合わせであれば、MOLSの定義より互いに直交するLatin squareから順序対を作ることになるため、それらは当然すべて異なる。

以上により、前述の方法で直交表が構築できることが分かった。

MOLSの限界とその他の手法

因子数の減らし方

前述した直交表の構築方法では、MOLSの完全集合の要素をすべて使用した。この場合、因子数が $q+1$ に限定されてしまい、生成できる直交表のパターンがあまりにも限定的である。

実は、因子数を $q+2$ 以上に上げることは難しいが、下げることは割と簡単である。単にMOLSの要素から必要な分だけ排除してやれば良いだけである。例えば、先ほど直交表 $2\text{-}(3, 4, 1)$ を構築したが、これを $2\text{-}(3, 3, 1)$ にしたければ単に $L_2$ を無視して最終列を削れば良いだけである。

ただし、この場合でも行数を削ることはできない。これは、因子数が3だろうが4だろうが、例えば最初の2行の順序対のパターンを出し尽くすのに必要な行数は確保する必要があるためである。

MOLSの限界

結局、MOLS( $q$ )から生成できる直交表のパターンは前述のとおり $2\text{-}(q, m, 1)$ ( $m = 2, \cdots, q+1$ )となる。これがMOLSを用いた直交表構築法の限界である。具体的には、以下のような直交表は生成することができない。

強さが3以上
水準数が素数冪でない
因子数が $q+2$ 以上

その他の直交表構築法

MOLSを用いる方法以外にも、直交表を構築する方法はあるようだ。詳しくは英語版Wikipedia[2]にいろいろ書いてある。例えば、Latin squareの代わりにLatin cubeとかLatin hypercubeと呼ばれるものを使う方法、アダマール行列なるものを使う方法、及び何やら難しい符号化の理論を使う方法があるらしい。いずれについても、現時点で私は理解できていない。

直交表を構築してみよう

では、実際に直交表を構築するプログラムを書いてみよう。前回の記事で作成したmols.pyを部品として使用し、直交表を構築するプログラム例を以下に示す。

#!/usr/bin/python3

import mols
import sys

def calcOrthogonalArray(ml, factor):
    order = len(ml[0])
    height = order**2
    oa = [[0 for i in range(factor)] for j in range(height)]
    for i in range(height):
        for j in range(factor):
            if j == 0:
                oa[i][j] = int(i/order) + 1
            elif j == 1:
                oa[i][j] = i%order + 1
            else:
                latSq = ml[j-2]
                oa[i][j] = latSq[int(i/order)][i%order]
    return oa

def main(argc, argv):
    if argc != 3:
        print("usage: " + argv[0] + " factor level")
        return

    factor = int(argv[1])
    level = int(argv[2])
    oa = []
    if factor == 1:
        # This is a special case
        oa = [[i+1] for i in range(level)]
    elif factor > level+1:
        print("error: factor should be less than or equal to level+1.")
    else:
        ml = mols.calcMOLS(level)
        oa = calcOrthogonalArray(ml, factor)

    for line in oa:
        print(line)

if __name__ == "__main__":
    main(len(sys.argv), sys.argv)

このプログラムの実行例を以下に示す。

shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 1 1
[1]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 1 3
[1]
[2]
[3]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 2 3
[1, 1]
[1, 2]
[1, 3]
[2, 1]
[2, 2]
[2, 3]
[3, 1]
[3, 2]
[3, 3]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 3 3
[1, 1, 1]
[1, 2, 2]
[1, 3, 3]
[2, 1, 2]
[2, 2, 3]
[2, 3, 1]
[3, 1, 3]
[3, 2, 1]
[3, 3, 2]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 4 3
[1, 1, 1, 1]
[1, 2, 2, 2]
[1, 3, 3, 3]
[2, 1, 2, 3]
[2, 2, 3, 1]
[2, 3, 1, 2]
[3, 1, 3, 2]
[3, 2, 1, 3]
[3, 3, 2, 1]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 5 3
error: factor should be less than or equal to level+1.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 6 5
[1, 1, 1, 1, 1, 1]
[1, 2, 2, 2, 2, 2]
[1, 3, 3, 3, 3, 3]
[1, 4, 4, 4, 4, 4]
[1, 5, 5, 5, 5, 5]
[2, 1, 2, 3, 4, 5]
[2, 2, 3, 5, 1, 4]
[2, 3, 5, 4, 2, 1]
[2, 4, 1, 2, 5, 3]
[2, 5, 4, 1, 3, 2]
[3, 1, 3, 4, 5, 2]
[3, 2, 5, 1, 4, 3]
[3, 3, 4, 2, 1, 5]
[3, 4, 2, 5, 3, 1]
[3, 5, 1, 3, 2, 4]
[4, 1, 4, 5, 2, 3]
[4, 2, 1, 4, 3, 5]
[4, 3, 2, 1, 5, 4]
[4, 4, 5, 3, 1, 2]
[4, 5, 3, 2, 4, 1]
[5, 1, 5, 2, 3, 4]
[5, 2, 4, 3, 5, 1]
[5, 3, 1, 5, 4, 2]
[5, 4, 3, 1, 2, 5]
[5, 5, 2, 4, 1, 3]
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./oa.py 6 6
/usr/lib/python3/dist-packages/apport/report.py:13: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses
  import fnmatch, glob, traceback, errno, sys, atexit, locale, imp, stat
Traceback (most recent call last):
  File "./oa.py", line 42, in <module>
    main(len(sys.argv), sys.argv)
  File "./oa.py", line 35, in main
    ml = mols.calcMOLS(level)
  File "/mnt/c/Users/shiny/Dropbox/programs/blog/2021/mols.py", line 17, in calcMOLS
    gf = GF(q, name="a", modulus="primitive")
  File "sage/structure/factory.pyx", line 367, in sage.structure.factory.UniqueFactory.__call__ (build/cythonized/sage/structure/factory.c:2162)
  File "/usr/lib/python3/dist-packages/sage/rings/finite_rings/finite_field_constructor.py", line 566, in create_key_and_extra_args
    raise ValueError("the order of a finite field must be a prime power")
ValueError: the order of a finite field must be a prime power

第1引数が因子数、第2引数が水準数である。因子数が $q+1$ 以下でなければならなかったり、水準数が素数冪でなければならない制限が見て取れる*1。

まとめ

本稿ではMOLSから直交表を構築する方法について説明した。また、この方法で生成できる直交表のバリエーションの限界について述べた。

現時点では良く分かっていないのだが、MOLSは有限射影幾何なるものと関連があるらしい。次回はそれについて調べてまとめてみたいと思う。ただ、見通しはないのでどうなるかは分からない。調べた結果、私が理解できて、かつ面白いと思ったら記事を書こうと思う。

参考

[1] Mutually orthogonal Latin squares - Wikipedia
[2] Orthogonal array - Wikipedia

*1:水準数が素数冪でない場合のエラーハンドリングがなく、sageの関数内で異常終了しているが、サンプルプログラムなのでご容赦願いたい。

2021-04-16

直交表の数理～ MOLSの構築編～

実験計画法

世の中には実験計画法というものがある。実験計画法そのものに対して私は特段知見があるわけではないが、これがソフトウェア等のテストパターンを効率的に削減するのに役立つということは以前からなんとなく知っていた。考え方としては、全部のテスト条件の組み合わせをテストするのはパターンが多すぎてやりきれないので、その中から可能な限り互いに異なるものを選別することで、テストのバリエーションを確保しつつ項目数を削減しようというわけだ。

その選別方法にはいくつか手法があるが、その1つが直交表を用いた方法である。直交表の詳細な定義は後述するが、これはある種の制約を持ったテーブルであり、その制約を満たすパターンを見つけるのが難しい。私は当初この生成方法に興味を持ってなんとなく調べていたのだが、その背後には組み合わせ数学、代数学、幾何学など様々な数学との関連があることが分かった。

しかし、実験計画法というのがどちらかといえば実用寄りの話であり、（単なる推測だが）数学的な話題が主眼ではないためか、あまり直交表の数理的側面について詳しく解説された日本語の記事*1をインターネット上で見つけることができなかった。

そこで、本稿から続くいくつかの記事の中で、直交表の数理的側面について自分が理解したことをまとめてみようと思う。手始めに、本稿では直交表の基本的な説明を行った後、直交表を生成するために必要なmutually orthogonal Latin squares (MOLS) と呼ばれるテーブル群の生成方法について説明する。

直交表とは

定義

直交表の定義を日本語版のWIkipedia[1]から引用する。

直交表

直交表（ちょっこうひょう）とは、どの2列をとっても、その水準のすべての組み合わせが同数回現れる配列(表)のことである。これを「その2列はバランスしている」あるいは「直交している」と呼ぶ。この性質を使って、実験計画法において因子と水準の割り付け表に用いられる。

これはちょっと定義がいまいちである。というのも、列数を2つに限定しているが、実際には3つ以上ということもあり得るからだ。英語版のWikipedia[2]ではこの点も含めた定義になっているので、そちらも引用する。

Orthogonal array

In mathematics, an orthogonal array is a "table" (array) whose entries come from a fixed finite set of symbols (typically, $\{1,2,...,n\}$ ), arranged in such a way that there is an integer $t$ so that for every selection of $t$ columns of the table, all ordered $t$ -tuples of the symbols, formed by taking the entries in each row restricted to these columns, appear the same number of times.

直交表を決めるパラメータ

直交表にはいくつかパターンがある。そのパターンを決めるためのパラメータが複数あるのだが、それを明示するために直交表を $t\text{-}(v,k,\lambda)$ という記号で表す。これは英語版Wikipedia[2]の流儀であり、他にも記法は存在するが、本稿ではこれを採用する。

それぞれのパラメータの意味を以下に示す。

$v$ : 水準数
$k$ : 因子数
$t$ : 強さ
$\lambda$ : 指数

このとき、テストケース数は $\lambda v^t$ となる。つまり、直交表は $\lambda v^t$ 行 $k$ 列のテーブルとなる。

これだけだと分かりづらいので例を挙げよう。あるソフトウェアのテストをすることを考える。テスト条件は以下のようであったとする。

OS：Windows、Linux、MacOS
CPU：Intel、AMD、ARM
データベース：MySQL、SQLite、Oracle DB
ブラウザ：IE、Chrome、firefox

ここではテスト条件の項目が4つあるが、これが因子数である。また、各因子について選択肢が3つずつあるが、これが水準数となる。また、4つある因子の各水準の組み合わせとして、例えば異なる2つの組み合わせを網羅したいのであれば、強さは2となる。さらに、直交表ではそれらの組み合わせは必ず同じ回数ずつ出現するが、その回数を表すのが指数となる。例えば各組み合わせを1回ずつテストしたいとすると、上記のテストのための直交表は $2\text{-}(3,4,1)$ となる。

実際の直交表がどうなるかも見ておこう。以下の例はTech Note[3]というWebサイトを参考に作成してみた。

#	OS	CPU	データベース	ブラウザ
1	1	1	1	1
2	1	2	2	2
3	1	3	3	3
4	2	1	2	3
5	2	2	3	1
6	2	3	1	2
7	3	1	3	2
8	3	2	1	3
9	3	3	2	1

ただし、各水準に適当に1～3の数値を割り当てて考えるものとする。例えば、WIndows: 1、Linux: 2、MacOS: 3という感じで考えればよい。

直交表を用いない場合、テスト条件の組み合わせ数は $3^4 = 81$ 通りとなる。それに比べると大幅にテストケース数を削減できていることが分かるだろう。

強さのパラメータは直交表構築手法に重大な影響を与える。強さ3以上の直交表の構築方法は現時点では難しくて分からないので、しばらくは強さ2で考えることにする。

また、指数は増やすことはできても、ある一定の値以下にすることが不可能な場合があり、実質的には調整可能なパラメータにならないように思われる。例えば[3]で例示されている $L_8(2^7)$ (本稿の記法では $2\text{-}(2,7,2)$ ) については、因子数が多いのでどう足掻いても $\lambda =1$ にすることはできない。

Mutually orthogonal Latin squares (MOLS)

強さ2の直交表を構築するために、まずはmutually orthogonal Latin squares (MOLS) というものを構築する必要がある。以下ではこれについて説明する。

定義

MOLSを理解するためには、先にLatin squareとは何かを知っておく必要がある。Latin squareの定義を英語版Wikipedia[4]より引用する。

Latin square

In combinatorics and in experimental design, a Latin square is an $n \times n$ array filled with $n$ different symbols, each occurring exactly once in each row and exactly once in each column.

これは例を見たほうが早い。以下に次数3のLatin squareの例を示す。

1	2	3
2	3	1
3	1	2

これを見ると、任意の行・任意の列について、1～3の数字がちょうど一回ずつ表れていることが分かる。

Latin squareは同じ次数に対して複数存在する。その中である条件を満たしたものを集めた集合がMOLSである。MOLSの定義を英語版Wikipedia[5]より引用する。

Mutually orthogonal Latin squares

In combinatorics, two Latin squares of the same size (order) are said to be orthogonal if when superimposed the ordered paired entries in the positions are all distinct. A set of Latin squares, all of the same order, all pairs of which are orthogonal is called a set of mutually orthogonal Latin squares.

要するに互いに「直交する」Latin squareの集合がMOLSである。引用した定義の中に記載されているが、Latin squre同士が直交するというのは独特の用語であり、ベクトルの直交なんかとはあまり関連がないので注意されたい。

以下に次数4のMOLSの例を示す。これをMOLS(4) のように表記する。

1	2	3	4
2	1	4	3
3	4	1	2
4	3	2	1

1	2	3	4
4	3	2	1
2	1	4	3
3	4	1	2

1	2	3	4
3	4	1	2
4	3	2	1
2	1	4	3

これらのうち2つを選ぶと、それらは互いに「直交」している。例えば、1つ目と2つ目の表について、同じセル同士の数字で順序対を形成すると、各セルに対応する順序対は全て異なる。具体的には以下のようになる。

(1, 1)	(2, 2)	(3, 3)	(4, 4)
(2, 4)	(1, 3)	(4, 2)	(3, 1)
(3, 2)	(4, 1)	(1, 4)	(2, 3)
(4, 3)	(3, 4)	(2, 1)	(1, 2)

標準形

MOLSにはある種の自由度がある。例えば、属する全ての表について、行や列の入れ替えを行ってもやはりMOLSとなる。このような本質的でない自由度を排除するため、MOLSの標準形というものを考える[5]。

MOLSの標準形の定義を以下に示す。

Standard form

any set of MOLS can be put into standard form, meaning that the first row of every square is identical and normally put in some natural order, and one square has its first column also in this order

要するに、標準形とは全ての表の最初の行が自然な順序で並んでおり、かつある1つの表については最初の列もその順序で並んでいるようなMOLSのことである。先ほど示したMOLS(4)の例は標準形になっている。

MOLS( $n$ )の要素数

MOLS( $n$ )の標準形の要素数はどうなるだろうか？ $n = 1$ の場合は自明なので、 $n \ge 2$ として考える。このとき、(2, 1)成分の値を考えると、標準形であることから1にはなり得ない。そのため、2〜 $n$ の $n-1$ 通りの値が考えられる。

ここで、(2, 1)成分の値が $m (2, 3, \cdots , n)$ になるようなLatin squareが2つ以上あると仮定する。すると、それらのうちの2つの表から得られる(2, 1)成分の順序対は $(m, m)$ となる。しかし、これは $(1, m)$ 成分の順序対と一致しているため、MOLSであるということに矛盾する。

そのため仮定は誤りで、(2, 1)成分の値が $m$ になるようなLatin squareは高々1つ存在する。すなわち、MOLS( $n$ )の要素数は $n-1$ 以下となる。

MOLSの完全集合

MOLS( $n$ )の要素数がちょうど $n-1$ になるとき、それをMOLSの完全集合と呼ぶ。完全集合は $n$ が素数冪のときには存在することが知られている。存在性の証明は単純明快で、完全集合を具体的に構築する手法 (詳細は後述) が存在する。

ちなみに、一般の $n$ について完全集合が存在するかどうかは不明らしい[5]。ちょっと気になったことを調べていると簡単に未解決問題にぶち当たるので、数学というのは案外分かっていないことが多いのだなぁと思った。

MOLSの構築方法

準備が整ったので、いよいよMOLSの構築方法について説明する。本稿では次数が素数冪の場合にMOLSの完全集合を求める方法について説明する。内容は[5]を参考にした。なお、記述をシンプルにするため、これ以降は行・列の番号は0始まりとする。

$p$ を素数、 $k$ を1以上の整数、 $q = p^k$ とする。このとき、MOLS( $q$ )を求める。有限体 $\mathbb{F}_q$ の乗法群は巡回群であるため (例えば青雪江[6]の命題1.11.38に証明がある) 、ある生成元 $\gamma$ が存在して、 $\mathbb{F}_q$ の元は以下のように書ける。

$\displaystyle{ a_0=0, a_1=1, a_2=\gamma, a_3=\gamma^2, \cdots , a_{q-1}=\gamma^{q-2} }$

なお、 $\gamma^{q-1} = 1$ である。

ここで、MOLSの要素に番号をつけて $L_r (r = 1, 2, \cdots , q-1)$ と表す。また、 $L_r$ の $(i, j)$ 成分を $L_r(i, j)$ と表記する。このとき、 $L_r(i, j)$ を以下のようにしてやればMOLSの完全集合が得られる。

$\displaystyle{ L_r(i, j) = a_j + a_r a_i }$
*2

正しくMOLSが構築できることの証明

上で示した方法はあまりにもシンプルなのだが、本当に正しいのだろうか?残念ながら英語版Wikipedia[5]に証明がないので、自力で考えてみた。

MOLSの定義に従い、以下を示せばよい。

任意の $r$ について $L_r$ はLatin squareであること
相異なる $r, s$ を任意に固定したとき、任意の $i, j$ について $L_r(i, j)$ と $L_s(i, j)$ の各セルからなる順序対は互いに異なること

1点目について、まず行 $i$ を固定する。ある $j, j'$ について $a_j + a_r a_i = a_{j'} + a_r a_i$ であるとすると、両辺から $a_r a_i$ を引いて $a_j = a_{j'}$ となる。よって $j = j'$ となる他ない。

次に列 $j$ を固定する。ある $i, i'$ について $a_j + a_r a_i = a_j + a_r a_{i'}$ であるとすると、両辺から $a_j$ を引いて $a_r a_i = a_r a_{i'}$ となる。 $a_r \ne 0$ なのでさらに両辺を $a_r$ で割ると $a_i = a_{i'}$ となる。よって $i = i'$ となる他ない。

以上により、同一行・同一列の値はすべて異なることが分かったので、 $L_r$ はLatin squareである。

2点目について、任意の $i, j$ に対して順序対 $(a_j + a_r a_i, a_j + a_s a_i)$ がすべて異なっていればよい。そこで、ある $(i, j)$ , $(i', j')$ ( $i \ne i' \vee j \ne j'$ )について2つの順序対 $(a_j + a_r a_i, a_j + a_s a_i)$ および $(a_{j'} + a_r a_{i'}, a_{j'} + a_s a_{i'})$ が等しいと仮定して矛盾を導く。仮定より以下の2つの等式が成立する。

$\displaystyle{ \begin{eqnarray} a_j + a_r a_i &=& a_{j'} + a_r a_{i'} \\ a_j + a_s a_i &=& a_{j'} + a_s a_{i'} \end{eqnarray} }$

上の式から下の式を引いて整理すると以下のようになる。

$\displaystyle{ \begin{eqnarray} a_r a_i - a_s a_i = a_r a_{i'} - a_s a_{i'} \\ a_r a_i - a_r a_{i'} = a_s a_i - a_s a_{i'} \\ a_r(a_i - a_{i'}) = a_s (a_i - a_{i'}) \end{eqnarray} }$

もし $i = i'$ とすると、同じ行に同じ値が2つ現れることになるため、これは $L_r, L_s$ がLatin squareであることに反する。よって $i \ne i'$ であり、これより $a_i \ne a_{i'}$ が分かる。 $\mathbb{F}_q$ は体なので零元以外での除算を考えることができる。そこで、上式の両辺を $a_{i} - a_{i'}$ で割ると $a_r = a_s$ が得られる。これは $r \ne s$ に矛盾する。よって順序対 $(a_j + a_r a_i, a_j + a_s a_i)$ および $(a_{j'} + a_r a_{i'}, a_{j'} + a_s a_{i'})$ は異なる。

以上により、前述の方法で得られるテーブル群はMOLSになっていることが分かった。

数値への変換

前述の方法で得られたテーブルの要素は $\mathbb{F}_q$ の元となっている。このままだと使いにくいので、これを数値に変換しておくとよいだろう。これは簡単で、以下のような写像 $f$ によって $\mathbb{F}_q$ の元を自然数に対応付ければよい。

$\displaystyle{ \begin{eqnarray} f: a_i \mapsto i+1 \end{eqnarray} }$

こうすることで自然とMOLSの標準形が得られる。

MOLSを構築してみよう

MOLSの構築方法が分かったので、実際に構築してみよう。以下に本手法のsageによる実装例を示す。

#!/usr/bin/env sage

import sys
from sage.all import *

def isOrthogonal(ls1, ls2):
    orderedPairs = set()
    for row in range(len(ls1)):
        for col in range(len(ls1)):
            orderedPair = (ls1[row][col], ls2[row][col])
            if orderedPair in orderedPairs:
                return False
            orderedPairs.add(orderedPair)
    return True

def calcMOLS(q):
    gf = GF(q, name="a", modulus="primitive")
    gamma = gf.gen()
    orderedGf = [0] + [gamma**(val) for val in range(q-1)]

    mols = []
    for r in range(1, q):
        latinSquare = [[0 for i in range(q)] for j in range(q)]
        for row in range(q):
            for col in range(q):
                latinSquare[row][col] = orderedGf[col] + orderedGf[r] * orderedGf[row]
        mols.append(latinSquare)

    # Hasmap to convert elements of GF(q) to corresponding integers
    gfToInt = {}
    for i, val in enumerate(orderedGf):
        gfToInt[val] = i+1

    for i, latSq in enumerate(mols):
        for j, line in enumerate(latSq):
            mols[i][j] = [gfToInt[val] for val in line]

    return mols

def main(argc, argv):
    if argc != 2:
        print("usage: " + argv[0] + " order")
        return

    q = int(argv[1])

    print("=== Calc MOLS ===")
    mols = calcMOLS(q)
    for i, latSq in enumerate(mols):
        print("Latin square#{}:".format(i+1))
        for line in latSq:
            print(line)

    print("=== Sanity check ===")
    for r in range(len(mols)-1):
        for s in range(r+1, len(mols)):
            if not isOrthogonal(mols[r], mols[s]):
                print("The result is wrong.")
                return
    print("Correct MOLS.")

if __name__ == "__main__":
    main(len(sys.argv), sys.argv)

実行例を以下に示す。

shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./mols.py 2
=== Calc MOLS ===
Latin square#1:
[1, 2]
[2, 1]
=== Sanity check ===
Correct MOLS.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./mols.py 3
=== Calc MOLS ===
Latin square#1:
[1, 2, 3]
[2, 3, 1]
[3, 1, 2]
Latin square#2:
[1, 2, 3]
[3, 1, 2]
[2, 3, 1]
=== Sanity check ===
Correct MOLS.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./mols.py 4
=== Calc MOLS ===
Latin square#1:
[1, 2, 3, 4]
[2, 1, 4, 3]
[3, 4, 1, 2]
[4, 3, 2, 1]
Latin square#2:
[1, 2, 3, 4]
[3, 4, 1, 2]
[4, 3, 2, 1]
[2, 1, 4, 3]
Latin square#3:
[1, 2, 3, 4]
[4, 3, 2, 1]
[2, 1, 4, 3]
[3, 4, 1, 2]
=== Sanity check ===
Correct MOLS.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./mols.py 5
=== Calc MOLS ===
Latin square#1:
[1, 2, 3, 4, 5]
[2, 3, 5, 1, 4]
[3, 5, 4, 2, 1]
[4, 1, 2, 5, 3]
[5, 4, 1, 3, 2]
Latin square#2:
[1, 2, 3, 4, 5]
[3, 5, 4, 2, 1]
[4, 1, 2, 5, 3]
[5, 4, 1, 3, 2]
[2, 3, 5, 1, 4]
Latin square#3:
[1, 2, 3, 4, 5]
[4, 1, 2, 5, 3]
[5, 4, 1, 3, 2]
[2, 3, 5, 1, 4]
[3, 5, 4, 2, 1]
Latin square#4:
[1, 2, 3, 4, 5]
[5, 4, 1, 3, 2]
[2, 3, 5, 1, 4]
[3, 5, 4, 2, 1]
[4, 1, 2, 5, 3]
=== Sanity check ===
Correct MOLS.
shinya@DESKTOP-Q1NSEEU:~/programs/blog/2021% ./mols.py 6
=== Calc MOLS ===
/usr/lib/python3/dist-packages/apport/report.py:13: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses
  import fnmatch, glob, traceback, errno, sys, atexit, locale, imp, stat
Traceback (most recent call last):
  File "./mols.py", line 63, in <module>
    main(len(sys.argv), sys.argv)
  File "./mols.py", line 48, in main
    mols = calcMOLS(q)
  File "./mols.py", line 17, in calcMOLS
    gf = GF(q, name="a", modulus="primitive")
  File "sage/structure/factory.pyx", line 367, in sage.structure.factory.UniqueFactory.__call__ (build/cythonized/sage/structure/factory.c:2162)
  File "/usr/lib/python3/dist-packages/sage/rings/finite_rings/finite_field_constructor.py", line 566, in create_key_and_extra_args
    raise ValueError("the order of a finite field must be a prime power")
ValueError: the order of a finite field must be a prime power

コマンドライン引数で次数 $q$ を受け取っているが、確かにサイズ $q-1$ のMOLSが構築できた。次数6のケースで異常終了しているが、6は素数冪ではないためこれは期待通りである。

余談だが、sageは本当に万能で惚れ惚れする。今回は自分のPC (Windows10のWSL2で動作するUbuntu20.04) にインストールしたsageを使用したが、最近ではそんなことしなくてもブラウザ上で実行できたりするようである[7]。興味のある方は試してみるとよい。

まとめ

本稿では直交表の基本的な事項について説明した後、直交表の生成に必要となるMOLSの構築方法について述べた。具体的には有限体の性質をうまく使ってMOLSが構築できることを説明した。また、説明した手法を実際にsageで動かして正しさを実証した。

次回はMOLSを用いて直交表を生成する方法について書きたいと思う。恐らく今回ほどのボリュームにはならないだろう。さらにその次は有限射影幾何との関連についても考えてみたいが、こちらは現時点ではさっぱり見通しがない。頑張って勉強しようと思う。

参考

[1] 直交表 - Wikipedia
[2] Orthogonal array - Wikipedia
[3] www.ipros.jp
[4] Latin square - Wikipedia
[5] Mutually orthogonal Latin squares - Wikipedia
[6]

代数学2　環と体とガロア理論

作者:雪江　明彦
発売日: 2010/12/07
メディア: 単行本（ソフトカバー）

[7] doc.sagemath.org

*1:英語を読めばいいというのはその通りで、自分も今回そうやって勉強したが、やはり取っ掛かりとして母語の記事があるというのは理解のスピードが違ってくる。

*2:参考にした英語版Wikipedia[5]は $i, j$ がこの式と逆になっているが、それだと標準形にならない。

2021-04-01

その仮説検定、本当に大丈夫ですか？〜 2つの誤りをコントロールしよう

統計学

統計検定2級に合格しました

最近、統計検定2級に合格した。素直に嬉しい。2級は大学教養レベルの内容ではあるが、それらを良く理解していないと合格できないので、良い勉強になった。

余談だが、以下のサイトは2級合格に必要な情報が網羅されており、大変参考になった。

bellcurve.jp

統計検定2級では、とにかく区間推定や仮説検定の問題がたくさん出る。その勉強の中で、仮説検定について私が今までよく理解出来ていなかった点があることに気づいたため、本稿ではその気付きについて紹介したい。

棄却できない帰無仮説の問題

状況設定

ある正規母集団から $n$ 個のサンプルを抽出し、標本平均 $\overline{X}$ の値を元に母平均 $\mu$ の値について仮説検定をすることを考えてみよう。帰無仮説は「 $\mu = 0$ 」とし、対立仮説は「 $\mu \ne 0$ 」とする。これには両側検定を用いるのが適切である。

以後、有意水準はすべて5%固定で考える。また、議論をメインの話題にフォーカスするため、母分散 $\sigma^2$ は既知とする。つまり、t分布は使わない。さらに、簡単のため $\sigma = 1$ の場合のみを考えることにする*1。

もし帰無仮説が正しい場合、 $\overline{X}$ は正規分布 $N(0, \sigma^2/n)$ に従うはずである。そのため、帰無仮説の採択域は以下のようになる。

$\displaystyle{ - Z_{0.975} \frac{\sigma}{\sqrt{n}} < \overline{X} < Z_{0.025} \frac{\sigma}{\sqrt{n}} }$

ただし、 $Z_{\alpha}$ は標準正規分布 $N(0, 1)$ の上側 $100\alpha$ %点を表す。

もしサンプル数が少なく、例えば $n=3$ だった場合、採択域は以下のようになる。

$\displaystyle{ - 1.13 < \overline{X} < 1.13 }$

ここで、実はこの正規母集団の真の分布が $N(1, \sigma^2)$ だったとしよう。つまり、帰無仮説 $\mu = 0$ は誤りだったというわけである。このとき、母平均 $\mu = 1$ を中心に $\overline{X}$ が分布するわけだが、 $n=3$ のケースだと $\mu$ が帰無仮説の採択域の中に入ってしまっているので、これだと結構な確率で帰無仮説が棄却できないことになる。つまり、本当は棄却しなければならない帰無仮説を棄却できず、判断を誤ることになる。

何が問題か？

どうしてこのようなことになってしまうのだろうか？清く正しく仮説検定をやったはずなのに、どうして判断を誤ってしまうのだろうか？それはずばり、仮定した分布と真の分布の切り分けを付けるための分解能が足りなかったことが原因である。

仮定した母平均に対して真の母平均がある程度近くにある場合、サンプル数が足りないと $\overline{X}$ の分散が大きくなってしまい、2つの確率分布は大きくオーバーラップする。その結果、帰無仮説が容易には棄却できないということになってしまうのである。サンプル数が多くなれば $\overline{X}$ の分散が小さくなっていくため、仮定した分布と真の分布がはっきりと別れるようになり、仮説検定により区別が付くようになる。

第一種の誤りと第二種の誤り

実は、このような概念は統計学の世界においては普通に知られていることである。それが仮説検定における「第一種の誤り」と「第二種の誤り」というものである。それぞれ以下のような意味で使われる。

第一種の誤り：帰無仮説が正しいのに誤って棄却すること
第二種の誤り：帰無仮説が間違っているのに誤って棄却しないこと

今回のケースで言うと、サンプル数が少ない場合、第二種の誤りを起こす確率が高くなってしまうということである。以下では第一種・第二種の誤りを起こす確率をそれぞれ $\alpha, \beta$ と表記する。

仮説検定では、第一種の誤りの許容範囲については、有意水準という形で明示的に指定する。今回は有意水準5%としているので、5%の確率で第一種の誤りが起こることは許容していることになる。それに対して、第二種の誤りについては明示的に述べられることがあまりないように思う。しかし、第二種の誤りの存在も認識して仮説検定のパラメータを選ばないと、仮説検定から有益な結果が得られなくなるので、どちらも意識してやる必要がある。

2つの誤りの関係

第一種・第二種の誤りをどちらも小さくすることが理想であるが、両者は一般的にトレードオフの関係にある。なぜなら、第一種の誤りを減らすということは「間違って帰無仮説を棄却したくない」ということになるのに対して、第二種の誤りを減らすということは「間違って帰無仮説を採択したくない」ということに等しく、両者が相反するからである。

では、全く打つ手が無いのかというと、そうではない。大切なのは、まず第一種の誤りと第二種の誤りの存在を正しく認識し、それらがどういうパラメータに影響を受けるかを知った上で、自分の考えている問題設定において、2つの誤りが許容できる範囲内に収まるようにコントロールすることである。これにより、仮説検定の有益性を格段に向上することができる。

2つの誤りの可視化

ここまで文章だけで説明してしまったので、第一種の誤りと第二種の誤りを可視化してみよう。可視化にはDesmosというサービスを利用した。状況設定として、仮説検定で $N(0, 1)$ を仮定したが、真の分布は $N(1, 1)$ だったというケースについてグラフを描いた。

www.desmos.com

赤線が仮定した分布、青線が真の分布である。また、赤色部分の面積が $\alpha$ 、青色部分の面積が $\beta$ の値を表している。 $\alpha$ は固定で5%としている。これを見ると、サンプル数を増やすことでグラフがシャープになって2つの分布のオーバーラップが小さくなり、 $\beta$ が減っていく様子が見て取れる。

俺の考えた最強の仮説検定

以上を踏まえて、正しい仮説検定のステップについて考えてみたいと思う。

それぞれの誤り確率を左右する要因を押さえる

第一種・第二種の誤りの発生確率を左右する要因として代表的なものを以下にまとめる。

要因	第一種の誤りへの影響	第二種の誤りへの影響
有意水準 $\alpha$	↑	↓
サンプル数 $n$	影響なし	↓
母分散 $\sigma^2$	影響なし	↑

表の見方について補足する。要因の列に示した値が増加したときに、第一種・第二種の誤りが増えるか減るかを上下の矢印で示した。例えば、有意水準 $\alpha$ が大きくなれば、許容する第一種の誤り確率が大きくなるため、第一種の誤りは増加する。

許容される誤り確率を設定しよう

次に、第一種・第二種の誤りがどれくらいの確率で発生することを許容できるかを考えよう。一般的には $\alpha = 0.05, \beta = 0.2$ くらいが目安のようだが、これはもちろん仮説検定毎に異なる。例えば、第一種の誤りを起こすことが人命に関わるようなケースでは、第一種の誤りを極力起こしたくないというようなことがあるだろう。

真の分布を推定しよう

これが仮説検定の泣き所である。実は、当然だが第二種の誤りは真の分布が分からないと分からない。しかし、そもそも真の分布が分からないから仮説検定をやっているわけなので、これは矛盾している。しかし、仮説検定は残念ながらこの矛盾を常に孕んだ状態で行わなければならない宿命にある。

この問題に対する私の考えを書いておく。真の分布は究極的には分からないので、何らかの推定をするしかない。これにはいろいろなやり方があるはずである。

例えば、製造したビールの容量が本当に500mlと言えるかを仮説検定するような場合、過去に誤って480mlだったことがあるというデータがあれば、平均が480mlの正規分布を真の分布として $\beta$ を求める手はあるだろう。

他には、例えば学習塾のとある講座を受けることで試験の点が上がったかどうかを仮説検定する場合について、保護者の皆様が5点以上の点数向上を求めているという期待される効果量があったとする。

このケースでは、「実は講座の効果として真に点数が2点向上する効果があった」というような場合、第二種の誤りを犯しても犯さなくても、保護者の期待を裏切っていることには変わりない。そのため、2点差を切り分けられる分解能を得るためにサンプルをたくさん用意することには、コストがかかるだけで意味がない。

一方、本当は点数が5点向上する効果があったのに第二種の誤りを犯すということは避けたい。そのため、受講者の平均点が5点向上した場合を元に真の分布を仮定し、 $\beta$ を算出するのが良いだろう。

このように、過去のデータ、求められている効果、及び専門家の知見などから、真の分布を仮定してやるのがよいと考えられる。

パラメータを決定しよう

最後に、先程決めた $\alpha, \beta$ の値を元に、仮説検定のための各種パラメータを設定する。母平均の推定の例で言えば、コントロールできるパラメータは有意水準 $\alpha$ とサンプル数 $n$ になる。 $\alpha$ は前のステップで既に決めてあるので、あとは求める $\beta$ の条件を満たすように $n$ を決定すればよい。例えば、先程の状況設定において $\alpha = 0.05, \beta = 0.2$ であれば、 $n \ge 8$ とすればよい。詳細は先程張ったDesmosのリンク先にて確かめてみて欲しい。

まとめ

本稿では仮説検定における第一種の誤りと第二種の誤りについて述べた。仮説検定では第一種の誤りだけでなく第二種の誤りも存在することを正しく認識し、それらをどのような値に納めたいかを検討した上で各種パラメータを決定する必要があるということを説明した。

統計学は非常に面白いので統計検定1級まで取りたい気持ちもありつつ、近頃は人生の時間の有限性に思いを馳せる日々である。自分と家族の幸せのために何をすべきか、よく考えていきたい。

参考

[1]

統計学入門 (基礎統計学Ⅰ)

発売日: 1991/07/09
メディア: 単行本

*1:このように状況を限定しなくても同様の議論ができるのだが、面倒なので簡単な状況に限定した。

2021-03-04

クーポンコレクター問題の拡張とその確率分布

確率論

以前、クーポンコレクター問題の確率分布について調べたことがある (詳細はこちら) 。これはこれで大きな成果だったのだが、いざこれを実用しようと思うと、少々物足りないことに気づいた。

例えば、6種類のクーポンを全て集めることを考える。以前求めた確率分布では、クーポンを集め始める前の時点における確率を考えることはできる。しかし、実際にクーポンが揃っていく速さはその時々であり、例えば途中で早々に2種類集まった場合、そこからどれくらいの確率で残り4つのクーポンが集まるのかを知りたくなる。

つまり、知りたいのは $k$ 種類のクーポンのうちすでに $l$ 種類のクーポンを取得している状況で、残りの $k-l$ 種類のクーポンを集めるために必要な試行回数の確率分布はどうなるのか？ということである。本稿ではこの疑問に対する回答について述べる。

考え方の方針

$n$ 回の試行で出現し得るクーポンの出方の総数は $k^n$ となる。そのため、すでに $l$ 種類のクーポンを持っている状態から始めて、ちょうど $n$ 回で全てのクーポンが揃うパターンの総数が分かれば確率が分かる。

考えを整理するために、いくつか場合分けをしてパターンを数え上げてみよう。以下ではクーポンは全部で $k$ 種類あり、そのうちすでに $l$ 種類のクーポンを持っているとする。

$n < k-l$ のとき

残りのクーポンの数より小さい試行回数でクーポンが集まることはあり得ないので、このケースは0通りである。

$n = k-l$ のとき

それぞれの試行で未取得のクーポンが毎回ダブリなく出る必要がある。このケースは残りのクーポンの並び替えの総数と等しく、 $(k-l)!$ 通りとなる。

$k-l < n \le k$ のとき

ここが山場である。このパターンでちょうど $n$ 回目にクーポンが全て集まるためには、 $n-1$ 回目までに $k-1$ 個のクーポンが出尽くしている必要がある。試行回数が残りのクーポンの数よりも多いため、最初の時点で持っているクーポンが出ることもあり得る。それも加味すると、 $n-1$ 回目までに出るクーポンは $k-l-1$ 種類以上、 $n-1$ 種類以下となる。

最初に持っている $l$ 種類のクーポンのうち $i$ 種類のクーポンが出るとすると、 $n-1$ 回目までには $k-l-1+i$ 種類のクーポンが出ることになる。

それぞれの試行に番号を付けて、これを $k-l-1+i$ 個の区別された箱に分類しているのだと考えると、このパターンの総数は以下のように第二種スターリング数を用いて書ける。

$\displaystyle{ \left(\textstyle{k-l \atop k-l-1}\right) \left(\textstyle{l \atop i}\right) (k-l-1+i)! \{\textstyle{n-1 \atop k-l-1+i}\} }$

第二種スターリング数自体は区別のない箱への分類の総数を表しているので、 $(k-l-1+i)!$ を掛けていることに注意されたい。また、すでに持っているクーポンから $i$ 種類のクーポンを選ぶ方法には自由度がある。同様に、未取得のクーポンから $k-l-1$ 種類のクーポンを選ぶ方法にも自由度がある。そのため、さらに $\left(\textstyle{k-l \atop k-l-1}\right) \left(\textstyle{l \atop i}\right)$ を掛けている。

$0 \le i \le n-k+l$ より、 $i$ による総和を取れば、求めるパターンの総数は以下のようになる。

$\displaystyle{ \sum_{i=0}^{n-k+l} \left(\textstyle{k-l \atop k-l-1}\right) \left(\textstyle{l \atop i}\right) (k-l-1+i)! \{\textstyle{n-1 \atop k-l-1+i}\} }$

クーポンの出方の例: $k=6, l=3, n=5$

この場合、 $0 \le i \le 2$ である。それぞれの $i$ についてクーポンの出方の具体例を見てみよう。

ここでは6種類のクーポンをAからFのアルファベットで表す。また、最初の時点ですでにA, B, Cの3つのクーポンを持っているとする。

・ $i=0$

すでに持っているクーポンは1種類も出ないため、以下のようなクーポンの出方となる。

D F F F E

・ $i=1$

すでに持っているクーポンは1種類だけ出るため、以下のようなクーポンの出方となる。

D A F A E

・ $i=2$

すでに持っているクーポンは2種類だけ出るため、以下のようなクーポンの出方となる。

B F A D E

$k < n$ のとき

これは1つ前のケースとほぼ同様である。ただし、 $n-1$ 回目までに出るクーポンは $k-l-1$ 種類以上、 $k-1$ 種類以下となる。ゆえに $0 \le i \le l$ となる。これより、このパターンの総数は以下のようになる。

$\displaystyle{ \sum_{i=0}^{l} \left(\textstyle{k-l \atop k-l-1}\right) \left(\textstyle{l \atop i}\right) (k-l-1+i)! \{\textstyle{n-1 \atop k-l-1+i}\} }$

クーポンの出方の例: $k=6, l=3, n=8$

この場合、 $0 \le i \le 3$ である。それぞれの $i$ についてクーポンの出方の具体例を見てみよう。

ここでも6種類のクーポンをAからFのアルファベットで表し、最初からA, B, Cの3つのクーポンを持っているとする。

・ $i=0$

すでに持っているクーポンは1種類も出ないため、以下のようなクーポンの出方となる。

D F F F D D F E

・ $i=1$

すでに持っているクーポンは1種類だけ出るため、以下のようなクーポンの出方となる。

D A F A A F D E

・ $i=2$

すでに持っているクーポンは2種類だけ出るため、以下のようなクーポンの出方となる。

B F F A D A B E

・ $i=3$

すでに持っているクーポンは3種類全て出るため、以下のようなクーポンの出方となる。

B F F C D A B E

4つのケースをまとめる

実は前述の4つのケースは以下のように一本の式にまとめられる。

$\displaystyle{ \sum_{i=0}^{\mathrm{min}(n-k+l, l)} \left(\textstyle{k-l \atop k-l-1}\right) \left(\textstyle{l \atop i}\right) (k-l-1+i)! \{\textstyle{n-1 \atop k-l-1+i}\} }$

これを $k^n$ で割れば求める確率分布が得られる。

$\displaystyle{ P(N = n) = \sum_{i=0}^{\mathrm{min}(n-k+l, l)} \frac{\left(\textstyle{k-l \atop k-l-1}\right) \left(\textstyle{l \atop i}\right) (k-l-1+i)! \{\textstyle{n-1 \atop k-l-1+i}\}}{k^n} }$

ただし、 $N$ はクーポンが全て集まるまでに必要な試行回数を表す確率変数であるとする。

計算してみよう

確率分布

前回同様、今回もpythonを使用して確率分布を計算してグラフを描いてみた。パラメータが多いので $k$ の値は6固定とした。ソースコードを以下に示す。

#!/usr/bin/python3.6
from sympy.functions.combinatorial.numbers import stirling
from sympy.functions.combinatorial.numbers import nC
import matplotlib.pyplot as plt
import math

def coupon(k, l, n):
    loopNum = min(n-k+l, l)+1

    numPattern = 0
    for i in range(loopNum):
        stir = stirling(n-1, k-l-1+i, kind=2)
        numPattern += (k-l) * nC(l, i) * math.factorial(k-l-1+i) * stir
    return float(numPattern) / k**n

def main():
    numCouponKind = 6
    numTrial = 30
    for l in range(0, numCouponKind):
        x = []
        y = []
        for n in range(1, numTrial+1):
            prob = coupon(numCouponKind, l, n)
            x.append(n)
            y.append(prob)
        plt.plot(x, y, label='l = {}'.format(l))
    plt.legend()
    plt.grid()
    plt.xlabel("n")
    plt.ylabel("probability")
    plt.savefig('coupon_lx_k6.png')
    #plt.show()

if __name__ == "__main__":
    main()

得られたグラフを以下に示す。

f:id:peng225:20210116010347p:plain — 6種類のクーポンのうちすでにl種類のクーポンを持っている場合の確率分布

$l$ が大きいほど分布が左に偏っており (つまり右に歪んでおり) 、最初に持っているクーポンが多ければ多いほど少ない試行回数でクーポンを全て収集できることが分かる。

累積分布

これまでの議論で確率分布は分かった。しかし、実際に自分がどれくらい試行を重ねればクーポンを揃えられるかというのは、累積分布を見た方が分かりやすいだろう。

私が日常生活で良く考えたくなるのは $k=3, 6$ のケースなので*1、これらについて累積分布関数の値がそれぞれ0.5, 0.9, 0.95以上になるために必要な試行回数を $l$ の値毎に示す。

・ $k=3$

l	累積0.5	累積0.9	累積0.95
0	5	9	11
1	4	8	10
2	2	6	8

・ $k=6$

l	累積0.5	累積0.9	累積0.95
0	13	23	27
1	12	22	26
2	11	21	24
3	10	19	23
4	7	17	21
5	4	13	17

クーポンを確実に集めるためには、思ったより多くの試行が必要という印象を受けた。クーポン集めは最後の方が異常にしんどいので、 $l$ が多少大きくても焼け石に水という感じだということも分かった。

まとめ

本稿ではクーポンコレクター問題の拡張として、最初からいくつかのクーポンを持っている状態でクーポンを全て集める際の試行回数について考えた。クーポンを確実に集めるために必要な試行回数は思いのほか大きく、クーポンを真面目に集めることがいかに厳しいことであるかを肌で感じた。

今後の人生で何かを集めたくなったら、大人しくメル○リを利用した方が良いかもしれない。

*1:子どものためにハッピーセットのおもちゃを集めることがある。

有限射影幾何を用いた直交表構築手法の証明

おさらい

証明

計算してみよう

ちなみに

まとめ

参考

行列を用いた直交表の構築手法

手法の概要

疑問点のオンパレードと解決の筋道

有限体上の射影空間

射影空間

有限体の場合

有限射影空間の完全代表系

有限射影空間の存在条件

有限射影空間上の点の数

有限射影空間と行列

行列の求め方

強さ2の場合

強さ3以上の場合

計算してみよう

次元の決め方

まとめ

参考

直交表の構築

基本の構築方法

直交表が得られていることの証明

準備

0列目と1列目

0列目と2列目以降

1列目と2列目以降

2列目以降同士

MOLSの限界とその他の手法

因子数の減らし方

MOLSの限界

その他の直交表構築法

直交表を構築してみよう

まとめ

参考

直交表とは

定義

直交表を決めるパラメータ

Mutually orthogonal Latin squares (MOLS)

定義

標準形

MOLS()の要素数

MOLSの完全集合

MOLSの構築方法

正しくMOLSが構築できることの証明

数値への変換

MOLSを構築してみよう

まとめ

参考

統計検定2級に合格しました

棄却できない帰無仮説の問題

状況設定

何が問題か？

第一種の誤りと第二種の誤り

2つの誤りの関係

2つの誤りの可視化

俺の考えた最強の仮説検定

それぞれの誤り確率を左右する要因を押さえる

許容される誤り確率を設定しよう

真の分布を推定しよう

パラメータを決定しよう

まとめ

参考

考え方の方針

のとき

のとき

のとき

クーポンの出方の例:

のとき

クーポンの出方の例:

4つのケースをまとめる

計算してみよう

確率分布

累積分布

まとめ

$G$ 行列を用いた直交表の構築手法

有限射影空間と $G$ 行列

$G$ 行列の求め方

次元 $d$ の決め方

MOLS( $n$ )の要素数

$n < k-l$ のとき

$n = k-l$ のとき

$k-l < n \le k$ のとき

クーポンの出方の例: $k=6, l=3, n=5$

$k < n$ のとき

クーポンの出方の例: $k=6, l=3, n=8$