非心χ2分布の分布関数を

　　Ｆ（ｘ，ν，λ）　　ν：自由度　　λ：非心パラメータ

とすると，（その５）に掲げたプログラムは，λが既知のもとで，

　　１）パーセント点ｘ　→　下側確率Ｆ

　　２）下側確率Ｆ　→　パーセント点ｘ

を求めるものであった．

　

　今回（その６）では攻略法を変えて，パーセント点ｘと下側確率Ｆが与えられたときに，非心パラメータλを求めるプログラムを紹介する．

　　３）ｘ，Ｆ　→　非心パラメータλ

　

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

　

【プログラム】

　

10000 '

10010 ' ** NON-CENTRAL CHI-SQUARE DISTRIBUTION **

10020 '

10030 *NCHI2.PERCENT.HGF.LAMBDA:

10040 PP=1-P0:UU=UUX:GOSUB *NCHI2.LAMBDA

10050 PROB$="LX":GOSUB *NEWTON.HGF.LAMBDA

10060 LAMBDA=X

10070 RETURN

10080 '

10090 ' ** NON-CENTRAL T-DISTRIBUTION **

10100 '

10110 *NT.PERCENT.HGF.LAMBDA:

10120 PP=1-P0:UU=UUT:GOSUB *NT.LAMBDA

10130 PROB$="LT":GOSUB *NEWTON.HGF.LAMBDA

10140 LAMBDA=X

10150 RETURN

10160 '

10170 ' ** NON-CENTRAL F-DISTRIBUTION **

10180 '

10190 *NF.PERCENT.HGF.LAMBDA:

10200 PP=1-P0:UU=UUF:GOSUB *NF.LAMBDA

10210 PROB$="LF":GOSUB *NEWTON.HGF.LAMBDA

10220 LAMBDA=X

10230 RETURN

10240 '

10250 ' *** ニュートン法 ***

10260 '

10270 *NEWTON.HGF.LAMBDA:

10280 EPSL=1E-09

10290 IMAX=50

10300 '

10310 Q0=P0

10320 D0=LAMBDA

10330 '

10340 SW=0

10350 ICNT=0

10360 WHILE SW=0

10370　 ICNT=ICNT+1:PRINT "*";

10380　 IF ICNT>IMAX THEN SW=1

10390　 DD=D0:GOSUB 10430:D0=DD

10400 WEND

10410 X=DD

10420 RETURN

10430 '

10440 ' ** 逐次計算 **

10450 '

10460 IF PROB$="LX" THEN LAMBDA=D0:GOSUB *NCHI2.PROB.HGF.LAMBDA

10470 IF PROB$="LT" THEN LAMBDA=D0:GOSUB *NT.PROB.HGF.LAMBDA

10480 IF PROB$="LF" THEN LAMBDA=D0:GOSUB *NF.PROB.HGF.LAMBDA

10490 '

10500 F=L0-Q0:G=L1

10510 '

10520 F2=F:D2=D0

10530 IF F1=0 THEN 10560

10540　GOSUB *BISECTION.LAMBDA

10550 'GOSUB *DAMPING.LAMBDA

10560 '

10570 H=L2:DI=G*G-F*H*2

10580 IF DI<0 THEN DELTA=-F/G:GOTO 10610

10590 IF G<0 THEN DELTA=F*2/(-G+SQR(DI)) ELSE DELTA=F*2/(-G-SQR(DI))

10600 '

10610 DD=D0+DELTA

10620 '

10630 EPSL=.00001

10640　IF ABS(DD-D0)<EPSL THEN SW=1

10650 'IF ABS((DD-D0)/D0)<EPSL THEN SW=1

10660 F1=F2:D1=D2

10670 RETURN

10680 '

10690 *BISECTION.LAMBDA:

10700 IF SGN(F1)=SGN(F2) THEN RETURN

10710　 DD=(D1+D2)/2:' bisection

10720　 DD=(F1*D2-F2*D1)/(F1-F2):' regula falsi

10730 RETURN 10620

10740 '

10750 *DAMPING.LAMBDA:

10760 IF ABS(F1)>ABS(F2) THEN RETURN

10770　 DD=D0-DELTA/2:' damping

10780 RETURN 10620

10790 '

10800 ' ** NON-CENTRAL CHI-SQUARE DISTRIBUTION **

10810 '

10820 *NCHI2.PROB.HGF.LAMBDA:

10830 GOSUB *CHI2.PROB.HGF.SUB

10840 GOSUB *POISSON.MIX.CHI2.LAMBDA

10850 RETURN

10860 '

10870 ' *** ポアソン混合 ***

10880 '

10890 *POISSON.MIX.CHI2.LAMBDA:

10900 EPSL=1E-09

10910 Z=UUX/2:P0#=P0:H1#=G1#

10920 L0#= EXP(-LAMBDA/2)*P0#

10930 L1#=-EXP(-LAMBDA/2)*P0#/2

10940 L2#= EXP(-LAMBDA/2)*P0#/4

10950 SS#=0:S0#=L0#:S1#=L1#:S2#=L2#

10960 J=1

10970 WHILE ABS((SS#-S0#)/S0#)>EPSL

10980　 A=A0:C=C0+J-1

10990　 H2#=C/Z*(H1#-1#)

11000　 '

11010　 O#=H2#:I#=H1#

11020　 L0#=L0#*LAMBDA/2/J/C*O#/I#*Z

11030　 '

11040　 O#=(-.5+J/LAMBDA)*H2#

11050　 I#=(-.5+(J-1)/LAMBDA)*H1#

11060　 L1#=L1#*LAMBDA/2/J/C*O#/I#*Z

11070　 '

11080　 O#=((-.5+J/LAMBDA)^2-J/LAMBDA/LAMBDA)*H2#

11090　 I#=((-.5+(J-1)/LAMBDA)^2-(J-1)/LAMBDA/LAMBDA)*H1#

11100　 L2#=L2#*LAMBDA/2/J/C*O#/I#*Z

11110　 '

11120　 H1#=H2#

11130　 SS#=S0#

11140　 S0#=S0#+L0#

11150　 S1#=S1#+L1#

11160　 S2#=S2#+L2#

11170　 J=J+1

11180 WEND

11190 L0=S0#:L1=S1#:L2=S2#

11200 RETURN

11210 '

11220 ' ** NON-CENTRAL T-DISTRIBUTION **

11230 '

11240 *NT.PROB.HGF.LAMBDA:

11250 PI#=3.141592653589732#

11260 UUN=-LAMBDA:GOSUB *NORMAL.PROB.HGF:T00=P0

11270 UUN=-LAMBDA:T01=-EXP(-UUN*UUN/2)/SQR(PI#*2)

11280 UUN=-LAMBDA:T02=-EXP(-UUN*UUN/2)/SQR(PI#*2)*UUN

11290 UUN= LAMBDA:'GOSUB *NORMAL.PROB.HGF

11300 UUN= LAMBDA:TTT=1-P0-.5

11310 '

11320 DFX=DF/2:DFY=.5:GOSUB *T.PROB.HGF.SUB

11330 GOSUB *POISSON.MIX.T1.LAMBDA:T10=L0:T11=L1:T12=L2

11340 '

11350 DFX=DF/2:DFY= 1:GOSUB *T.PROB.HGF.SUB

11360 GOSUB *POISSON.MIX.T2.LAMBDA:T20=L0:T21=L1:T22=L2

11370 '

11380 IF UUT>0 THEN L0=T00+T10+T20:L1=T01+T11+T21:L2=T02+T12+T22

11390 IF UUT=0 THEN L0=T00　　+T20:L1=T01　　+T21:L2=T02　　+T22

11400 IF UUT<0 THEN L0=T00-T10+T20:L1=T01-T11+T21:L2=T02-T12+T22

11410 RETURN

11420 '

11430 ' *** ポアソン混合 ***

11440 '

11450 *POISSON.MIX.T1.LAMBDA:

11460 EPSL=1E-09

11470 Z=DF/(DF+UUT*UUT):P0#=P0:H1#=G1#

11480 LAMBDA2=LAMBDA*LAMBDA

11490 L0#= EXP(-LAMBDA2/2)*P0#/2

11500 L1#=-EXP(-LAMBDA2/2)*P0#/2*LAMBDA

11510 L2#= EXP(-LAMBDA2/2)*P0#/2*(LAMBDA2-1)

11520 SS#=0:S0#=L0#:S1#=L1#:S2#=L2#

11530 J=1

11540 WHILE ABS((SS#-S0#)/S0#)>EPSL

11550　 A=A0+J-1:B=B0:C=C0

11560　 H2#=-((C-A-1)*H1#+(B-C)*1#)/(A-B+1)/(1-Z)

11570　 '

11580　 O#=H2#:I#=H1#

11590　 L0#=L0#*LAMBDA2/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

11600　 '

11610　 O#=(-LAMBDA+J*2/LAMBDA)*H2#

11620　 I#=(-LAMBDA+(J-1)*2/LAMBDA)*H1#

11630　 L1#=L1#*LAMBDA2/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

11640　 '

11650　 O#=((-LAMBDA+J*2/LAMBDA)^2-1-J*2/LAMBDA)*H2#

11660　 I#=((-LAMBDA+(J-1)*2/LAMBDA)^2-1-(J-1)*2/LAMBDA/LAMBDA)*H1#

11670　 L2#=L2#*LAMBDA2/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

11680　 '

11690　 H1#=H2#

11700　 SS#=S0#

11710　 S0#=S0#+L0#

11720　 S1#=S1#+L1#

11730　 S2#=S2#+L2#

11740　 J=J+1

11750 WEND

11760 L0=S0#:L0=.5-L0

11770 L1=S1#:L1=-L1

11780 L2=S2#:L2=-L2

11790 RETURN

11800 '

11810 ' *** ポアソン混合 ***

11820 '

11830 *POISSON.MIX.T2.LAMBDA:

11840 EPSL=1E-09

11850 Z=DF/(DF+UUT*UUT):P0#=P0:H1#=G1#

11860 LAMBDA2=LAMBDA*LAMBDA

11870 L0#= EXP(-LAMBDA2/2)*P0#*LAMBDA/SQR(PI#*2)

11880 L1#=-EXP(-LAMBDA2/2)*P0#*LAMBDA/SQR(PI#*2)*LAMBDA

11890 L2#= EXP(-LAMBDA2/2)*P0#*LAMBDA/SQR(PI#*2)*(LAMBDA2-1)

11900 SS#=0:S0#=L0#:S1#=L1#:S2#=L2#

11910 J=1

11920 WHILE ABS((SS#-S0#)/S0#)>EPSL

11930　 A=A0+J-1:B=B0:C=C0

11940　 H2#=-((C-A-1)*H1#+(B-C)*1#)/(A-B+1)/(1-Z)

11950　 '

11960　 O#=H2#:I#=H1#

11970　 L0#=L0#*LAMBDA2/2/(J+.5)*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

11980　 '

11990　 O#=(-LAMBDA+J*2/LAMBDA)*H2#

12000　 I#=(-LAMBDA+(J-1)*2/LAMBDA)*H1#

12010　 L1#=L1#*LAMBDA2/2/(J+.5)*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

12020　 '

12030　 O#=((-LAMBDA+J*2/LAMBDA)^2-1-J*2/LAMBDA/LAMBDA)*H2#

12040　 I#=((-LAMBDA+(J-1)*2/LAMBDA)^2-1-(J-1)*2/LAMBDA/LAMBDA)*H1#

12050　 L2#=L2#*LAMBDA2/2/(J+.5)*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

12060　 '

12070　 H1#=H2#

12080　 SS#=S0#

12090　 S0#=S0#+L0#

12100　 S1#=S1#+L1#

12110　 S2#=S2#+L2#

12120　 J=J+1

12130 WEND

12140 L0=S0#:L0=TTT-L0

12150 L1=S1#:L1=-L1

12160 L2=S2#:L2=-L2

12170 RETURN

12180 '

12190 ' ** NON-CENTRAL F-DISTRIBUTION **

12200 '

12210 *NF.PROB.HGF.LAMBDA:

12220 GOSUB *F.PROB.HGF.SUB

12230 GOSUB *POISSON.MIX.F.LAMBDA

12240 RETURN

12250 '

12260 ' *** ポアソン混合 ***

12270 '

12280 *POISSON.MIX.F.LAMBDA:

12290 EPSL=1E-09

12300 Z=DF2/(DF2+DF1*UUF):P0#=P0:H1#=G1#

12310 L0#= EXP(-LAMBDA/2)*P0#

12320 L1#=-EXP(-LAMBDA/2)*P0#/2

12330 L2#= EXP(-LAMBDA/2)*P0#/4

12340 SS#=0:S0#=L0#:S1#=L1#:S2#=L2#

12350 J=1

12360 WHILE ABS((SS#-S0#)/S0#)>EPSL

12370　 A=A0+J-1:B=B0:C=C0

12380　 H2#=-((C-A-1)*H1#+(B-C)*1#)/(A-B+1)/(1-Z)

12390　 '

12400　 O#=H2#:I#=H1#

12410　 L0#=L0#*LAMBDA/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

12420　 '

12430　 O#=(-.5+J/LAMBDA)*H2#

12440　 I#=(-.5+(J-1)/LAMBDA)*H1#

12450　 L1#=L1#*LAMBDA/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

12460　 '

12470　 O#=((-.5+J/LAMBDA)^2-J/LAMBDA/LAMBDA)*H2#

12480　 I#=((-.5+(J-1)/LAMBDA)^2-(J-1)/LAMBDA/LAMBDA)*H1#

12490　 L2#=L2#*LAMBDA/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

12500　 '

12510　 H1#=H2#

12520　 SS#=S0#

12530　 S0#=S0#+L0#

12540　 S1#=S1#+L1#

12550　 S2#=S2#+L2#

12560　 J=J+1

12570 WEND

12580 L0=S0#:L0=1-L0

12590 L1=S1#:L1=-L1

12600 L2=S2#:L2=-L2

12610 RETURN

12620 '

12630 ' ** NON-CENTRAL T-DISTRIBUTION **

12640 '

12650 *NT.LAMBDA:

12660 GOSUB *NORMAL.PERCENT1

12670 D=1-1/DF/4+1/DF/DF/32

12680 AA2=SQR(1+UUT*UUT*(1-D*D))

12690 LAMBDA=AA2*ABS(UUN)+UUT*D

12700 RETURN

12710 '

12720 ' ** NON-CENTRAL F-DISTRIBUTION **

12730 '

12740 *NF.LAMBDA:

12750 GOSUB *NORMAL.PERCENT1

12760 AA1=DF1*(UUF-1)+1+UUN*UUN

12770 AA2=(DF1*(UUF-1)+1)^2-DF1*UUF*(1+DF1/DF2*UUF)*UUN*UUN*2

12780 DI=AA1*AA1-AA2

12790 LAMBDA=AA1+SQR(DI)

12800 RETURN

12810 '

12820 ' **　NON-CENTRAL CHI-SQUARE DISTRIBUTION **

12830 '

12840 *NCHI2.LAMBDA:

12850 GOSUB *NORMAL.PERCENT1

12860 AA1=UUX-DF+1+UUN*UUN

12870 AA2=(UUX-DF+1)^2-UUX*UUN*UUN*2

12880 DI=AA1*AA1-AA2

12890 LAMBDA=AA1+SQR(DI)

12900 RETURN

　

　分布関数Ｆ（ｘ，ν，λ）をλについて微分しているだけなので，大した説明は入らないと思われるが，１つだけ注意点を述べておく．

　

　非心ｔ分布，非心Ｆ分布では，不完全ベータ関数の計算を

　　Ｂx（ｐ，ｑ）／Β（ｐ，ｑ）＝１−Ｂ1-x（ｑ，ｐ）／Β（ｑ，ｐ）

と変数変換して行っているので，１１７７０〜８０行，１２１５０〜６０行，１２５９０〜６０行では負号がつくことに注意されたい．

　なお，Β（ｐ，ｑ）＝Β（ｑ，ｐ）であるが，Ｂx（ｐ，ｑ）≠Ｂx（ｑ，ｐ）である．

　

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

　

【まとめ】

　

　これまで，超幾何関数を用いた確率計算プログラムを開発しながら，このシリーズを書いてきたが，所期に予定していた到達点：

　　ａ）正規分布・χ2分布・ｔ分布・Ｆ分布について

　　　　　　パーセント点　←→　下側確率の双方向の計算

　　ｂ）非心χ2分布・非心ｔ分布・非心Ｆ分布について

　　　　　　パーセント点　←→　下側確率の双方向の計算

　　ｃ）非心χ2分布・非心ｔ分布・非心Ｆ分布について

　　　　　　パーセント点，下側確率　→　非心度の計算

は今回でオールクリアである．

　

　これらのプログラムを組み合わせて用いれば，小数点自由度に対しても，統計数値表にある主要な計算結果を得ることが可能となる．ａ，ｂについては前回（その５）で掲げたが，使用の際には，このプログラムｃをａ，ｂにマージする必要がある．

　

　ソフトを利用する方法ばかりを覚え考えないで使うと，考える力は奪われ，頭脳は確実に退化するだろう．コンピュータは魔物であるが，しかし，考えながら使えばこれほど強力なものもない．このシリーズではどのように計算するのかという原理を具体的（concrete）かつ明示的（explicit）に書いたので，ぜひとも有効に利用して欲しい．

　

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