超幾何関数を用いた確率分布の計算（その４）

　「理数嫌い」や「理系離れ」が叫ばれている昨今，小生のみならず，この状態を非常に危険なものと憂いているひとは少なくないだろう．教育の家元たる文部省が詰め込み教育に歯止めをかけようとして「ゆとり教育」を標榜するあまり，分数のできない大学生が大量発生し，さらにπ＝３が国家規格となるようでは文部省のもくろみは完全な失敗だったと断言できよう．

　ゆとり教育は過当競争社会を是正するために考え出されたものなのであろうが，教育にはアメとムチの両方が必要で，アメをしゃぶらせておくだけでムチで打つことを忘れてしまったツケが回ってきたのである．このままでは巻き返しはますます困難になるだろうから，小手先だけの対応でなく，抜本的な発想の転換が迫られてしかるべきであり，再度，国家百年の計という視点からの基本理念の変革が必要と思われるが，文部省主導の教育体制（家元制度）に対して，正直いって小生はあまり期待していない．文部省の失策は，市中引き回しのうえ打ち首・獄門に価する．にもかかわらず，一向にその責任を全うしようとはしないし，逆に罪を現場になすりつけているフシさえ見受けられるからである．

　現場を任された教師達が主導すべきと思うが，それは教師達に一層の奮励努力を強いることになり，酷な要望ともいえるだろう．しかし，現場には「文部省のお役人や教育庁の取りまき連中は何もわかっていない」とする気骨ある教師も少なくないにちがいない．酷と知りつつも，教育のプロを自認する教師達にエールを送りたいのである．

　さて，ホームページは科学教育を行う媒体としても有用なものと思われるが，小生は数年前より科学の面白さを伝えることを目的として，ＨＰに科学記事を連載する活動を続けている．現在のところ，このＨＰには１万件を超えるアクセスがあるが，この数は私的に開設したホームページとしては破格の値であろう．このＨＰはＵＲＬ

　　http://www.mcc.pref.miyagi.jp/people/ikuro/

で参照可能であり，今回のコラムで記念すべき第５０話に至るが，いまとなってはＨＰを利用した科学教育の実践は小生の日々の研究課題といってもよい．このような些細な活動が「理数嫌い」や「理系離れ」に歯止めをかける後方支援として役立ってくれれば幸甚なのであるが，・・・．

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

　当初，ＨＰの今月のコラム（閑話休題）は，１話平均原稿用紙２０枚程度の記事でスタートしたのだが，次第に冗長になり，最近では１００枚を越すことも珍しくはなくなった．これは大いに反省すべき点である．所期の目的に立ち戻るためにも，第５０話では要点をかい摘んで述べ，できるだけ短報にとどめたい．

　非線形方程式ｆ（ｘ）＝０の解をニュートン法で計算する場合，たとえば，２次までの微分を用いたニュートン法では

　　ｆ（ｘ），ｆ’（ｘ），ｆ”（ｘ）

の値を求めることが必要になる．このことから，超幾何関数を用いた確率分布の計算を正攻法で解こうとすると，最低でも３回は超幾何関数値の計算ルーチンにアクセスして，正規分布・χ2分布では合流型超幾何関数

　　1F1(a,c,x),1F1(a+1,c+1,x),1F1(a+2,c+2,x)

ｔ分布・Ｆ分布ではガウス型超幾何関数

　　2F1(a,b,c,x),2F1(a+1,b+1,c+1,x),2F1(a+2,b+2,c+2,x)

の値を計算することが要求される．

　前回掲げたプログラムは，超幾何関数のパラメータに関する漸化式を駆使して，アクセス回数を２回に抑えたものであった．たとえば，正規分布のサブルーチン

1170 ' ** NORMAL DISTRIBUTION **

1190 *NORMAL.PROB.HGF:

1200 PI#=3.141592653589732#

1230 Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF:G1#=G#

1240 P0=.5#+G1#*UUN*EXP(-UUN*UUN/2)/SQR(PI#*2)

1270 Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF:G2#=G#*(A0/C0-1)+G1#

1280 P10=G1#*(1-UUN*UUN)+G2#*UUN*UUN

1290 P1=P10*EXP(-UUN*UUN/2)/SQR(PI#*2)

1310 G3#=G1#*A0/Z-G2#*(C0+Z)/Z

1320 P20=-P10*UUN+G2#*UUN*(1-UUN*UUN)-G1#*UUN*2+G3#*UUN*UUN*UUN+G2#*UUN*2

1330 P2=P20*EXP(-UUN*UUN/2)/SQR(PI#*2)

では，１２３０行と１２７０行で２回のアクセスを行っているが，今回の改良プログラムではもう一歩からめ手から攻めて，いずれの分布においてもアクセス回数を１回に減らすことに成功した．

（改良プログラム）

1010 ' ** NORMAL DISTRIBUTION **

1030 *NORMAL.PROB.HGF:

1040 PI#=3.141592653589732#

1070 Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF :G1#=G#

1080 P0=.5#+G1#*UUN*EXP(-UUN*UUN/2)/SQR(PI#*2)

1110 A=A0-1:C=C0:G0#=1#

1120 G#=(G1#-G0#)*C0/Z

1130 G2#=G#*(A0/C0-1)+G1#

1140 P10=G1#*(1-UUN*UUN)+G2#*UUN*UUN

1150 P1=P10*EXP(-UUN*UUN/2)/SQR(PI#*2)

1170 G3#=G1#*A0/Z-G2#*(C0+Z)/Z

1180 P20=-P10*UUN+G2#*UUN*(1-UUN*UUN)-G1#*UUN*2+G3#*UUN*UUN*UUN+G2#*UUN*2

1190 P2=P20*EXP(-UUN*UUN/2)/SQR(PI#*2)

1220 ' ** CHI-SQUARE DISTRIBUTION **

1240 *CHI2.PROB.HGF:

1270 Z=UUX/2:GOSUB *CONFLUENT.HGF :G1#=G#

1280 Z=1+DF/2:GOSUB *LOG.GAMMA:G2#=G#

1290　P0=G1#/EXP(G2#)*(UUX/2)^(DF/2)*EXP(-UUX/2)

1320 A=A0-1:C=C0:G0#=1#

1330 G#=(G1#-G0#)*C0/Z

1340 G3#=G#*(A0/C0-1)+G1#

1350 P10=DF/2*(UUX/2)^(DF/2-1)*G1#/2-(UUX/2)^(DF/2)*G1#/2+(UUX/2)^(DF/2)*G3#/2

1360　P1=P10/EXP(G2#)*EXP(-UUX/2)

1380 G4#=G1#*A0/Z-G3#*(C0+Z)/Z

1390 P20=-P10/2+(DF/2*(DF/2-1)*(UUX/2)^(DF/2-2)-DF/2*(UUX/2)^(DF/2-1))*G1#/4+(DF/2*(UUX/2)^(DF/2-1)-(UUX/2)^(DF/2))*G3#/4+DF/2*(UUX/2)^(DF/2-1)*G3#/4+(UUX/2)^(DF/2)*G4#/4

1400 P2=P20/EXP(G2#)*EXP(-UUX/2)

1430 ' ** T-DISTRIBUTION **

1460 DFX=DF/2:DFY=.5

1470 A0=DFX+DFY:B0=1:C0=DFX+1

1490 Z=DF/(DF+UUT*UUT):GOSUB *GAUSS.HGF :G1#=G#

1500 P0=G1#/DFX*Z^DFX*(1-Z)^DFY

1510 Z=DFX+DFY:GOSUB *LOG.GAMMA:G2#=G#

1520 Z=DFX　　:GOSUB *LOG.GAMMA:G3#=G#

1530 Z=DFY　　:GOSUB *LOG.GAMMA:G4#=G#

1540　P0=P0*EXP(G2#-G3#-G4#)

1580 Z=DF/(DF+UUT*UUT)

1590 A=A0+1:B=B0-1:C=C0:G0#=1#

1600 G#=-(G1#*(C0-A0-1)+G0#*(B0-C0))/(A0-B0+1)/(1-Z)

1610 G5#=(G#-G1#)*A0/Z

1620 ZZ=-DF*UUT*2/(DF+UUT*UUT)^2

1630 P10=(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G1#+Z^DFX*(1-Z)^DFY*G5#

1640　P1=-P10*ZZ/DFX*EXP(G2#-G3#-G4#)/2

1660 A=A0+1:B=B0-1:C=C0:G0#=1#

1670 G6#=-G#*A0*(C0-B0*2+(B0-A0-1)*Z)/Z/(1-Z)-G0#*A0*(B0-C0)/Z/(1-Z)-G5#*C0/Z

1680 ZZZ=-DF*2*(DF-UUT*UUT*3)/(DF+UUT*UUT)^3

1690 P20=P10*ZZZ+((DFX*(DFX-1)*Z^(DFX-2)*(1-Z)^DFY-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)+Z^DFX*DFY*(DFY-1)*(1-Z)^(DFY-2))*G1#+(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G5#+Z^DFX*(1-Z)^DFY*G6#)*ZZ*ZZ

1700 P2=-P20/DFX*EXP(G2#-G3#-G4#)/2

1730 ' ** F-DISTRIBUTION **

1760 DFX=DF2/2:DFY=DF1/2

1770 A0=DFX+DFY:B0=1:C0=DFX+1

1790 Z=DF2/(DF2+DF1*UUF):GOSUB *GAUSS.HGF :G1#=G#

1800 P0=G1#/DFX*Z^DFX*(1-Z)^DFY

1810 Z=DFX+DFY:GOSUB *LOG.GAMMA:G2#=G#

1820 Z=DFX　　:GOSUB *LOG.GAMMA:G3#=G#

1830 Z=DFY　　:GOSUB *LOG.GAMMA:G4#=G#

1840　P0=P0*EXP(G2#-G3#-G4#)

1870 Z=DF2/(DF2+DF1*UUF)

1880 A=A0+1:B=B0-1:C=C0:G0#=1#

1890 G#=-(G1#*(C0-A0-1)+G0#*(B0-C0))/(A0-B0+1)/(1-Z)

1900 G5#=(G#-G1#)*A0/Z

1910 ZZ=-DF2*DF1/(DF2+DF1*UUF)^2

1920 P10=(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G1#+Z^DFX*(1-Z)^DFY*G5#

1930　P1=-P10*ZZ/DFX*EXP(G2#-G3#-G4#)

1950 A=A0+1:B=B0-1:C=C0:G0#=1#

1960 G6#=-G#*A0*(C0-B0*2+(B0-A0-1)*Z)/Z/(1-Z)-G0#*A0*(B0-C0)/Z/(1-Z)-G5#*C0/Z

1970 ZZZ=-DF2*DF1*DF1*2/(DF2+DF1*UUF)^3

1980 P20=P10*ZZZ+((DFX*(DFX-1)*Z^(DFX-2)*(1-Z)^DFY-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)+Z^DFX*DFY*(DFY-1)*(1-Z)^(DFY-2))*G1#+(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G5#+Z^DFX*(1-Z)^DFY*G6#)*ZZ*ZZ

1990 P2=-P20/DFX*EXP(G2#-G3#-G4#)

　時間のかかる計算を極力回避することによって，確率分布の計算はさらに高速化されたが，そのトリックは，超幾何関数がそれぞれ

　　1F1(1,c,x),2F1(a,1,c,x)

の形をしていることである．説明は割愛するが，この形はパラメータに関する漸化式にとって，都合のよい形だとだけ付記しておく．

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

　実は，特殊形

　　1F1(1,c,x),2F1(a,1,c,x)

は連分数展開にとっても都合のよい形となっている．πの分数近似：22/7,355/113,･･･などは連分数から求められたものであるらしいが，最良の分数近似と聞いていたので，小生は連分数に対して過大な期待をもっていた．ところが，超幾何関数に対しては，ベキ級数版

2690 ' *** 合流型超幾何関数 ***

2710 *CONFLUENT.HGF:' [Kummer]

2750 WHILE ABS(G0#-G#)>EPSL

2760　 T#=T#*A/C*Z/K

2860 ' *** ガウス型超幾何関数 ***

2920 WHILE ABS(G0#-G#)>EPSL

2930　 T#=T#*A*B/C*Z/K

の方が，連分数版

2110 ' *** 合流型超幾何関数の連分数展開 ***

2130 *CONFLUENT.HGF1:' 1F1(1,c,x)

2160 FOR ID=NZ TO 1 STEP -1

2170　 EVEN#=ID/(C+2*ID-1)/(C+2*ID-2)

2180　 ODD#=-(C+ID-1)/(C+2*ID)/(C+2*ID-1)

2190　 G#=ODD#*Z/(1+G#)

2200　 G#=EVEN#*Z/(1+G#)

2260 ' *** ガウス型超幾何関数の連分数展開 ***

2280 *GAUSS.HGF1:' 2F1(a,1,c,x)

2310 FOR ID=NZ TO 1 STEP -1

2320　 ODD#=-(A+ID-1)*(C+ID-2)/(C+2*ID-2)/(C+2*ID-3)

2330　 EVEN#=-ID*(C-A+ID-1)/(C+2*ID-1)/(C+2*ID-2)

2340　 G#=EVEN#*Z/(1+G#)

2350　 G#=ODD#*Z/(1+G#)

より，よい結果を与えることがわかった．連分数展開式では項数（変数名：NZ）を大きくとっても誤差が小さくなるとは限らない．どうやら，逐次近似の場合，級数の形式の方がよいらしい．

　このことに関して，畏友・阪本氏より，

　　ａ）計算機の場合，表現可能な整数の範囲が限られているので，分子，分母がともに大きくなるのは困る．

　　ｂ）理論上の収束性がよいというのと，計算機での計算の容易性は相関関係がないと思った方が良い．

　　ｃ）同じ計算機でも，IBMなどの大型機は１６進数の浮動小数点，その他（パソコンなど）は２の二進なので，非常にやっかいだ．

　　ｄ）Mathematicaでは，どのように大きい数でも表現できるが，特殊関数などの組み込み関数は従来型の浮動小数点をつかっているようだ．

などのコメントを頂いた．連分数展開に対して，氏はハナから懐疑的であったようである．

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

　大急ぎで説明してきたが，逆にいえば，本稿のようなラフな説明でわかる人は大したものである．結論として，正規分布・χ2分布・ｔ分布・Ｆ分布は，いずれも超幾何関数で表すことができるが，この一連のコラムで，超幾何関数を用いた確率計算の有用性が実証されたことになる．

　一方，非心χ2分布・非心ｔ分布・非心Ｆ分布は，直接的に超幾何関数で表すことはできないが，それぞれχ2分布・ｔ分布・Ｆ分布の加重和として間接的には表せるので，ここで述べた方法は非心χ2分布・ｔ分布・Ｆ分布に対する計算においても有用であるはずである．有用性が判明され次第，ＨＰにアップロードしてみたい．

　とはいったものの，この種の数値計算プログラムを作る工程では，微分積分などの工程に多くの労力が注がれるので，その部分の検算だけでもってヘロヘロになってしまう．小生，突発性難聴に罹り，低音部の聴力低下，高音部はエコーがかかったように二度聞こえる，持続性の耳鳴りという症状に冒された経験をもっている．このような精神的緊張と肉体的負担を考慮に入れると，数値計算プログラムを作成するのに，もっとも安価な方法はMathematicaなどの数式処理ソフトウェアを用いることであろう．

　数式処理ソフトであれば，同じ結果を得るのに少ない行数でバグも少なくて済むから，それに越したことはない．多大な時間と労力を費やし，デバッグに腐心し，プログラミング言語と悪戦苦闘するのはもはや時代錯誤であって，決して他人に勧められることではない．これから数値計算を始めようとする研究者・技術者には，Mathematicaなどの数式処理ソフトをお勧めする次第である．

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