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

　N88-BASICの精度は特に指定しない限り７桁（単精度）ですが，倍精度指定すると１６桁の精度で計算することができます．通常のデータならば単精度で十分正確な値が得られますが，精度を考慮するならば倍精度で計算するのが安全です．

　倍精度の指定には定数や変数の後ろに＃印を付けるか，DEFDBL A-Zのように倍精度型宣言します．しかし，変数を倍精度形式にすると計算速度は遅くなりますし，FOR～NEXTループの制御変数（ループカウンタ）には整数型か単精度型で用いなければなりませんから，倍精度にするとループの管理がおかしくなってしまいます．すなわち，高精度計算にとってFOR～NEXTループは厄介者になっているのです．

　一方，WHILE～WENDは高精度計算することも含んで作られた命令ですから，今回のコラムでは超幾何関数を用いた確率分布の計算プログラム（その７）のFOR～NEXTループをWHILE～WENDループに変換して，超幾何関数の連分数展開式の精度について再検討してみることにしました．

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

　FOR～NEXTは繰り返し回数がわかっているときに用いる反復命令ですが，繰り返し回数がまえもって決まっていないときは，WHILE～WENDを使います．WENDはEnd of Whileを意味するBASICの予約語です．しかし，WHILE～WENDの効用はそればかりではありません．

１）条件が成立したときにループを終了する場合

　FOR～NEXTのなかで，IF～THEN GOTO～文を働かせ，条件によってループを抜け出すことが平気で行われていますが，この書き方では弊害をもたらすことがよくあります．WHILE～WENDの処理機構には，FOR～NEXTとIF～THEN GOTO～を同時に処理する機能がありますから，WHILE～WENDではGOTO文を使わないで，スマートに仕上げることができます．

２）高精度計算する場合

　FOR～NEXTは最もとっつきやすい反復命令なのですが，前述したようにループカウンタを倍精度にするとループの管理がおかしくなって高精度計算できません．したがって，精密計算の効果を得るには，WHILE～WENDの方が望ましいのです．以下では，実例をみてみましょう．

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

【超幾何関数の連分数展開】

　以下に掲げるプログラムは，合流型超幾何関数1F1(1,c,x)，ガウス型超幾何関数2F1(a,1,c,x)の連分数展開式です．ベキ級数展開では，収束の様子を探りながら繰り返しを続けるか，打ち切るかを決めていましたが，収束状況の判定時間が惜しいので，連分数展開ではあらかじめ何回繰り返せばよいのかを決めておき（変数名：NZ），FOR～NEXTループを用いて，無条件にその回数だけ繰り返しています．

（FOR～NEXTループ）

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

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

17520 FOR ID=NZ TO 1 STEP -1

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

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

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

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

17580 G#=-Z/C/(1#+G#)

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

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

17670 FOR ID=NZ TO 1 STEP -1

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

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

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

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

　しかし，このプログラムではすべての変数を倍精度形式（DEFDBL A-Z）にするとループカウンタ（変数名：ID）の管理がおかしくなってしまいますから，その問題を解消するためにWHILE～WENDループに変換してみることにします．

（WHILE～WENDループ）

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

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

17520 'FOR ID=NZ TO 1 STEP -1

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

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

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

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

17580 G#=-Z/C/(1#+G#)

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

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

17670 'FOR ID=NZ TO 1 STEP -1

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

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

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

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

　このプログラムでは，連分数展開式は％点の小さいところで精度がよいのですが，大きいところでは，ベキ級数展開よりもかなり悪くなってしまいました．計算項数を多くすれば（NZ:20→50），精度が向上するのではないかと考えられましたが，時間がかかるうえに誤差も小さくはなりませんでした．今回の改良によって倍精度型宣言が可能になりましたから，変数をすべて倍精度形式（DEFDBL A-Z）にして計算し直したところ，ベキ級数展開の場合と遜色のない精度が得られました．

　しかし，一般的にいって，連分数展開の精度はベキ級数展開に及ばず，精度がよいとする従来の見解には懐疑的にならざるを得ません．理論上の収束性がよいというのと計算機での計算の容易性は相関関係がないと思った方が良さそうです．なお，超幾何関数の連分数展開式については，一松信「特殊関数入門」森北出版を参考にしました．

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

【対数ガンマ関数】

　ガンマ関数Γ（ｘ）は引数ｘが大きくなるとすぐにオーバーフローエラーを起こしてしまいますので，コンピュータを使うには危険です．たとえば，３５とか３６とかという引数では飽和現象を起こします．そこで，対数ガンマ関数ｌｎ（Γ（ｘ））を求めることによって飽和現象を回避します．

　以下の対数ガンマ関数のサブルーチンは，統計数値表（日本規格協会）に掲載されていたFORTRANプログラム（ハートの近似式）をそっくりそのままBASICに書き直したものです．その際，FORTRANのDO～I=1,7はBASICのFOR I=1 TO 7～NEXTに変換されます．

1010 ' *** 対数ガンマ関数 ***

1040 A1#= .918938533204673#

1050 B1#= .000766345188#

1060 B2#=-.00059409561052#

1070 B3#= .0007936431104845#

1080 B4#=-.00277777775657725#

1090 B5#= .08333333333331692#

1100 IF (Z-8)<0 THEN X=Z+8:FLG=-1:GOTO 1120

1130 A#=(X-.5#)*LOG(X)-X+A1#

1140 B#=((((B1#*Y+B2#)*Y+B3#)*Y+B4#)*Y+B5#)/X

1160 IF FLG>0 THEN RETURN

1210　　A#=A#*(X-FI)

1230 L.GAMMA=L.GAMMA-LOG(A#)

　ここで気になるのは，やはりFOR～NEXTループです．WHILE～WENDループに書き直してみましょう．

1300　 A#=A#*(X-FI)

1330 L.GAMMA=L.GAMMA-LOG(A#)

　また，対数ガンマ関数であることを意識して，次のように書き直すこともできます．このほうがスマートかもしれません．

1300　 A#=A#+LOG(X-FI)

1330 L.GAMMA=L.GAMMA-A#

　この近似式は単精度の範囲で数値表の値とよく一致しているのですが，ヘイスティングスやコリンジの近似式よりも精度は劣っているようです．

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