（その５）では，ｘ，ｙ双方に誤差があるときの取り扱い法として，粟屋の２次元フィッティングをプログラムなしに解説した．また，その原法であるデミング法についてはプログラムと計算例を示したが，決して他人に勧められるような代物ではなかった．

　

　一般的な曲線：ｆ（ｘ，ｙ）＝０に対する粟屋法の近似プログラムは「耕太郎」にサポートされているのだが，そこでは微分でなく差分が使われていること，また，陰関数の描画手法を用いているため，計算と出力に要する時間のかかるものとなっている．そこで，（その６）では必要な部分のみを切り取って粟屋法による楕円近似専用プログラムと円近似専用プログラムを作成し，紹介することにした．

　

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

　

［１］楕円近似プログラム（粟屋法）

　

1000 '

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

1020 '　　　　　　　　　　　　　2003/01/15 (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:L2$="AWAYA METHOD"

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-P(4))*COS(P(3))+(Y-P(5))*SIN(P(3)))/P(1))^2+(((X-P(4))*SIN(P(3))-(Y-P(5))*COS(P(3)))/P(2))^2-1"

1430 RETURN

1440 '

1450 ' *** DEFINED EQUATION ***

1460 '

1470 *DEFINE.FORMULA:

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

1490 DEF FNX(X,Y)=((X-P(4))*COS(P(3))+(Y-P(5))*SIN(P(3)))*COS(P(3))*2/P(1)^2+((X-P(4))*SIN(P(3))+(Y-P(5))*COS(P(3)))*SIN(P(3))*2/P(2)^2

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

1510 RETURN

1520 '

1530 *DEFINE.FORMULA2:

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

1550 RETURN

1560 '

1570 *J1:

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

1590 RETURN

1600 *J2:

1610 DEF FNJ(X,Y)=-2*((X-P(4))*SIN(P(3))-(Y-P(5))*COS(P(3)))^2/P(2)^3

1620 RETURN

1630 *J3:

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

1650 RETURN

1660 *J4:

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

1680 RETURN

1690 *J5:

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

1710 RETURN

1720 '

1730 *D:

1740 DATA 25

1750 DATA 2,0

1760 DATA 7,-2

1770 DATA 14,-3

1780 DATA 18,-4

1790 DATA 23,-4

1800 DATA 31,-3

1810 DATA 38,-1

1820 DATA 41,4

1830 DATA 40,8

1840 DATA 38,12

1850 DATA 34,13

1860 DATA 31,16

1870 DATA 26,17

1880 DATA 19,19

1890 DATA 13,20

1900 DATA 7,20

1910 DATA 2,19

1920 DATA -3,18

1930 DATA -6,17

1940 DATA -9,14

1950 DATA -10,12

1960 DATA -9,8

1970 DATA -7,5

1980 DATA -3,3

1990 DATA 0,1

2000 '

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

2020 '

2030 *READ.DATA:

2040 CLS 3

2050 RESTORE *D

2060 READ N1

2070 'DFILE$="DATA01.TXT"

2080 'OPEN DFILE$ FOR INPUT AS #1

2090 'INPUT #1,N1

2100 '

2110 DIM XX1(NV,N1),YH(N1),WT(N1),XWT(N1),YWT(N1)

2120 DIM Y1(N1),E1(N1)

2130 DIM Z1(N1)

2140 DIM JACOBI(M,N1)

2150 FOR J=1 TO N1:WT(J)=1:XWT(J)=1:YWT(J)=1:NEXT J

2160 DIM X1(N1),XX(N1),YY(N1)

2170 DIM X11(N1),Y11(N1),DX(N1),DY(N1)

2180 DIM F10(M,N1),FX0(N1),FY0(N1)

2190 '

2200 FOR I=1 TO N1

2210　 FOR J=1 TO NV

2220　　 READ XX1(J,I)

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

2240　 NEXT J

2250 NEXT I

2260 'CLOSE #1

2270 '

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

2290 FOR I=1 TO N1

2300 PRINT "*";

2310 Y1(I)=0

2320 X(0)=Y1(I)

2330 FOR J=1 TO NV

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

2350 NEXT J

2360　Z1(I)=1

2370 'Z1(I)=X(1)*X(1)+X(2)*X(2)

2380 GOSUB *CALCULATE

2390 NEXT I

2400 '

2410 FOR J=0 TO NV

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

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

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

2450　 AVE(J)=AVE(J)/N1

2460 NEXT J

2470 RETURN

2480 '

2490 ' *** CALCULATE ***

2500 '

2510 *CALCULATE:

2520 FOR J=0 TO NV

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

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

2550　 AVE(J)=AVE(J)+X(J)

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

2570 NEXT J

2580 RETURN

2590 '

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

2610 '

2620 *EXECUTION:

2630 FOR I=1 TO N1

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

2650　 X1(I)=X(1):Y1(I)=X(2)

2660 NEXT I

2670 FOR I=1 TO N1

2680　 XX(I)=X1(I)

2690　 YY(I)=Y1(I)

2700 NEXT I

2710 '

2720 GOSUB *GAUSS.NEWTON

2730 GOSUB *FINAL

2740 RETURN

2750 '

2760 ' *** 最適化計算 ***

2770 '

2780 *GAUSS.NEWTON:

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

2800 LMAX=200

2810 LOOP=0

2820 GOSUB *RESIDUAL:SO=SS

2830 '

2840 GOSUB *START.2DFIT

2850 GOSUB *HESSIAN

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

2870 FOR I=1 TO N1:X11(I)=X1(I):Y11(I)=Y1(I):NEXT I

2880 GOSUB *MARQUARDT

2890 GOSUB *UPDATE

2900 GOSUB *TRANSFORM

2910 '

2920 IF SN<SO THEN 2970

2930 FOR I=1 TO M:DP(I)=DP(I)/2:NEXT I

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

2950 FOR I=1 TO N1:X1(I)=X11(I):Y1(I)=Y11(I):NEXT I

2960 GOSUB *REFORM:GOSUB *TRANSFORM:' [DAMPING]

2970 GOSUB *CONVERGENCE:IF SW=0 THEN RETURN

2980 SO=SN:IF LOOP<LMAX THEN 2840

2990 RETURN

3000 '

3010 *START.2DFIT:

3020 FOR J=1 TO M

3030　 FOR K=1 TO M:HESSE(J,K)=0:NEXT K

3040　 W(J)=0:DP(J)=0

3050 NEXT J

3060 RETURN

3070 '

3080 *HESSIAN:

3090 FOR I=1 TO N1

3100　 X=X1(I):Y=Y1(I)

3110　 XX=XX(I):YY=YY(I)

3120　 GOSUB *APPROX

3130　 W0=1/(FX*FX/XWT(I)+FY*FY/YWT(I))

3140　 FOR J=1 TO M

3150　　 FOR K=1 TO M:HESSE(J,K)=HESSE(J,K)+W0*F1(J)*F1(K):NEXT K

3160　　 W(J)=W(J)+W0*(F0+(XX-X)*FX+(YY-Y)*FY)*F1(J):' [AWAYA]

3170　 NEXT J

3180 NEXT I

3190 RETURN

3200 '

3210 *UPDATE:

3220 NZ=M:GOSUB *INVERSE

3230 FOR K=1 TO M

3240　 FOR J=1 TO M:DP(K)=DP(K)+BZ(K,J)*W(J):NEXT J

3250 NEXT K

3260 *REFORM:

3270 FOR I=1 TO N1

3280　 X=X1(I):Y=Y1(I)

3290　 XX=XX(I):YY=YY(I)

3300　 GOSUB *APPROX

3310　 W0=1/(FX*FX/XWT(I)+FY*FY/YWT(I))

3320　 SUM=0:FOR K=1 TO M:SUM=SUM+F1(K)*DP(K):NEXT K

3330　 L0=-W0*(F0+(XX-X)*FX+(YY-Y)*FY-SUM):' [AWAYA]

3340　 DX(I)=L0*FX/XWT(I)+XX-X:' [AWAYA]

3350　 DY(I)=L0*FY/YWT(I)+YY-Y:' [AWAYA]

3360 NEXT I

3370 RETURN

3380 '

3390 *TRANSFORM:

3400 FOR I=1 TO N1

3410　 X1(I)=X11(I)+DX(I)

3420　 Y1(I)=Y11(I)+DY(I)

3430 NEXT I

3440 FOR I=1 TO M

3450　 P(I)=P1(I)-DP(I)

3460 NEXT I

3470 GOSUB *RESIDUAL:SN=SS

3480 RETURN

3490 '

3500 ' ** MARQUARDT **

3510 '

3520 *MARQUARDT:

3530 FOR I=1 TO M

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

3550 NEXT I

3560 RETURN

3570 '

3580 ' ** SUM of RESIDUALS **

3590 '

3600 *RESIDUAL:

3610　 SS=0

3620　 FOR I=1 TO N1

3630　　 X=X1(I):Y=Y1(I)

3640　　 XX=XX(I):YY=YY(I)

3650　　 GOSUB *APPROX

3660　　 W0=1/(FX*FX/XWT(I)+FY*FY/YWT(I))

3670　　 SS=SS+W0*(F0+(XX-X)*FX+(YY-Y)*FY)^2　:' [AWAYA]

3680　　'SS=SS+XWT(I)*(XX-X)^2+YWT(I)*(YY-Y)^2:' [DEMING]

3690　 NEXT I

3700 RETURN

3710 '

3720 ' *** CHECK of CONVERGENCE ***

3730 '

3740 *CONVERGENCE:

3750 EPS=.0001

3760 SW=0

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

3780 GOSUB *TEMPORARY

3790 RETURN

3800 '

3810 ' *** TEMPORARY OUTPUT ***

3820 '

3830 *TEMPORARY:

3840 LOOP=LOOP+1

3850 CLS 3

3860 PRINT "LOOP=";LOOP

3870 PRINT

3880 PRINT "PARAMETERS: "

3890 FOR I=1 TO M

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

3910 NEXT I

3920 PRINT

3930 PRINT "RESIDUAL=";SS

3940 RETURN

3950 '

3960 ' *** FINAL OUTPUT ***

3970 '

3980 *FINAL:

3990 CLS 3

4000 PRINT "LOOP=";LOOP

4010 PRINT

4020 PRINT "PARAMETERS: "

4030 FOR I=1 TO M

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

4050 NEXT I

4060 PRINT

4070 PRINT "RESIDUAL=";SN

4080 '

4090 SEE=0

4100 FOR I=1 TO N1

4110　 X=X1(I):Y=Y1(I)

4120　 XX=XX(I):YY=YY(I)

4130　 GOSUB *APPROX

4140　 SEE=SEE+(F0+(XX-X)*FX+(YY-Y)*FY)^2:' [AWAYA]

4150　'SEE=SEE+(XX-X)^2+(YY-Y)^2　　　　 :' [DEMING]

4160 NEXT I

4170 '

4180 'GOSUB *GAUSS.NEWTON.ERROR

4190 GOSUB *PROPAGATION

4200 PRINT

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

4220 WHILE INKEY$="":WEND

4230 RETURN

4240 '

4250 ' *** PROPAGATION of ERROR ***

4260 '

4270 *PROPAGATION:

4280 SS=0

4290 FOR I=1 TO N1

4300　 X=X1(I):Y=Y1(I)

4310　 XX=XX(I):YY=YY(I)

4320　 GOSUB *APPROX

4330　　 FX0(I)=FX:FY0(I)=FY

4340　　 FOR J=1 TO M:F10(J,I)=F1(J):NEXT J

4350　 W0=1/(FX*FX/XWT(I)+FY*FY/YWT(I))

4360　 SS=SS+W0*(F0+(XX-X)*FX+(YY-Y)*FY)^2　:' [AWAYA]

4370　'SS=SS+XWT(I)*(XX-X)^2+YWT(I)*(YY-Y)^2:' [DEMING]

4380 NEXT I

4390 '

4400 FOR J1=1 TO M

4410　 FOR J2=J1 TO M

4420　　 DD=0

4430　　 FOR I=1 TO N1

4440　　　 FX=FX0(I):FY=FY0(I)

4450　　　 W0=1/(FX*FX/XWT(I)+FY*FY/YWT(I))

4460　　　 SUM1=0:SUM2=0

4470　　　 FOR KK=1 TO M

4480　　　　 F1=F10(KK,I)

4490　　　　 SUM1=SUM1+BZ(J1,KK)*F1

4500　　　　 SUM2=SUM2+BZ(J2,KK)*F1

4510　　　 NEXT KK

4520　　　 DD=DD+(W0*FX)^2*SUM1*SUM2/XWT(I)+(W0*FY)^2*SUM1*SUM2/YWT(I)

4530　　 NEXT I

4540　　 COV(J1,J2)=DD*SS/(N1-M):COV(J2,J1)=COV(J1,J2)

4550　 NEXT J2

4560　 SD(J1)=SQR(COV(J1,J1))

4570　 SE(J1)=SD(J1)/SQR(N1)

4580 NEXT J1

4590 RETURN

4600 '

4610 ' *** DIFFERENCE APPROXIMATION ***

4620 '

4630 *APPROX:

4640 GOSUB *DEFINE.FORMULA

4650 F0=FNZ(X,Y)

4660 FX=FNX(X,Y)

4670 FY=FNY(X,Y)

4680 '

4690 FOR J=1 TO M

4700　 GOSUB *DEFINE.FORMULA2

4710　 F1(J)=FNJ(X,Y)

4720 NEXT J

4730 RETURN

4740 '

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

4760 '

4770 *INVERSE:' [GAUSS-JORDAN]

4780 FOR I=1 TO NZ

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

4800　 BZ(I,I)=1

4810 NEXT I

4820 FOR I=1 TO NZ

4830　 IZ=I

4840　 IF AZ(I,I)=0 THEN 4850 ELSE 4910

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

4860　 IF AZ(IZ,I)=0 THEN 4850

4870　 FOR J=1 TO NZ

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

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

4900　 NEXT J

4910　 AII=AZ(I,I)

4920　 FOR J=1 TO NZ

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

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

4950　 NEXT J

4960　 FOR K=1 TO NZ

4970　　 AKI=AZ(K,I):IF K=I THEN 5020

4980　　 FOR J=1 TO NZ

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

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

5010　　 NEXT J

5020　 NEXT K

5030 NEXT I

5040 RZ=0

5050 RETURN

5060 '

5070 ' *** SCATTER DIAGRAM ***

5080 '

5090 *SCATTER.DIAGRAM:

5100 CLS 3:WIDTH 80,25:COLOR 6

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

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

5130 IF DMAX<XMAX(1) THEN DMAX=XMAX(1)

5140 IF DMIN>XMIN(1) THEN DMIN=XMIN(1)

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

5160 '

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

5180 IF DMAX<XMAX(2) THEN DMAX=XMAX(2)

5190 IF DMIN>XMIN(2) THEN DMIN=XMIN(2)

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

5210 GOSUB *DRAW.AXIS

5220 GOSUB *DRAW.ELLIPSE:'GOSUB *GSAVE

5230 GOSUB *TEMP

5240 '

5250 CLS 3

5260 RETURN

5270 '

5280 ' *** 楕円 ***

5290 '

5300 *DRAW.ELLIPSE:

5310 'DRAW.COLOR=4

5320 WUO=UO:WUL=UL

5330 WVO=VO:WVL=VL

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

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

5360 '

5370 GOSUB *CONTOUR

5380 '

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

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

5410 RETURN

5420 '

5430 ' *** 輪郭 ***

5440 '

5450 *CONTOUR:

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

5470 AAA=P(1)

5480 BBB=P(2)

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

5500 '

5510 G=0

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

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

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

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

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

5570 NEXT ANGLE

5580 RETURN

5590 '

5600 ' *** 座標軸 ***

5610 '

5620 *DRAW.AXIS:

5630 CLS 3

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

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

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

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

5680 RR=3:GOSUB *DOT.PLOT

5690 GOSUB *TITLE.BACK2

5700 RETURN

5710 '

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

5730 '

5740 *X.SCALE2:

5750 CU=(SX2-SX1)/DU

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

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

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

5790 FOR K=0 TO DU

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

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

5820 NEXT K

5830 RETURN

5840 '

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

5860 '

5870 *X.REF2:

5880 REF.COLOR=6

5890 UW=UL-UO

5900 FOR K=0 TO DU STEP SU

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

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

5930　 FL=LEN(F$)

5940　 FOR L=1 TO FL

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

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

5970　 NEXT L

5980 NEXT K

5990 RETURN

6000 '

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

6020 '

6030 *Y.SCALE2:

6040 CV=(SY2-SY1)/DV

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

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

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

6080 FOR K=0 TO DV

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

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

6110 NEXT K

6120 RETURN

6130 '

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

6150 '

6160 *Y.REF2:

6170 REF.COLOR=6

6180 VW=VL-VO

6190 FOR K=0 TO DV STEP SV

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

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

6220　 FL=LEN(F$)

6230　 FOR L=1 TO FL

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

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

6260　 NEXT L

6270 NEXT K

6280 RETURN

6290 '

6300 ' *** プロット ***

6310 '

6320 *DOT.PLOT:

6330 WUO=UO:WUL=UL

6340 WVO=VO:WVL=VL

6350 WX=WUL-WUO

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

6370 AX=(SX2-SX1)/WX

6380 WY=WVL-WVO

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

6400 AY=(SY1-SY2)/WY

6410 '

6420 FOR I=1 TO N1

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

6440 XX=AX*X+BX

6450 YY=AY*Y+BY

6460 CIRCLE(XX,YY),RR,1

6470 PAINT(XX,YY),CLR,1

6480 CIRCLE(XX,YY),RR,CLR

6490 NEXT I

6500 RETURN

6510 '

6520 ' *** タイトル ***

6530 '

6540 *TITLE.BACK2:

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

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

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

6580 L1$="elliptic approximation by "

6590 'L2$="GAUSS-NEWTON METHOD"

6600 L$=L1$+L2$

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

6620 '

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

6640 FOR K=1 TO M

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

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

6670 NEXT K

6680 RETURN

6690 '

6700 *TEMP:

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

6720 RETURN

6730 '

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

6750 '

6760 *LINER.SCALE:

6770 DIM INTERVAL(10)

6780 RESTORE *REF.INTERVAL

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

6800 MAX.INTERVAL=10

6810 '

6820 ARGUMENT=DMAX-DMIN

6830 GOSUB *SCALE.RTN

6840 INDEX=CMP

6850 INTERVAL=INTERVAL(INDEX)*EP

6860 LO=INT(DMIN/INTERVAL)

6870 HI=-INT(-DMAX/INTERVAL)

6880 WHILE (HI-LO)>MAX.INTERVAL

6890　 INDEX=INDEX+1

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

6910　 INTERVAL=INTERVAL(INDEX)*EP

6920　 LO=INT(DMIN/INTERVAL)

6930　 HI=-INT(-DMAX/INTERVAL)

6940 WEND

6950 SMIN=LO*INTERVAL

6960 SMAX=HI*INTERVAL

6970 DIF=HI-LO

6980 STP=1

6990 ERASE INTERVAL

7000 RETURN

7010 '

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

7030 '

7040 *SCALE.RTN:

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

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

7070 'EP=10^RE

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

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

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

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

7120 RETURN

7130 '

7140 ' ** 参照間隔 **

7150 '

7160 *REF.INTERVAL:

7170 DATA .1,.2,.5,.5,.5,1,1,1,1,1

7180 '

7190 ' *** REPORT ***

7200 '

7210 *REPORT:

7220 CLS 3

7230 PRINT:PRINT

7240 PRINT "DATA No.=";N1

7250 PRINT

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

7270 FOR I=1 TO NV

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

7290 NEXT I

7300 PRINT:PRINT

7310 '

7320 PRINT "[ DEFINED EQUATION ]"

7330 PRINT "Y= ";FORM$

7340 '

7350 GOSUB *CEF.CONFIDENCE

7360 GOSUB *TEMP.STOP

7370 '

7380 CLS 3

7390 PRINT:PRINT "残差"

7400 PRINT " No."," OBSERVED X"," EXPECTED X"," OBSERVED Y"," EXPECTED Y"

7410 FOR I=1 TO N1

7420　 PRINT I,XX(I),X1(I),YY(I),Y1(I)

7430　 IF (I MOD 20)=0 THEN PRINT

7440 NEXT I

7450 PRINT

7460 GOSUB *TEMP.STOP

7470 RETURN

7480 '

7490 *CEF.CONFIDENCE:

7500 PRINT

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

7520 DF1=1:DF2=N1-M

7530 GOSUB *CEF.CONFIDENCE.SUB

7540 PRINT:PRINT

7550 RETURN

7560 '

7570 *CEF.CONFIDENCE.SUB:

7580 FOR J=1 TO M

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

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

7610　 IF SD(J)=0 THEN PRINT:GOTO 7640

7620　 GOSUB *CEF.ERROR

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

7640 NEXT J

7650 RETURN

7660 '

7670 *TEMP.STOP:

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

7690 WHILE INKEY$="":WEND

7700 RETURN

7710 '

7720 ' *** ERROR of COEFFICIENT ***

7730 '

7740 *CEF.ERROR:

7750 PP=.05 :GOSUB *F.PERCENT

7760 UUT=SQR(UUF)

7770 RETURN

7780 '

7790 ' ** NORMAL DISTRIBUTION **

7800 '

7810 *NORMAL.PERCENT:

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

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

7840 IF PP>.5 THEN UUN=-UUN

7850 RETURN

7860 '

7870 ' ** T-DISTRIBUTION **

7880 '

7890 *T.PERCENT:

7900 GOSUB *NORMAL.PERCENT

7910 UUN=ABS(UUN)

7920 UU2=UUN*UUN

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

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

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

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

7970 IF PP>.5 THEN UUT=-UUT

7980 RETURN

7990 '

8000 ' ** F-DISTRIBUTION **

8010 '

8020 *F.PERCENT:

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

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

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

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

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

8080 GOSUB *NORMAL.PERCENT

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

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

8110 RETURN

8120 '

8130 ' ** CHI-SQUARE DISTRIBUTION **

8140 '

8150 *CHI2.PERCENT:

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

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

8180 GOSUB *NORMAL.PERCENT

8190 UUX=2/(9*DF)

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

8210 UUX=UUX^3*DF

8220 RETURN

8230 '

8240 ' ** VRAM SAVE **

8250 '

8260 *GSAVE:

8270 FL$="b:circle"

8280 DEF SEG=&HA000

8290 BSAVE FL$+".CHR",&H0,&HFA0

8300 DEF SEG=&HA200

8310 BSAVE FL$+".ATR",&H0,&HFA0

8320 DEF SEG=&HA800

8330 BSAVE FL$+".BLU",&H0,&H7D00

8340 DEF SEG=&HB000

8350 BSAVE FL$+".RED",&H0,&H7D00

8360 DEF SEG=&HB800

8370 BSAVE FL$+".GRN",&H0,&H7D00

8380 RETURN

8390 '

8400 ' ** 係数変換 **

8410 '

8420 *P2A:

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

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

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

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

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

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

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

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

8510 RETURN

8520 '

8530 *A2P:

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

8550 DIF=A1-A3

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

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

8580 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)))

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

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

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

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

8630 RETURN

8640 '

8650 ' *** DEFINED EQUATION ***

8660 '

8670 *DEFINE.FORMULA3:

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

8690 RETURN

8700 '

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

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

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

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

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

8760 '

8770 ' ** 初期値 **

8780 '

8790 *INIT.VAL:

8800 GOSUB *LEAST.SQUARES

8810 GOSUB *A2P

8820 RETURN

8830 '

8840 ' ** XT*X **

8850 '

8860 *LEAST.SQUARES:

8870 FOR M1=1 TO M

8880　 FOR M2=M1 TO M

8890　 S=0

8900　　 FOR I=1 TO N1

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

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

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

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

8950　　 NEXT I

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

8970　 NEXT M2

8980 NEXT M1

8990 NZ=M

9000 '

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

9020 '

9030 GOSUB *INVERSE

9040 '

9050 ' ** XT*Y **

9060 '

9070 FOR J=1 TO M

9080　 S=0

9090　 FOR I=1 TO N1

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

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

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

9130　 NEXT I

9140　 W(J)=S

9150 NEXT J

9160 '

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

9180 '

9190 FOR J=1 TO M

9200　 S=0

9210　 FOR K=1 TO M

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

9230　 NEXT K

9240　 A(J)=S

9250 NEXT J

9260 RETURN

　

　このプログラムを実行すると，線形法で得られた初期値のほうが残差２乗和が小さく，かえって無駄な時間を費やしてしまった．しかし，微分に変えたため，計算自体は短時間で終了するようになった．

　　P(1)=　２２．８９±０．０４

　　P(2)=　１１．４１±０．０２

　　P(3)=−０．１５３±０．００２

　　P(4)=　１５．７２±０．０３

　　P(5)=　　７．９９±０．０２

　

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

　

［２］円近似プログラム（粟屋法）

　

　変更箇所だけ掲載すると

　

1010 ' **** circular approximation ****

1160 M=3:' No. of parameters

　

1420 FORM$="(X-P(1))^2+(Y-P(2))^2-P(3)^2"

1430 RETURN

1440 '

1450 ' *** DEFINED EQUATION ***

1460 '

1470 *DEFINE.FORMULA:

1480 DEF FNZ(X,Y)=(X-P(1))^2+(Y-P(2))^2-P(3)^2

1490 DEF FNX(X,Y)=(X-P(1))*2

1500 DEF FNY(X,Y)=(Y-P(2))*2

1510 RETURN

1520 '

1530 *DEFINE.FORMULA2:

1540 ON J GOSUB *J1,*J2,*J3

1550 RETURN

1560 '

1570 *J1:DEF FNJ(X,Y)=-(X-P(1))*2:RETURN

1580 *J2:DEF FNJ(X,Y)=-(Y-P(2))*2:RETURN

1590 *J3:DEF FNJ(X,Y)=-P(3)*2　　:RETURN

　

1960 'Z1(I)=1

1970　Z1(I)=X(1)*X(1)+X(2)*X(2)

　

4710 DMIN=P(1)-P(3):DMAX=P(1)+P(3)

4760 DMIN=P(2)-P(3):DMAX=P(2)+P(3)

　

5010 '

5020 ' *** 輪郭 ***

5030 '

5040 *CONTOUR:

5050 CCC=P(3)

5060 X0=P(1):Y0=P(2)

5070 '

5080 G=0

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

5100　 GX=CCC*COS(ANGLE)+X0

5110　 GY=CCC*SIN(ANGLE)+Y0

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

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

5140 NEXT ANGLE

5150 RETURN

　

6150 L1$="circular approximation by "

　

7960 '

7970 ' ** 係数変換 **

7980 '

7990 *P2A:

8000 P1=P(1):P2=P(2):P3=P(3)

8010 A(1)=P1*2

8020 A(2)=P2*2

8030 A(3)=P3*P3-P1*P1-P2*P2

8040 RETURN

8050 '

8060 *A2P:

8070 A1=A(1):A2=A(2):A3=A(3)

8080 P(1)=A1/2

8090 P(2)=A2/2

8100 P(3)=SQR(P(1)*P(1)+P(2)*P(2)+A3)

8110 RETURN

8120 '

8130 ' *** DEFINED EQUATION ***

8140 '

8150 *DEFINE.FORMULA3:

8160 ON L GOSUB *L1,*L2,*L3

8170 RETURN

8180 '

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

8200 *L2:DEF FNJ(X)=X(2):RETURN

8210 *L3:DEF FNJ(X)=1　 :RETURN

　

　最終収束値は，

　　ｘ0＝１２０．２７±３．２２

　　ｙ0＝　８５．６９±０．３８

　　ｒ ＝　９３．４５±３．２０

であって，「耕太郎」の計算結果

　　ｘ0＝１１８．７５±３．１３

　　ｙ0＝　８５．５８±０．３７

　　ｒ ＝　９１．９６±３．１１

と若干異なるが，それは微分を差分に変えたためと思われる．もちろん微分を使ったほうが正確である．

　

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