（その２）に掲載したプログラムは，その後，いろいろな使われ方をしていることがわかった．ところが，先日，このプログラムに対する苦情が寄せられてきた．データを処理したところ，半径P(1)がマイナスの方向に発散していき，エラーで止まってしまうというものである．

　

　最小２乗法ソフト「耕太郎」本体の方には，たとえば，半径のように負にならないパラメータに対して，非負なる制約条件をつけたフィッティング機能が付いているし，初期値がかなりずれていても安定に収束させるための計算アルゴリズムも装備しているのだが，（その２）に載せたプログラムではそれらのサブルーチンを省いている．

　

　当初はそのためのエラーかと思っていたら，どうもそうではないらしいことが判明した．添付データを見たところ，そのデータは楕円の弧のほんの一部分だけにすぎなかったのである．メールを交換してみると，楕円のデータは全周にわたってある場合と中心線より片側に集中する場合があり，現場で特に必要としているのは後者の方であるということがわかった．

　

　しかるに，ＨＰに掲載した非線形最小２乗法のプログラムの初期値サブルーチンは，ほぼ全周にわたってデータ点が採取されているときの設定になっているので，そのサブルーチンを使う限り，苦情にあげたようなエラーは避けられない．

　

　対策としては，

　　(1)パラメータの値を初期値サブルーチンを使わないで，直接入力する

　　(2)初期値に関わらないでパラメータを計算できる線形最小２乗計算を用いて計算してから，その値を初期値として非線形最小２乗計算に移る

であるが，(1)のように自動計算してくれないようではユーザーの手を煩わせることになり，あまり推奨できない．最も忌み嫌われるプログラムなのである．

　

　そこで，今回のコラムではデータが片側あるいは弧の一部に集中する場合であっても，手入力なしに自動計算できるようにプログラムを書き換えてみることにした．

　

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

　

［１］あてはめ例

　

　

　なお，この計算例では，P(1)＜P(2)となりましたが，プログラム上はどちらかを長半径にどちらかを短半径に特定しているわけではなく，大きい方が長半径ということであって，P(1)＜P(2)となっても実害は生じません．

　

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

　

［２］プログラムの解説

　

　中心が原点ではなく，しかも，傾きのある２次元楕円は，５つのパラメータをもつ曲線

　　ｆ（ｘ，ｙ）＝ａ1ｘ^2 ＋ａ2ｘｙ＋ａ3ｙ^2＋ａ4ｘ＋ａ5ｙ−１＝０

の形に書けます．この形であれば，多変量の線形最小２乗法が適用可能です．非線形のプログラムでも線形最小２乗法は解けるのですが，その逆は不可能です．

　

　一方，楕円の基本情報，すなわち，長半径・短半径をａ，ｂ，長軸とｘ軸のなす角度をθ，楕円の中心を（ｘ0，ｙ0）として，楕円モデル式を

　　｛（ｘ−ｘ0）ｃｏｓθ＋（ｙ−ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ−ｘ0）ｓｉｎθ−（ｙ−ｙ0）ｃｏｓθ｝^2／ｂ^2＝１

とすると，未定係数が三角関数の中に含まれているため，どのように変形してもパラメータに関して線形とはなりません．そのため，非線形最小２乗法によるあてはめが必要になります．

　

　しかし，このモデル式では，楕円を媒介変数表示の形に表せるので，描画が短時間で済むようになりますし，現場からの要求である楕円基本情報

　　ａ，ｂ，θ，ｘ0，ｙ0

とその推定誤差

　　Δａ，Δｂ，Δθ，Δｘ0，Δｙ0

を求めることにも応えることができるのです．

　

　したがって，まず最初に

　　ｆ（ｘ，ｙ）＝ａ1ｘ^2 ＋ａ2ｘｙ＋ａ3ｙ^2＋ａ4ｘ＋ａ5ｙ−１

に対して，多変量線形最小２乗法によるあてはめを行い，そのあとで，

　　ａ1，ａ2，ａ3，ａ4，ａ5　→　ａ，ｂ，θ，ｘ0，ｙ0

の変換式を利用して，

　　｛（ｘ−ｘ0）ｃｏｓθ＋（ｙ−ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ−ｘ0）ｓｉｎθ−（ｙ−ｙ0）ｃｏｓθ｝^2／ｂ^2＝１

の形に直してから，パラメータの推定誤差を求め，描画を行うという２ステップの計算方法が考えられます．

　

　そのため，これまでの初期値サブルーチン（2530行〜2750行）を廃止して，主に8520行以下を新規追加しました．パラメータの計算は，線形最小２乗法が済んだ時点でほぼ収束しているのですが，さらに時間をかけて，パラメータの誤差計算に移ります．計算に時間を要するのですが，非線形最小２乗法をあてはめのためにではなく，主にパラメータの誤差値を求めるために用いるという贅沢な方法ともいえるでしょう．

　

　このような書き換えによって，弧の一部しかないデータの場合でも安定に収束するようになり，計算が発散するのを未然に防ぐことができるようになりました．なお，プログラムの簡単な書き換えにより，楕円面のあてはめ問題も解くことができるようになります．

　

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

　

［３］プログラム

　

1000 '

1010 ' **** elliptic approximation ****

1020 '　　　　　　　　　　　　　2002/08/23 (C) ｻﾄｳ ｲｸﾛｳ

1030 SCREEN 3,0:CONSOLE ,,0,1

1040 CLS 3:WIDTH 80,25:COLOR 6

1050 GOSUB *INITIALIZE

1060 GOSUB *READ.DATA

1070 GOSUB *INIT.VAL

1080 GOSUB *EXECUTION

1090 GOSUB *SCATTER.DIAGRAM

1100 GOSUB *REPORT

1110 END

1120 '

1130 ' *** 初期設定 ***

1140 '

1150 *INITIALIZE:

1160 M=5:' No. of parameters

1170 DIM P(M),W(M)

1180 DIM A(M)

1190 DIM P0(M),P1(M),DP(M)

1200 DIM AZ(M,M),BZ(M,M)

1210 DIM HESSE(M,M)

1220 DIM SD(M),SE(M)

1230 DIM COV(M,M)

1240 '

1250 NV=2:' No. of independent variables

1260 DIM X(NV)

1270 DIM AVE(NV),DEV(NV),SEM(NV)

1280 DIM XMAX(NV),XMIN(NV)

1290 FOR J=0 TO NV

1300　 XMAX(J)=-1E+10:XMIN(J)=1E+10

1310 NEXT J

1320 '

1330 RX=.7:RY=.8:BOX=1:KIND=1:POWER=0

1340 IX=.7:IY=.8:MARK=1:CLR=2:DRAW.COLOR=4

1350 TXT.COLOR=6:GRP.COLOR=5:REF.COLOR=6

1360 SCALE.KIND=1:PLOT=1

1370 PI=3.14159

1380 'JX=.709:JY=1

1390 JX=1　 :JY=1

1400 RX=IX*JX:RY=IY*JY

1410 '

1420 FORM$="(((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))/P(1))^2+(((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))/P(2))^2-1"

1430 RETURN

1440 '

1450 ' *** DEFINED EQUATION ***

1460 '

1470 *DEFINE.FORMULA:

1480 DEF FNY(X)=(((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))/P(1))^2+(((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))/P(2))^2-1

1490 RETURN

1500 '

1510 *DEFINE.FORMULA2:

1520 ON J GOSUB *J1,*J2,*J3,*J4,*J5

1530 RETURN

1540 '

1550 *J1:

1560 DEF FNJ(X)=-2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))^2/P(1)^3

1570 RETURN

1580 *J2:

1590 DEF FNJ(X)=-2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))^2/P(2)^3

1600 RETURN

1610 *J3:

1620 DEF FNJ(X)=2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))*(-(X(1)-P(4))*SIN(P(3))+(X(2)-P(5))*COS(P(3)))/P(1)^2+2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))/P(2)^2

1630 RETURN

1640 *J4:

1650 DEF FNJ(X)=-2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))*COS(P(3))/P(1)^2-2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))*SIN(P(3))/P(2)^2

1660 RETURN

1670 *J5:

1680 DEF FNJ(X)=-2*((X(1)-P(4))*COS(P(3))+(X(2)-P(5))*SIN(P(3)))*SIN(P(3))/P(1)^2+2*((X(1)-P(4))*SIN(P(3))-(X(2)-P(5))*COS(P(3)))*COS(P(3))/P(2)^2

1690 RETURN

1980 '

1990 ' *** データ読み込み ***

2000 '

2010 *READ.DATA:

2020 CLS 3

2030 RESTORE *D

2040 READ N1

2050 'DFILE$="DATA01.TXT"

2060 'OPEN DFILE$ FOR INPUT AS #1

2070 'INPUT #1,N1

2080 '

2090 DIM XX1(NV,N1),YH(N1),WT(N1)

2100 DIM Y1(N1),E1(N1)

2110 DIM Z1(N1)

2120 DIM JACOBI(M,N1)

2130 FOR J=1 TO N1:WT(J)=1:NEXT J

2140 '

2150 FOR I=1 TO N1

2160　 FOR J=1 TO NV

2170　　 READ XX1(J,I)

2180　　'INPUT #1,XX1(J,I)

2190　 NEXT J

2200 NEXT I

2210 'CLOSE #1

2220 '

2230 PRINT "ただいまデータの取り込み中"

2240 FOR I=1 TO N1

2250 PRINT "*";

2260 Y1(I)=0

2270 Z1(I)=1

2280 X(0)=Y1(I)

2290 FOR J=1 TO NV

2300　 X(J)=XX1(J,I)

2310 NEXT J

2320 GOSUB *CALCULATE

2330 NEXT I

2340 '

2350 FOR J=0 TO NV

2360　 SS=DEV(J)-AVE(J)*AVE(J)/N1

2370　 DEV(J)=SQR(SS/(N1-1))

2380　 SEM(J)=DEV(J)/SQR(N1)

2390　 AVE(J)=AVE(J)/N1

2400 NEXT J

2410 RETURN

2420 '

2430 ' *** CALCULATE ***

2440 '

2450 *CALCULATE:

2460 FOR J=0 TO NV

2470　 IF XMAX(J)<X(J) THEN XMAX(J)=X(J)

2480　 IF XMIN(J)>X(J) THEN XMIN(J)=X(J)

2490　 AVE(J)=AVE(J)+X(J)

2500　 DEV(J)=DEV(J)+X(J)*X(J)

2510 NEXT J

2520 RETURN

2530 '

2540 ' ** 初期値 **

2550 '

2560 '*INIT.VAL:

2570 X0=0:Y0=0

2580 FOR I=1 TO N1

2590　 X=XX1(1,I):Y=XX1(2,I)

2600　 X0=X0+X

2610　 Y0=Y0+Y

2620 NEXT I

2630 X0=X0/N1:Y0=Y0/N1

2640 P(4)=X0:P(5)=Y0

2650 '

2660 MAX=-1E+10:MIN=1E+10

2670 FOR I=1 TO N1

2680　 X=XX1(1,I):Y=XX1(2,I)

2690　 D=SQR((X-X0)^2+(Y-Y0)^2)

2700　 IF MAX<D THEN MAX=D:XMAX=X-X0:YMAX=Y-Y0

2710　 IF MIN>D THEN MIN=D

2720 NEXT I

2730 P(1)=MAX:P(2)=MIN

2740 IF XMAX=0 THEN P(3)=PI/2 ELSE P(3)=ATN(YMAX/XMAX)

2750 RETURN

2760 '

2770 ' *** 非線形最小２乗法 ***

2780 '

2790 *EXECUTION:

2800 PRINT:PRINT "ただいま計算中"

2810 GOSUB *GAUSS.NEWTON

2820 GOSUB *FINAL

2830 RETURN

2840 '

2850 ' *** GAUSS-NEWTON METHOD ***

2860 '

2870 *GAUSS.NEWTON:

2880 LMAX=50

2890 LOOP=0

2900 GOSUB *RESIDUAL:SO=SS

2910 GOSUB *TEMPORARY

2920 '

2930 GOSUB *JACOBIAN

2940 GOSUB *HESSIAN

2950 FOR I=1 TO M:P1(I)=P(I):NEXT I

2960 GOSUB *MARQUARDT

2970 GOSUB *UPDATE

2980 GOSUB *TRANSFORM

2990 '

3000 IF SN<SO THEN 3020

3010 FOR I=1 TO M:DP(I)=DP(I)/2:NEXT I:GOSUB *TRANSFORM:' [DAMPING]

3020 GOSUB *CONVERGENCE:IF SW=0 THEN RETURN

3030 SO=SN:IF LOOP<LMAX THEN 2930

3040 RETURN

3050 '

3060 ' ** X **

3070 '

3080 *JACOBIAN:

3090 FOR I=1 TO N1

3100　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

3110　 FOR J=1 TO M

3120　　 GOSUB *DEFINE.FORMULA2

3130　　 JACOBI(J,I)=FNJ(X)

3140　 NEXT J

3150 NEXT I

3160 GOSUB *JACOBIAN.SUB

3170 RETURN

3180 '

3190 ' ** Y **

3200 '

3210 *JACOBIAN.SUB:

3220 FOR I=1 TO N1

3230　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

3240　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

3250　 E1(I)=Y1(I)-YH

3260　 YH(I)=YH

3270 NEXT I

3280 '

3290 ' ** XT*Y **

3300 '

3310 FOR J=1 TO M

3320　 S=0

3330　 FOR I=1 TO N1

3340　　 S=S+JACOBI(J,I)*E1(I)*WT(I)

3350　 NEXT I

3360　 W(J)=S

3370 NEXT J

3380 RETURN

3390 '

3400 ' ** XT*X **

3410 '

3420 *HESSIAN:

3430 FOR M1=1 TO M

3440　 FOR M2=M1 TO M

3450　　 S=0

3460　　 FOR I=1 TO N1

3470　　　 S=S+JACOBI(M1,I)*JACOBI(M2,I)*WT(I)

3480　　 NEXT I

3490　　 HESSE(M1,M2)=S:HESSE(M2,M1)=S

3500　 NEXT M2

3510 NEXT M1

3520 NZ=M

3530 RETURN

3540 '

3550 ' ** MARQUARDT **

3560 '

3570 *MARQUARDT:

3580 FOR I=1 TO M

3590　 FOR J=1 TO M:AZ(I,J)=HESSE(I,J):NEXT J

3600 NEXT I

3610 RETURN

3620 '

3630 ' ** UPDATE **

3640 '

3650 *UPDATE:

3660 GOSUB *INVERSE:' (XT*X)^(-1)

3670 GOSUB *UPDATE.SUB

3680 RETURN

3690 '

3700 ' ** (XT*X)^(-1)*XT*Y **

3710 '

3720 *UPDATE.SUB:

3730 FOR J=1 TO M

3740　 S=0

3750　 FOR K=1 TO M

3760　　 S=S+BZ(J,K)*W(K)

3770　 NEXT K

3780　 DP(J)=S

3790 NEXT J

3800 RETURN

3810 '

3820 ' ** TRANSFORM of PARAMETERS **

3830 '

3840 *TRANSFORM:

3850 FOR I=1 TO M

3860　 P(I)=P1(I)+DP(I)

3870 NEXT I

3880 GOSUB *RESIDUAL:SN=SS

3890 RETURN

3900 '

3910 ' ** SUM of RESIDUALS **

3920 '

3930 *RESIDUAL:

3940 SS=0

3950 FOR K=1 TO N1

3960　 FOR V=1 TO NV:X(V)=XX1(V,K):NEXT V

3970　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

3980　 SS=SS+WT(K)*(Y1(K)-YH)^2

3990 NEXT K

4000 RETURN

4010 '

4020 ' *** CHECK of CONVERGENCE ***

4030 '

4040 *CONVERGENCE:

4050 'EPS=.0001

4060 EPS=.001

4070 SW=0

4080 IF ABS((SN-SO)/SO)>EPS THEN SW=1

4090 GOSUB *TEMPORARY

4100 RETURN

4110 '

4120 ' *** TEMPORARY OUTPUT ***

4130 '

4140 *TEMPORARY:

4150 LOOP=LOOP+1

4160 CLS 3

4170 PRINT "LOOP=";LOOP

4180 PRINT

4190 PRINT "PARAMETERS: "

4200 FOR I=1 TO M

4210　 PRINT "P(";I;")=";P(I)

4220 NEXT I

4230 PRINT

4240 PRINT "RESIDUALS=";SS

4250 RETURN

4260 '

4270 ' *** FINAL OUTPUT ***

4280 '

4290 *FINAL:

4300 CLS 3

4310 PRINT "LOOP=";LOOP

4320 PRINT

4330 PRINT "PARAMETERS: "

4340 FOR I=1 TO M

4350　 PRINT "P(";I;")=";P(I)

4360 NEXT I

4370 PRINT

4380 PRINT "RESIDUALS=";SN

4390 '

4400 GOSUB *PROPAGATION

4410 PRINT

4420 PRINT "何かキーを押してください"

4430 WHILE INKEY$="":WEND

4440 RETURN

4450 '

4460 ' ** INVERSE MATRIX ( AZ==>BZ ) **

4470 '

4480 *INVERSE:' [GAUSS-JORDAN]

4490 FOR I=1 TO NZ

4500　 FOR J=1 TO NZ:BZ(I,J)=0:NEXT J

4510　 BZ(I,I)=1

4520 NEXT I

4530 FOR I=1 TO NZ

4540　 IZ=I

4550　 IF AZ(I,I)=0 THEN 4560 ELSE 4620

4560　 IZ=IZ+1:IF IZ>NZ THEN RZ=1:RETURN

4570　 IF AZ(IZ,I)=0 THEN 4560

4580　 FOR J=1 TO NZ

4590　　 SWAP AZ(I,J),AZ(IZ,J)

4600　　 SWAP BZ(I,J),BZ(IZ,J)

4610　 NEXT J

4620　 AII=AZ(I,I)

4630　 FOR J=1 TO NZ

4640　　 AZ(I,J)=AZ(I,J)/AII

4650　　 BZ(I,J)=BZ(I,J)/AII

4660　 NEXT J

4670　 FOR K=1 TO NZ

4680　　 AKI=AZ(K,I):IF K=I THEN 4730

4690　　 FOR J=1 TO NZ

4700　　　 AZ(K,J)=AZ(K,J)-AKI*AZ(I,J)

4710　　　 BZ(K,J)=BZ(K,J)-AKI*BZ(I,J)

4720　　 NEXT J

4730　 NEXT K

4740 NEXT I

4750 RZ=0

4760 RETURN

4770 '

4780 ' *** 係数の誤差変換 ***

4790 '

4800 *PROPAGATION:

4810 FOR IS=1 TO M

4820　 SD(IS)=SQR(ABS(BZ(IS,IS)*SN/(N1-M)))

4830　 SE(IS)=SD(IS)/SQR(N1)

4840　 '

4850　 FOR JS=IS TO M

4860　　 COV(IS,JS)=BZ(IS,JS)*SN/(N1-M):COV(JS,IS)=COV(IS,JS)

4870　 NEXT JS

4880 NEXT IS

4890 RETURN

4900 '

4910 ' *** SCATTER DIAGRAM ***

4920 '

4930 *SCATTER.DIAGRAM:

4940 CLS 3:WIDTH 80,25:COLOR 6

4950 'DMIN=XMIN(1):DMAX=XMAX(1)

4960 'PX=SQR(P(1)*P(1)*COS(P(3))*COS(P(3))+P(2)*P(2)*SIN(P(3))*SIN(P(3)))

4970 'DMIN=P(4)-PX:DMAX=P(4)+PX

4980 IF P(1)>P(2) THEN LARGE=P(1) ELSE LARGE=P(2)

4990 DMIN=P(4)-LARGE:DMAX=P(4)+LARGE

5000 GOSUB *LINER.SCALE:UO=SMIN:UL=SMAX:DU=DIF:SU=STP

5010 '

5020 'DMIN=XMIN(2):DMAX=XMAX(2)

5030 'PY=SQR(P(1)*P(1)*SIN(P(3))*SIN(P(3))+P(2)*P(2)*COS(P(3))*COS(P(3)))

5040 'DMIN=P(5)-PY:DMAX=P(5)+PY

5050 DMIN=P(5)-LARGE:DMAX=P(5)+LARGE

5060 GOSUB *LINER.SCALE:VO=SMIN:VL=SMAX:DV=DIF:SV=STP

5070 GOSUB *DRAW.AXIS

5080 GOSUB *DRAW.ELLIPSE:'GOSUB *GSAVE

5090 GOSUB *TEMP

5100 '

5110 CLS 3

5120 RETURN

5130 '

5140 ' *** 楕円 ***

5150 '

5160 *DRAW.ELLIPSE:

5170 'DRAW.COLOR=4

5180 WUO=UO:WUL=UL

5190 WVO=VO:WVL=VL

5200 WINDOW(WUO,-WVL)-(WUL,-WVO)

5210 VIEW(SX1,SY1)-(SX2,SY2)

5220 '

5230 GOSUB *CONTOUR

5240 '

5250 WINDOW(0,0)-(639,399)

5260 VIEW(0,0)-(639,399)

5270 RETURN

5280 '

5290 ' *** 輪郭 ***

5300 '

5310 *CONTOUR:

5320 C=COS(P(3)):S=SIN(P(3))

5330 AAA=P(1)

5340 BBB=P(2)

5350 X0=P(4):Y0=P(5)

5360 '

5370 G=0

5380 FOR ANGLE=0 TO 2*PI+.01 STEP PI/180

5390　 GX=AAA*COS(ANGLE)*C-BBB*SIN(ANGLE)*S+X0

5400　 GY=AAA*COS(ANGLE)*S+BBB*SIN(ANGLE)*C+Y0

5410　 IF G=0 THEN PSET(GX,-GY),DRAW.COLOR:G=1

5420　 LINE-(GX,-GY),DRAW.COLOR

5430 NEXT ANGLE

5440 RETURN

5450 '

5460 ' *** 座標軸 ***

5470 '

5480 *DRAW.AXIS:

5490 CLS 3

5500 SX1=59:SY1=360*(1-RY):SX2=580*RX+59:SY2=360

5510 LINE(SX1,SY1)-(SX2,SY2),GRP.COLOR,B

5520 RR=3:GOSUB *X.SCALE2:GOSUB *X.REF2

5530 RR=3:GOSUB *Y.SCALE2:GOSUB *Y.REF2

5540 RR=3:GOSUB *DOT.PLOT

5550 GOSUB *TITLE.BACK2

5560 RETURN

5570 '

5580 ' *** 等分目盛り (X) ***

5590 '

5600 *X.SCALE2:

5610 CU=(SX2-SX1)/DU

5620 IF KIND=0 THEN KARA=-RR:MADE=0

5630 IF KIND=1 THEN KARA=0　:MADE=RR

5640 IF KIND=2 THEN KARA=-RR:MADE=RR

5650 FOR K=0 TO DU

5660　 LINE(CU*K+SX1,SY2+KARA)-(CU*K+SX1,SY2+MADE),GRP.COLOR

5670　 LINE(CU*K+SX1,SY1-KARA)-(CU*K+SX1,SY1-MADE),GRP.COLOR

5680 NEXT K

5690 RETURN

5700 '

5710 ' *** 参照値の記入 (X) ***

5720 '

5730 *X.REF2:

5740 REF.COLOR=6

5750 UW=UL-UO

5760 FOR K=0 TO DU STEP SU

5770　 F=UW*K/DU+UO:F$=STR$(F)

5780　 IF F>=0 THEN F$=MID$(F$,2)

5790　 FL=LEN(F$)

5800　 FOR L=1 TO FL

5810　　 GX=CU*K+59-FL*8/2+(L-1)*8:GY=23*16

5820　　 PUT(GX,GY),KANJI(ASC(MID$(F$,L,1))),PSET,REF.COLOR,0

5830　 NEXT L

5840 NEXT K

5850 RETURN

5860 '

5870 ' *** 等分目盛り (Y) ***

5880 '

5890 *Y.SCALE2:

5900 CV=(SY2-SY1)/DV

5910 IF KIND=0 THEN KARA=0　:MADE=RR

5920 IF KIND=1 THEN KARA=-RR:MADE=0

5930 IF KIND=2 THEN KARA=-RR:MADE=RR

5940 FOR K=0 TO DV

5950　 LINE(SX1+KARA,SY2-CV*K)-(SX1+MADE,SY2-CV*K),GRP.COLOR

5960　 LINE(SX2-KARA,SY2-CV*K)-(SX2-MADE,SY2-CV*K),GRP.COLOR

5970 NEXT K

5980 RETURN

5990 '

6000 ' *** 参照値の記入 (Y) ***

6010 '

6020 *Y.REF2:

6030 REF.COLOR=6

6040 VW=VL-VO

6050 FOR K=0 TO DV STEP SV

6060　 F=VW*K/DV+VO:F$=STR$(F)

6070　 IF F>=0 THEN F$=MID$(F$,2)

6080　 FL=LEN(F$)

6090　 FOR L=1 TO FL

6100　　 GX=48-FL*8+(L-1)*8:GY=360-CV*K-8

6110　　 PUT(GX,GY),KANJI(ASC(MID$(F$,L,1))),PSET,REF.COLOR,0

6120　 NEXT L

6130 NEXT K

6140 RETURN

6150 '

6160 ' *** プロット ***

6170 '

6180 *DOT.PLOT:

6190 WUO=UO:WUL=UL

6200 WVO=VO:WVL=VL

6210 WX=WUL-WUO

6220 BX=(WUL*SX1-WUO*SX2)/WX

6230 AX=(SX2-SX1)/WX

6240 WY=WVL-WVO

6250 BY=(WVL*SY2-WVO*SY1)/WY

6260 AY=(SY1-SY2)/WY

6270 '

6280 FOR I=1 TO N1

6290 X=XX1(1,I):Y=XX1(2,I)

6300 XX=AX*X+BX

6310 YY=AY*Y+BY

6320 CIRCLE(XX,YY),RR,1

6330 PAINT(XX,YY),CLR,1

6340 CIRCLE(XX,YY),RR,CLR

6350 NEXT I

6360 RETURN

6370 '

6380 ' *** タイトル ***

6390 '

6400 *TITLE.BACK2:

6410 LX=INT((580*RX+59)/8)

6420 LY=INT(360*(1-RY)/16)

6430 LZ=INT((290*RX+59)/8)

6440 L1$="elliptic approximation by "

6450 L2$="GAUSS-NEWTON METHOD"

6460 L$=L1$+L2$

6470 LOCATE LZ-LEN(L$)\2,LY-3:PRINT L$

6480 '

6490 LOCATE LX+3,LY:PRINT "PARAMETERS: "

6500 FOR K=1 TO M

6510　 V$="P"+MID$(STR$(K),2)+"="

6520　 LOCATE LX+3,LY+K+1:PRINT V$;:PRINT P(K)

6530 NEXT K

6540 RETURN

6550 '

6560 *TEMP:

6570 LOCATE 0,24:PRINT "何かキーを押してください";:WHILE INKEY$="":WEND

6580 RETURN

6590 '

6600 ' *** 平均律スケールの作成 ***

6610 '

6620 *LINER.SCALE:

6630 DIM INTERVAL(10)

6640 RESTORE *REF.INTERVAL

6650 FOR I=1 TO 10:READ INTERVAL(I):NEXT I

6660 MAX.INTERVAL=10

6670 '

6680 ARGUMENT=DMAX-DMIN

6690 GOSUB *SCALE.RTN

6700 INDEX=CMP

6710 INTERVAL=INTERVAL(INDEX)*EP

6720 LO=INT(DMIN/INTERVAL)

6730 HI=-INT(-DMAX/INTERVAL)

6740 WHILE (HI-LO)>MAX.INTERVAL

6750　 INDEX=INDEX+1

6760　 IF INDEX>10 THEN INDEX=1:EP=EP*10

6770　 INTERVAL=INTERVAL(INDEX)*EP

6780　 LO=INT(DMIN/INTERVAL)

6790　 HI=-INT(-DMAX/INTERVAL)

6800 WEND

6810 SMIN=LO*INTERVAL

6820 SMAX=HI*INTERVAL

6830 DIF=HI-LO

6840 STP=1

6850 ERASE INTERVAL

6860 RETURN

6870 '

6880 ' ** ラベリング関数 **

6890 '

6900 *SCALE.RTN:

6910 'RE=INT(LOG(ARGUMENT)/LOG(10))

6920 RE=INT(LOG(CDBL(ARGUMENT))/LOG(10#))

6930 'EP=10^RE

6940 EP=VAL("1E"+STR$(RE)):'[EXPONENT PART]

6950 MP=ARGUMENT/EP　　　 :'[MANTISSA PART]

6960 FMP=INT(MP)　:'[FLOORING MANTISSA PART]

6970 CMP=-INT(-MP):'[CEILING MANTISSA PART]

6980 RETURN

6990 '

7000 ' ** 参照間隔 **

7010 '

7020 *REF.INTERVAL:

7030 DATA .1,.2,.5,.5,.5,1,1,1,1,1

7040 '

7050 ' *** REPORT ***

7060 '

7070 *REPORT:

7080 CLS 3

7090 PRINT:PRINT

7100 PRINT "DATA No.=";N1

7110 PRINT

7120 PRINT " "," MEAN"," S.D."," MIN"," MAX"

7130 FOR I=1 TO NV

7140 PRINT "X("+MID$(STR$(I),2)+")",AVE(I),DEV(I),XMIN(I),XMAX(I)

7150 NEXT I

7160 PRINT:PRINT

7170 '

7180 PRINT "[ DEFINED EQUATION ]"

7190 PRINT "Y= ";FORM$

7200 '

7210 GOSUB *CEF.CONFIDENCE

7220 GOSUB *TEMP.STOP

7230 '

7240 CLS 3

7250 PRINT:PRINT "残差"

7260 PRINT " No."," OBSERVED Y"," EXPECTED Y"," RESIDUAL"

7270 FOR I=1 TO N1

7280　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

7290　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

7300　 PRINT I,Y1(I),YH,Y1(I)-YH

7310　 IF (I MOD 20)=0 THEN PRINT

7320 NEXT I

7330 PRINT

7340 GOSUB *TEMP.STOP

7350 RETURN

7360 '

7370 *CEF.CONFIDENCE:

7380 PRINT

7390 PRINT "回帰係数とその９５％信頼区間"

7400 DF1=1:DF2=N1-M

7410 GOSUB *CEF.CONFIDENCE.SUB

7420 PRINT:PRINT

7430 RETURN

7440 '

7450 *CEF.CONFIDENCE.SUB:

7460 FOR J=1 TO M

7470　 P$="P("+MID$(STR$(J),2)+")= "

7480　 PRINT P$;:PRINT P(J);

7490　 IF SD(J)=0 THEN PRINT:GOTO 7520

7500　 GOSUB *CEF.ERROR

7510　 PRINT TAB(30);"S.E.=";SE(J)　　　　　　　　　　　　　　　　　　　　　　　　　　　　　;TAB(50);P(J)-SE(J)*UUT;" - ";P(J)+SE(J)*UUT

7520 NEXT J

7530 RETURN

7540 '

7550 *TEMP.STOP:

7560 PRINT "何かキーを押してください"

7570 WHILE INKEY$="":WEND

7580 RETURN

7590 '

7600 ' *** ERROR of COEFFICIENT ***

7610 '

7620 *CEF.ERROR:

7630 PP=.05 :GOSUB *F.PERCENT

7640 UUT=SQR(UUF)

7650 RETURN

7660 '

7670 ' ** NORMAL DISTRIBUTION **

7680 '

7690 *NORMAL.PERCENT:

7700 XXX=-LOG(4*PP*(1-PP))

7710 UUN=SQR(XXX*(2.06118-5.72622/(XXX+11.6406)))

7720 IF PP>.5 THEN UUN=-UUN

7730 RETURN

7740 '

7750 ' ** T-DISTRIBUTION **

7760 '

7770 *T.PERCENT:

7780 GOSUB *NORMAL.PERCENT

7790 UUN=ABS(UUN)

7800 UU2=UUN*UUN

7810 AA1=(UU2+1)/DF/4

7820 AA2=((5*UU2+16)*UU2+3)/96/DF/DF

7830 AA3=(((3*UU2+19)*UU2+17)*UU2-15)/384/DF/DF/DF

7840 UUT=UUN*(1+AA1+AA2+AA3)

7850 IF PP>.5 THEN UUT=-UUT

7860 RETURN

7870 '

7880 ' ** F-DISTRIBUTION **

7890 '

7900 *F.PERCENT:

7910 IF DF1=1 THEN DF=DF2:PP=PP/2:GOSUB *T.PERCENT　　　　　　　　　　　　　　　　　　　　　　　　:UUF=UUT^2:PP=PP*2:RETURN

7920 IF DF2=1 THEN DF=DF2:PP=(1-PP)/2:GOSUB *T.PERCENT　　　　　　　　　　　　　　　　　　　　　　:UUF=1/UUT^2:PP=1-PP*2:RETURN

7930 IF DF1=2 THEN UUF=DF2/2*(PP^(-2/DF2)-1):RETURN

7940 IF DF2=2 THEN UUF=2/DF1/((1-PP)^(-2/DF1)-1):RETURN

7950 AA=2/9/DF1:BB=2/9/DF2:CC=1-AA:DD=1-BB

7960 GOSUB *NORMAL.PERCENT

7970 UUF=CC*DD+UUN*SQR(CC*CC*BB+DD*DD*AA-AA*BB*UUN*UUN)

7980 UUF=(UUF/(DD*DD-BB*UUN*UUN))^3

7990 RETURN

8000 '

8010 ' ** CHI-SQUARE DISTRIBUTION **

8020 '

8030 *CHI2.PERCENT:

8040 IF DF=1 THEN PP=PP/2:GOSUB *NORMAL.PERCENT　　　　　　　　　　　　　　　　　　　　　　　　　:UUX=UUN*UUN:PP=PP*2:RETURN

8050 IF DF=2 THEN UUX=-2*LOG(PP):RETURN

8060 GOSUB *NORMAL.PERCENT

8070 UUX=2/(9*DF)

8080 UUX=1-UUX+UUN*SQR(UUX)

8090 UUX=UUX^3*DF

8100 RETURN

8270 '

8280 ' ** 係数変換 **

8290 '

8300 *P2A:

8310 C=COS(P(3)):S=SIN(P(3))

8320 K=1-((P(4)*C+P(5)*S)/P(1))^2-((P(4)*S-P(5)*C)/P(2))^2

8330 A1=(C/P(1))^2+(S/P(2))^2

8340 A2=(1/P(1)^2-1/P(2)^2)*C*S*2

8350 A3=(S/P(1))^2+(C/P(2))^2

8360 A4=-(P(4)*C+P(5)*S)*C*2/P(1)^2-(P(4)*S-P(5)*C)*S*2/P(2)^2

8370 A5=-(P(4)*C+P(5)*S)*S*2/P(1)^2+(P(4)*S-P(5)*C)*C*2/P(2)^2

8380 A(1)=A1/K:A(2)=A2/K:A(3)=A3/K:A(4)=A(4)/K:A(5)=A5/K

8390 RETURN

8400 '

8410 *A2P:

8420 A1=A(1):A2=A(2):A3=A(3):A4=A(4):A5=A(5)

8430 DIF=A1-A3

8440 IF DIF=0 THEN P(3)=PI/4 ELSE P(3)=1/2*ATN(A2/(A1-A3))

8450 C=COS(P(3)):S=SIN(P(3))

8460 K=1/(1+(C*C-S*S)/4*((A4*C+A5*S)^2/(A1*C*C-A3*S*S)-(A4*S-A5*C)^2/(A1*S*S-A3*C*C)))

8470 P(1)=SQR((C*C-S*S)/(A1*C*C-A3*S*S)/K)

8480 P(2)=SQR(-(C*C-S*S)/(A1*S*S-A3*C*C)/K)

8490 P(4)=-(P(1)*P(1)*C*(A4*C+A5*S)+P(2)*P(2)*S*(A4*S-A5*C))/2*K

8500 P(5)=-(P(1)*P(1)*S*(A4*C+A5*S)-P(2)*P(2)*C*(A4*S-A5*C))/2*K

8510 RETURN

8520 '

8530 ' *** DEFINED EQUATION ***

8540 '

8550 *DEFINE.FORMULA3:

8560 ON L GOSUB *L1,*L2,*L3,*L4,*L5

8570 RETURN

8580 '

8590 *L1:DEF FNJ(X)=X(1)*X(1):RETURN

8600 *L2:DEF FNJ(X)=X(1)*X(2):RETURN

8610 *L3:DEF FNJ(X)=X(2)*X(2):RETURN

8620 *L4:DEF FNJ(X)=X(1)　　 :RETURN

8630 *L5:DEF FNJ(X)=X(2)　　 :RETURN

8640 '

8650 ' ** 初期値 **

8660 '

8670 *INIT.VAL:

8680 GOSUB *LEAST.SQUARES

8690 GOSUB *A2P

8700 RETURN

8710 '

8720 ' ** XT*X **

8730 '

8740 *LEAST.SQUARES:

8750 FOR M1=1 TO M

8760　 FOR M2=M1 TO M

8770　 S=0

8780　　 FOR I=1 TO N1

8790　　　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

8800　　　 L=M1:GOSUB *DEFINE.FORMULA3:L1=FNJ(X)

8810　　　 L=M2:GOSUB *DEFINE.FORMULA3:L2=FNJ(X)

8820　　　 S=S+WT(I)*L1*L2

8830　　 NEXT I

8840　 AZ(M1,M2)=S:AZ(M2,M1)=S

8850　 NEXT M2

8860 NEXT M1

8870 NZ=M

8880 '

8890 ' ** (XT*X)^(-1) **

8900 '

8910 GOSUB *INVERSE

8920 '

8930 ' ** XT*Y **

8940 '

8950 FOR J=1 TO M

8960　 S=0

8970　 FOR I=1 TO N1

8980　　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

8990　　 L=J:GOSUB *DEFINE.FORMULA3:L1=FNJ(X)

9000　　 S=S+WT(I)*L1*Z1(I)

9010　 NEXT I

9020　 W(J)=S

9030 NEXT J

9040 '

9050 ' ** (XT*X)^(-1)*XT*Y **

9060 '

9070 FOR J=1 TO M

9080　 S=0

9090　 FOR K=1 TO M

9100　　 S=S+BZ(J,K)*W(K)

9110　 NEXT K

9120　 A(J)=S

9130 NEXT J

9140 RETURN

9150 '

9160 *D:

9170 DATA 125

9180 DATA -165,-7.27991432

9190 DATA -164,-8.939199724

9200 DATA -163,-10.59978079

9210 DATA -162,-12.24757487

9220 DATA -161,-13.88258195

9230 DATA -160,-15.49711288

9240 DATA -159,-17.10654597

9250 DATA -158,-18.70319207

9260 DATA -157,-20.29344468

9270 DATA -156,-21.85682763

9280 DATA -155,-23.42789975

9290 DATA -154,-24.95421738

9300 DATA -153,-26.49332201

9310 DATA -152,-28.01195049

9320 DATA -151,-29.52548114

9330 DATA -150,-31.03901178

9340 DATA -149,-32.52057494

9350 DATA -148,-33.98805544

9360 DATA -147,-35.45043812

9370 DATA -146,-36.90003379

9380 DATA -145,-38.33684248

9390 DATA -144,-39.74678151

9400 DATA -143,-41.15801621

9410 DATA -142,-42.55007042

9420 DATA -141,-43.93573113

9430 DATA -140,-45.29452219

9440 DATA -139,-46.64821541

9450 DATA -138,-47.98912165

9460 DATA -137,-49.31724089

9470 DATA -136,-50.61849048

9480 DATA -135,-51.91464224

9490 DATA -134,-53.198007

9500 DATA -133,-54.46858477

9510 DATA -132,-55.71229289

9520 DATA -131,-56.95729667

9530 DATA -130,-58.17672647

9540 DATA -129,-59.38976277

9550 DATA -128,-60.57592942

9560 DATA -127,-61.75060474

9570 DATA -126,-62.91249307

9580 DATA -125,-64.0615944

9590 DATA -124,-65.18382608

9600 DATA -123,-66.30095993

9610 DATA -122,-67.39891329

9620 DATA -121,-68.47768616

9630 DATA -120,-69.53598288

9640 DATA -119,-70.58278826

9650 DATA -118,-71.61680666

9660 DATA -117,-72.63164456

9670 DATA -116,-73.61961281

9680 DATA -115,-74.59608974

9690 DATA -114,-75.55338617

9700 DATA -113,-76.49789561

9710 DATA -112,-77.4091419

9720 DATA -111,-78.31529036

9730 DATA -110,-79.19586483

9740 DATA -109,-80.05596315

9750 DATA -108,-80.89817664

9760 DATA -107,-81.72120965

9770 DATA -106,-82.52506216

9780 DATA -105,-83.30204502

9790 DATA -104,-84.06114305

9800 DATA -103,-84.8010606

9810 DATA -102,-85.51540416

9820 DATA -101,-86.20287807

9830 DATA -100,-86.87246715

9840 DATA -99,-87.52287574

9850 DATA -98,-88.14771035

9860 DATA -97,-88.7392818

9870 DATA -96,-89.31296844

9880 DATA -95,-89.86108108

9890 DATA -94,-90.39001324

9900 DATA -93,-90.87928875

9910 DATA -92,-91.34428593

9920 DATA -91,-91.79010263

9930 DATA -90,-92.19755835

9940 DATA -89,-92.57814442

9950 DATA -88,-92.92805866

9960 DATA -87,-93.25239893

9970 DATA -86,-93.5447717

9980 DATA -85,-93.79748784

9990 DATA -84,-94.01953215

10000 DATA -83,-94.20960898

10010 DATA -82,-94.36771833

10020 DATA -81,-94.47977753

10030 DATA -80,-94.56116492

10040 DATA -79,-94.60419133

10050 DATA -78,-94.60246326

10060 DATA -77,-94.55468505

10070 DATA -76,-94.46984152

10080 DATA -75,-94.34663701

10090 DATA -74,-94.16589104

10100 DATA -73,-93.93909492

10110 DATA -72,-93.66244649

10120 DATA -71,-93.33465009

10130 DATA -70,-92.95570571

10140 DATA -69,-92.50513721

10150 DATA -68,-92.0111099

10160 DATA -67,-91.44675412

10170 DATA -66,-90.81206988

10180 DATA -65,-90.10576151

10190 DATA -64,-89.33042034

10200 DATA -63,-88.47835721

10210 DATA -62,-87.54317862

10220 DATA -61,-86.50440843

10230 DATA -60,-85.38381844

10240 DATA -59,-84.154539

10250 DATA -58,-82.82296362

10260 DATA -57,-81.35582915

10270 DATA -56,-79.7685139

10280 DATA -55,-78.04054172

10290 DATA -54,-76.15912563

10300 DATA -53,-74.08460897

10310 DATA -52,-71.83237006

10320 DATA -51,-69.34996528

10330 DATA -50,-66.62460762

10340 DATA -49,-63.59106646

10350 DATA -48,-60.21996563

10360 DATA -47,-56.42689403

10370 DATA -46,-52.11594919

10380 DATA -45,-47.11960455

10390 DATA -44,-41.19749855

10400 DATA -43,-33.88800377

10410 DATA -42,-24.11061926

10420 DATA -41,-7.279912252

　

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

　

　データ読み込みサブルーチン（1980行〜2210行）では，プログラム中にDATA文で書き込まれたデータをREAD文により読み込んで計算を行っています．いわば，データをプログラム中に内蔵化する方法です．

　

　しかし，この計算例のように大量のデータを扱わなければならない場合には効率が悪く，外部メモリにデータを格納しておき（ファイル名を例えばdata01.txtとする），必要なときにそこからデータを呼び出して，データ処理とデータ入力を切り離した方が都合がよくなります．以下，データ読み込みサブルーチンの変更方法を示します．

　

1980 '

1990 ' *** データ読み込み ***

2000 '

2010 *READ.DATA:

2020 CLS 3

2030 'RESTORE *D

2040 'READ N1

2050 DFILE$="DATA01.TXT"

2060 OPEN DFILE$ FOR INPUT AS #1

2070 INPUT #1,N1

2080 '

2090 DIM XX1(NV,N1),YH(N1),WT(N1)

2100 DIM Y1(N1),E1(N1)

2110 DIM Z1(N1)

2120 DIM JACOBI(M,N1)

2130 FOR J=1 TO N1:WT(J)=1:NEXT J

2140 '

2150 FOR I=1 TO N1

2160　 FOR J=1 TO NV

2170　　'READ XX1(J,I)

2180　　 INPUT #1,XX1(J,I)

2190　 NEXT J

2200 NEXT I

2210 CLOSE #1

　

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