今年に入ってから「酔歩」をテーマとするシリーズに取り組んできたのですが，ｄ次元格子上の数え上げプログラムを設計しているうちに，それが「適合度検定」にも応用できるに違いないというアイデアがひらめきました．これまで適合度検定の「Fisher's exact probability」の計算プログラムには何度かチャレンジしながら挫折を繰り返していただけに，思わぬところで光明をみたような気がしました．

　

　酔歩は物理数学上の問題，かたや統計検定上の問題であって，ほとんど接点がないかのように思えるのですが，「多項分布」に従うことが両者のほとんど唯一の共通点になっています．今回のコラムでは，酔歩問題の番外編として，フィッシャーの直接確率の計算法を取り上げますが，酔歩問題の予期せざる効用といってよいでしょう．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【１】χ2検定による適合度の近似計算

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

　適合度検定では，母集団がｋ個のカテゴリーに分類されていて，各カテゴリーの母比率が

　　ｐ1，ｐ2，・・・，ｐk　　　（ｐ1＋ｐ2＋・・・＋ｐk＝１）

であることを検定する問題を考えます．ｎ個の無作為標本を抽出したとき，各カテゴリーの観測度数（Ｏ）と期待度数（Ｅ）が，以下のような１×ｋ分割表の形式であったとしましょう．

　

　　カテゴリー　　Ｃ1　　　Ｃ2　　・・・　Ｃk　　　計

　　観測度数　　　ｎ1　　　ｎ2　　・・・　ｎk　　　ｎ

　　期待度数　　　ｎｐ1　　ｎｐ2　・・・　ｎｐk　　ｎ

　

　適合度検定ではχ2検定あるいは尤度比検定が用いられることが多いわけですが，χ2検定については，１９００年，イギリスの統計学者カール・ピアソンが実測度数と期待度数の差の２乗を期待度数で割った量の和がカイ２乗分布に従うことを示したことに始まります．

　

　ピアソンの適合度検定では，理論値と実測値の一致の度合いのよしあしをみるために，各階級ごとに

　　（Ｏ−Ｅ）^2 ／Ｅ　　　（Ｏ：観測度数　Ｅ：期待度数）

を加えて

　　χ2 ＝Σ（Ｏ−Ｅ）^2 ／Ｅ

を作って，このχ2 統計量が自由度＝階級数−１のカイ２乗分布をすることを用いて理論分布と実測値の適合度を検定します．ピアソンのχ2 検定は，分布全体にわたって適合していればχ2 統計量が小さくなるという素朴で見やすい尺度になっています．

　

　ピアソンの適合度検定の欠点は，

ａ）ｎが大きいとモデルが棄却されやすい（適合しないと判定される）．

ｂ）ｎが小さいとモデルが棄却されにくい（適合すると判定される）．

ｃ）χ2 値の自由度がｎに依存せず，もっぱら階級数によって決まる．→このとき，区間数や区間の分割点をどのように定めたらよいかという問題が生ずる．

ｄ）漸近分布がχ2 分布でよく近似されるためには各セルの期待度数≧５でなければ用いることができない（近似的にχ2 分布に従うことの証明においてスターリングの公式を使用しているため）．→５以下の項目は適当に合併する．→合併する場合，自由度が減る．

ｅ）理論度数は小さいのに観測度数が大きい少数の特異点が全体の当てはめに大きく影響を与える．

などがあげられます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【２】適合度検定の exact probability

　

　χ2 検定法はよく用いられますが，あくまでも近似計算であって，期待度数Ｅが小さいとき不正確であるなど理論的にいろいろな問題点のあることが指摘されています．そこで，χ2検定や尤度比検定には頼らずに，正確なｐ値（exact probability）を算出することを考えてみます．

　

　最初に，２項分布となるケース：メンデルの遺伝学におけるエンドー豆の例を引いて説明しますが，

　

　　　　　　観測数　　期待数

　　優性　　5474　　　7342*3/4=5493

　　劣性　　1850　　　7324*1/4=1831

　　　計　　7324　　　7324

　

　優性：劣性が３：１になるという仮説では，現在の度数よりもさらに偏ったものが起こる確率を計算します．したがって，正確なｐ値は

　exact probability＝ΣnＣk(3/4)^k(1/4)^(n-k)

のように２項分布の和の形で求められることになります．

　

　２項分布の場合は簡単に計算ができて，実際に計算すると上の例では

　　EXACT PROBABILITY= .608168

となり，χ2検定や尤度比検定の近似ｐ値（それぞれP= .614516，P= .614293）にほぼ一致しました．検定結果を掲げておきます．

　

(1) DATA SUMMARY

Category　　　No.　　　　　 Prob.

1　　　　　　 5474　　　　　.75

2　　　　　　 1850　　　　　.25

(2) PEARSON CHI-SQUARE TEST

CHI^2= .26288　　( N.S )　 ( P= .614516 )

(3) modified CHI-SQUARE TEST

modified CHI^2= .261083　　( N.S )　 ( P= .615708 )

(4) LIKELIHOOD RATIO TEST

LIKELIHOOD RATIO= .263218　　( N.S )　 ( P= .614293 )

(5) EXACT PROBABILITY= .608168

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　次に，２項分布モデルを拡張して，多項分布モデルを扱ってみることにします．２種類の形質の表現型の頻度が９：３：３：１となる例が入手できなかったので，以下の血液型のデータを用いてみます．

　

　　　　　　　観測数　　期待数

　　Ａ型　　　62　　　　200*4/10=80

　　Ｏ型　　　84　　　　200*3/10=60

　　Ｂ型　　　30　　　　200*2/10=40

　　ＡＢ型　　24　　　　200*1/10=20

　　　計　　　200　　　 200

　

　日本人の血液型の分布はＡ：Ｏ：Ｂ：ＡＢ＝４：３：２：１といわれていますが，ある２００人の集団の血液型を調べたところ，６２：８４：３０：２４であったという例です．Ｏ型の割合がとくに大きいようですが，日本人としては特異な集団といえるかどうかを検定してみましょう．

　

(1) DATA SUMMARY

Category　　　No.　　　　　 Prob.

1　　　　　　 62　　　　　　.4

2　　　　　　 84　　　　　　.3

3　　　　　　 30　　　　　　.2

4　　　　　　 24　　　　　　.1

(2) PEARSON CHI-SQUARE TEST

CHI^2= 16.95　　( P<0.01 )　 ( P= .0010516 )

(3) modified CHI-SQUARE TEST

modified CHI^2= 16.0829　　( P<0.01 )　 ( P= .0014559 )

(4) LIKELIHOOD RATIO TEST

LIKELIHOOD RATIO= 16.4112　　( P<0.01 )　 ( P= 1.28573E-03 )

　検定結果より，この集団は日本人の平均的な血液型分布とは異なるといえます．

　

　さて，観測度数ｎ1，ｎ2，ｎ3，ｎ4が出現する確率は，４項分布（一般には多項分布）

　　n!/n1!n2!n3!n4!(4/10)^n1(3/10)^n2(2/10)^n3(1/10)^n4

で与えられるわけですが，正確なｐ値を求めるには２項分布の場合と同様に，現在の度数よりもさらに偏ったものが起こる確率を計算します．

　

　ここで，カテゴリー数がｋの場合のexact probabilityを求めるプロシージャをまとめておきます．

　

（１）カテゴリー　　Ｃ1　　　Ｃ2　　・・・　Ｃk　　　計

　　　観測度数　　　ｎ1　　　ｎ2　　・・・　ｎk　　　ｎ

　　　期待度数　　　ｎｐ1　　ｎｐ2　・・・　ｎｐk　　ｎ

のi=1,2,･･･,kに対して，与えられた事象が起こる確率

　　n!/n1!n2!･･･!nk!(p1)^n1(p2)^n2･･･(pk)^nk

を計算する．

（２）ｎ1＋ｎ2＋・・・＋ｎk＝ｎの条件のもとで，多項分布のすべてのパターンの確率を計算し，現在の度数よりもさらに偏ったものが起こったとき，その確率を足し合わせる．

（３）ｐ値を帰無仮説の有意水準αと比較して，それ以下であったら帰無仮説を棄却する．

　

　２項分布の場合の数え上げプログラムは容易に作ることができます．しかし，３項分布になったとたん，この問題は突如として難しくなります．２項分布の場合と違って，多項分布ではまず過不足のない正確な数え上げプログラムを設計すること自体が大変になるのです．

　

　また，ｎが少し大きくなると，すべてのパターンの確率を足し合わせて正確なｐ値を求めるのに大変な時間がかかってしまいます．そのため，探索時間ができるだけ短くて済むような工夫が必要になりますが，これが簡単ではないのです．このような問題はとかくとり漏らしやすいもので，探索範囲を狭めようとすると見逃されてしまうものがあるやもしれないという不安と恐れ，そしてストレスがつきまといます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【３】フィッシャー直接法の計算プログラム

　

1000 '

1010 ' ** exact test **

1020 '

1030 DEFDBL A-G,L-Z

1040 '

1050 M=2

1060 DIM N(M),O(M)

1070 DIM P(M),Q(M)

1080 N(1)=5474:P(1)=.75

1090 N(2)=1850:P(2)=.25

1100 NN=7324:DIM G(NN)

1110 GOSUB *EXACT.P

1120 END

1130 '

1140 M=3

1150 DIM N(M),O(M)

1160 DIM P(M),Q(M)

1170 N(1)=5 :P(1)=.2

1180 N(2)=1 :P(2)=.3

1190 N(3)=4 :P(3)=.5

1200 NN=10:DIM G(NN)

1210 GOSUB *EXACT.P

1220 END

1230 '

1240 M=4

1250 DIM N(M),O(M)

1260 DIM P(M),Q(M)

1270 N(1)=8 :P(1)=.25

1280 N(2)=7 :P(2)=.25

1290 N(3)=5 :P(3)=.25

1300 N(4)=0 :P(4)=.25

1310 NN=20:DIM G(NN)

1320 GOSUB *EXACT.P

1330 END

　

2000 '

2010 ' ** exact test **

2020 '

2030 *EXACT.P:

2040 IF M>4 THEN RETURN

2050 IF NN^(M-1)>10000 THEN RETURN

2060 '

2070 UNIFORM=1

2080 FOR I=1 TO M-1

2090　 IF P(I)<>P(I+1) THEN UNIFORM=0

2100 NEXT I

2110 '

2120 GGG=0:G(0)=0

2130 FOR I=1 TO NN

2140　 GGG=GGG+LOG(I)

2150　 G(I)=GGG

2160 NEXT I

2170 '

2180 S0=G(NN)

2190 FOR I=1 TO M

2200　 S0=S0-G(N(I))+N(I)*LOG(P(I))

2210 NEXT I

2220 '

2230 IF M=2 THEN GOSUB *NOMIAL2

2240 IF M=3 THEN GOSUB *NOMIAL3

2250 IF M=4 THEN GOSUB *NOMIAL4

2260 PRINT SS

2270 RETURN

2280 '

2290 ' ** ２項分布 **

2300 '

2310 *NOMIAL2:

2320 SS=0

2330 MAXI=NN

2340 MINI=-INT(-NN/2)

2350 FOR I=MAXI TO MINI STEP -1

2360　 J=NN-I

2370　 GOSUB *TIE2

2380　 GOSUB *PROB2

2390 NEXT I

2400 RETURN

2410 '

2420 *TIE2:

2430 C=2:F=2

2440 IF I=J THEN C=1

2450 RETURN

2460 '

2470 *PROB2:

2480 O(1)=I:O(2)=J

2490 IF UNIFORM=1 THEN GOSUB *UNIFORM:GOTO 2530

2500 '

2510 Q(1)=P(1):Q(2)=P(2):GOSUB *PROB.SUB

2520 Q(1)=P(2):Q(2)=P(1):GOSUB *PROB.SUB

2530 RETURN

2540 '

2550 ' ** ３項分布 **

2560 '

2570 *NOMIAL3:

2580 SS=0

2590 MAXI=NN

2600 MINI=-INT(-NN/3)

2610 FOR I=MAXI TO MINI STEP -1

2620　 IF (NN-I)<I THEN MAXJ=NN-I ELSE MAXJ=I

2630　 MINJ=-INT(-(NN-I)/2)

2640　 FOR J=MAXJ TO MINJ STEP -1

2650　　 K=NN-I-J

2660　　 GOSUB *TIE3

2670　　 GOSUB *PROB3

2680　 NEXT J

2690 NEXT I

2700 RETURN

2710 '

2720 *TIE3:

2730 C=6:F=6

2740 IF I=J AND J=K THEN C=1:GOTO 2760

2750 IF I=J OR J=K THEN C=3

2760 RETURN

2770 '

2780 *PROB3:

2790 O(1)=I:O(2)=J:O(3)=K

2800 IF UNIFORM=1 THEN GOSUB *UNIFORM:GOTO 2880

2810 '

2820 Q(1)=P(1):Q(2)=P(2):Q(3)=P(3):GOSUB *PROB.SUB

2830 Q(1)=P(1):Q(2)=P(3):Q(3)=P(2):GOSUB *PROB.SUB

2840 Q(1)=P(2):Q(2)=P(1):Q(3)=P(3):GOSUB *PROB.SUB

2850 Q(1)=P(2):Q(2)=P(3):Q(3)=P(1):GOSUB *PROB.SUB

2860 Q(1)=P(3):Q(2)=P(1):Q(3)=P(2):GOSUB *PROB.SUB

2870 Q(1)=P(3):Q(2)=P(2):Q(3)=P(1):GOSUB *PROB.SUB

2880 RETURN

2890 '

2900 ' ** ４項分布 **

2910 '

2920 *NOMIAL4:

2930 SS=0

2940 MAXI=NN

2950 MINI=-INT(-NN/4)

2960 FOR I=MAXI TO MINI STEP -1

2970　 IF (NN-I)<I THEN MAXJ=NN-I ELSE MAXJ=I

2980　 MINJ=-INT(-(NN-I)/3)

2990　 FOR J=MAXJ TO MINJ STEP -1

3000　　 IF (NN-I-J)<J THEN MAXK=NN-I-J ELSE MAXK=J

3010　　 MINK=-INT(-(NN-I-J)/2)

3020　　 FOR K=MAXK TO MINK STEP -1

3030　　　 L=NN-I-J-K

3040　　　 GOSUB *TIE4

3050　　　 GOSUB *PROB4

3060　　 NEXT K

3070　 NEXT J

3080 NEXT I

3090 RETURN

3100 '

3110 *TIE4:

3120 C=24:F=24

3130 IF I=J AND J=K AND K=L THEN C=1:GOTO 3200

3140 IF I=J AND J=K THEN C=4　　　　:GOTO 3200

3150 IF J=K AND K=L THEN C=4　　　　:GOTO 3200

3160 IF I=J AND K=L THEN C=6　　　　:GOTO 3200

3170 IF I=J THEN C=12　　　　　　　 :GOTO 3200

3180 IF J=K THEN C=12　　　　　　　 :GOTO 3200

3190 IF K=L THEN C=12

3200 RETURN

3210 '

3220 *PROB4:

3230 O(1)=I:O(2)=J:O(3)=K:O(4)=L

3240 IF UNIFORM=1 THEN GOSUB *UNIFORM:GOTO 3500

3250 '

3260 Q(1)=P(1):Q(2)=P(2):Q(3)=P(3):Q(4)=P(4):GOSUB *PROB.SUB

3270 Q(1)=P(1):Q(2)=P(2):Q(3)=P(4):Q(4)=P(3):GOSUB *PROB.SUB

3280 Q(1)=P(1):Q(2)=P(3):Q(3)=P(2):Q(4)=P(4):GOSUB *PROB.SUB

3290 Q(1)=P(1):Q(2)=P(3):Q(3)=P(4):Q(4)=P(2):GOSUB *PROB.SUB

3300 Q(1)=P(1):Q(2)=P(4):Q(3)=P(2):Q(4)=P(3):GOSUB *PROB.SUB

3310 Q(1)=P(1):Q(2)=P(4):Q(3)=P(3):Q(4)=P(2):GOSUB *PROB.SUB

3320 Q(1)=P(2):Q(2)=P(1):Q(3)=P(3):Q(4)=P(4):GOSUB *PROB.SUB

3330 Q(1)=P(2):Q(2)=P(1):Q(3)=P(4):Q(4)=P(3):GOSUB *PROB.SUB

3340 Q(1)=P(2):Q(2)=P(3):Q(3)=P(1):Q(4)=P(4):GOSUB *PROB.SUB

3350 Q(1)=P(2):Q(2)=P(3):Q(3)=P(4):Q(4)=P(1):GOSUB *PROB.SUB

3360 Q(1)=P(2):Q(2)=P(4):Q(3)=P(1):Q(4)=P(3):GOSUB *PROB.SUB

3370 Q(1)=P(2):Q(2)=P(4):Q(3)=P(3):Q(4)=P(1):GOSUB *PROB.SUB

3380 Q(1)=P(3):Q(2)=P(1):Q(3)=P(2):Q(4)=P(4):GOSUB *PROB.SUB

3390 Q(1)=P(3):Q(2)=P(1):Q(3)=P(4):Q(4)=P(2):GOSUB *PROB.SUB

3400 Q(1)=P(3):Q(2)=P(2):Q(3)=P(1):Q(4)=P(4):GOSUB *PROB.SUB

3410 Q(1)=P(3):Q(2)=P(2):Q(3)=P(4):Q(4)=P(1):GOSUB *PROB.SUB

3420 Q(1)=P(3):Q(2)=P(4):Q(3)=P(1):Q(4)=P(2):GOSUB *PROB.SUB

3430 Q(1)=P(3):Q(2)=P(4):Q(3)=P(2):Q(4)=P(1):GOSUB *PROB.SUB

3440 Q(1)=P(4):Q(2)=P(1):Q(3)=P(2):Q(4)=P(3):GOSUB *PROB.SUB

3450 Q(1)=P(4):Q(2)=P(1):Q(3)=P(3):Q(4)=P(2):GOSUB *PROB.SUB

3460 Q(1)=P(4):Q(2)=P(2):Q(3)=P(1):Q(4)=P(3):GOSUB *PROB.SUB

3470 Q(1)=P(4):Q(2)=P(2):Q(3)=P(3):Q(4)=P(1):GOSUB *PROB.SUB

3480 Q(1)=P(4):Q(2)=P(3):Q(3)=P(1):Q(4)=P(2):GOSUB *PROB.SUB

3490 Q(1)=P(4):Q(2)=P(3):Q(3)=P(2):Q(4)=P(1):GOSUB *PROB.SUB

3500 RETURN

3510 '

3520 *PROB.SUB:

3530 LP=G(NN)

3540 FOR H=1 TO M

3550　 LP=LP-G(O(H))+O(H)*LOG(Q(H))

3560 NEXT H

3570 IF LP>S0 THEN RETURN

3580 SS=SS+EXP(LP)*C/F

3590 RETURN

3600 '

3610 *UNIFORM:

3620 LP=G(NN)+NN*LOG(P(1))

3630 FOR H=1 TO M

3640　 LP=LP-G(O(H))

3650 NEXT H

3660 IF LP>S0 THEN RETURN

3670 SS=SS+EXP(LP)*C

3680 RETURN

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　３項分布の場合で説明しますが，このプログラムで最も工夫した点は２６１０行，２６４０行のループです．ｉ＋ｊ＋ｋ＝ｎとなる（ｉ，ｊ，ｋ）の組を順次探索しながら累積確率を求めているのですが，探索に時間がかかるため，対称性を利用して探索範囲を絞り込み，ｎ≧ｉ≧ｊ≧ｋ≧０となるような最小範囲の探索で済むようにしています．

　

　これによって，計算量を約１／６に減らすことができますから，その分，あとで因子Ｃを掛けることにします．同順位がない場合，Ｃにはその組合せ数

　　Ｃ＝３！／１！１！１！＝６

を割り当てます．

　一方，同順位の数があった場合，たとえば，ワンペアでは

　　Ｃ＝３！／２！１！＝３

スリーカードでは，

　　Ｃ＝３！／３！＝１

となり，Ｃにはそれぞれ３と１が割り付けられます．

　

　また，プログラムの２８２０〜２８７０行では順列（辞書式）を生成していますが，どうみてもスマートとはいえません．当初はスマートな順列生成プログラムを組むことを考えていたのですが，どうやったところで効率は大変悪そうなので，結局，泥臭く計算することにしました．プログラムがダサくなってしまいますが，計算時間を実用的な範囲内に抑えるためには，泥臭さは避けられないと考えたのです．→【補】順列生成プログラム

　

　なお，２８２０〜２８７０行の確率を辞書式に並べることによって，実行結果が昇順になることが期待され，そうすれば計算時間を短くする工夫ができるのですが，そうはならなかったことを付記しておきます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【４】実行例

　

　血液型の例題についての計算では，EXACT PROBABILITY= 0.00089605となりました．また，

　　N(1)=5 :P(1)=.2

　　N(2)=1 :P(2)=.3

　　N(3)=4 :P(3)=.5

の例では，p=0.069981

　

　　N(1)=8 :P(1)=.25

　　N(2)=7 :P(2)=.25

　　N(3)=5 :P(3)=.25

　　N(4)=0 :P(4)=.25

の例では，p=0.019660と計算されて，それぞれの正確値p=0.0700,p=0.020に一致する結果が得られました．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【５】あとがき

　

　ここに掲げたプログラムでは，ｎ：標本数，ｋ：カテゴリー数として，

　　ｋ＞４，あるいは，ｎ^(k-1)＞１００００

のときには，計算を実行しないようにしてあります．計算に長時間を要し，実用的ではないと思われるからです．そのような制限はあるのですが，ともかく，長い間，未解決であった問題にケリがついてほっとしています．

　

　次に考えるとすれば，ｒ×ｃ分割表に対するフィッシャーの直接法プログラムということになるでしょう．以前，２×２分割表に対するフィッシャーの直接法プログラムを作ったことがあるのですが，たった２×２分割表であっても標本数が大きければ多くの計算所要時間を要します．実際にプログラムを作ってみるとわかるのですが，周辺和を固定して得られる全ての多項分布の数は急激に増加しますから，フィッシャーの直接確率（exact probability）を計算することがいかに大変な労力を要するのか，おわかり頂けるでしょう．

　

　今回のコラムを書くにあたって，フィッシャーの直接法を２×２分割表より大きい一般のｒ×ｃ分割表へ拡張する方法について調べてみたところ，高速ネットワークアルゴリズム：

　　Mehta CR, Patel NR(1983): A network algorithm for performing Fisher's exact test in rxc contingency tables,

J American Statistical Association 78, 427-437

があることがわかりました．

　

　ただし，ｎ：標本数，ｒ：行数，ｃ：列数とすると，

　　ｎ／（ｒ−１）（ｃ−１）＞５

であれば，一般的なｒ×ｃ分割表に対しての検定には長時間を要し，実用的なものではないとのことでした．それでも，このアルゴリズムは多くの問題点をクリアしているわけですから大いに興味あるところです．早速，文献を取り寄せて調べてみたいと考えています．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】順列生成プログラム

　

　たとえば（1,2,3,4）の順列とは，(2,1,3,4)とか(4,3,2,1)とか，全部で４！＝２４組あるわけですが，それらをすべて作り出すにはどのようにすればよいのでしょうか？

　

　４つの数であればまだよいのですが，ｎ個の数字の順列ならば，ｎを与えてｎ重のfor-nextループを作ることが基本かと思われますが，それには再帰的手続きが必要となります．しかし，BASICでは再帰ができませんし，たとえ再帰的処理が可能な言語であってもスタックを大量に消費して，実行時間はかなり遅くなるでしょう．結局，私はいい解決法を思いつきませんでした．

　

　そこで，知人の知恵を拝借したのですが，Ｓ氏の談によると

　(1)Mathematicaでは順列を列挙する関数があるが，実行時間は決して良いとはいえない．たとえば，十個の順列

　　Permutations[{a, b, c, d, e, f, g, h, i, j}]

をMathematicaで求めようとすると，結局，仮想メモリを使い尽くし答えを返してはくれない．また，11個の順列の個数はマシンの整数の上限をこえてしまう．ほかの言語の場合，もっとひどいことになることが予想される．

　

　(2)swapコマンドでは，２数の交換しかできない．しかるにｎ＞２の場合の順列を作りたい．そのような場合，Knuth博士の「The Art of Computer Programing」（とくに第三巻）に，互換と順列に関連するいくつかの定理が提示してある．第三巻の記述では，二項の交替から順列を生成することに関する議論があるが，実用に耐える理論かどうかは判断がつかない．

ということでした．

　

　別のアプローチをすべきであろうと思われるのですが，実用的な計算時間を考慮してスマートな順列生成を断念し，泥臭く計算することにしたのです．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝