（その２）では

　　(1)確率構造の入った超楕円体の投影問題

　　(2)確率構造の入った超直方体の投影問題

を扱いました．

　

　データがｎ次元正規分布に従うとの仮定の下では，確率楕円を描くことは有効になりますが，ロボットアームなど工学上の問題を扱うときは，確率構造はかえってじゃまになります．たとえば，１００％のデータが入るような確率楕円では半径が∞になってしまうのですが，工学で求められている楕円はすべてのデータ点が含まれる有限の楕円であり，両者は大きく乖離してしまいます．

　

　そこで，工学への応用を考えて

　　(3)確率構造の入らない超楕円体・超直方体の投影問題

にも触れておくことにしました．この問題は(1)(2)より易しく，(1)(2)をきちんと押さえておけば簡単に解決できる問題です．

　

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

　

【１】確率構造の入らない投影問題

　

　ここでは，すべてのデータ点を含む最小体積の超楕円体，超直方体を描くことを考えます．その際，ω＝Σωi^2≦１をｎ次元単位超球の表現とすると，それに外接するｎ次元立方体は｜ωi｜≦１，内接する双対立方体はΣ｜ωi｜≦１で表されます．

　

［１］超楕円体の場合

　

　これまで説明してきたｎ次元楕円は

　　Σｈijｘiｘj＝ｘ’Ｈｘ＝ｃ^2

　　ｘ’＝［ｘ1，ｘ2，ｘ3，・・・，ｘn］

なのですが，座標変換

　［ｘ1，ｘ2，ｘ3，・・・，ｘn］’＝Ｕ［Ｘ1，Ｘ2，Ｘ3，・・・，Ｘn］’

　　Ｘ’＝［Ｘ1，Ｘ2，Ｘ3，・・・，Ｘn］

によって，

　　ｘ’Ｈｘ＝Ｘ’Ｕ’ＨＵＸ＝Ｘ’Ｕ~ＨＵＸ＝Ｘ’ΛＸ＝ｃ^2

　　Λ＝ｄｉａｇ（λ1，・・・，λn）

のように，標準形に変換することができます．

　

　変換後でも右辺の値は変わりません．右辺の値がｃ^2のとき，

　　ｘ’Ｈｘ＝ｃ^2

は変換後も

　　Ｘ’ΛＸ＝ｃ^2

となるというわけです．

　

　したがって，超楕円体の場合，変換後の値

　　Ｘ’ΛＸ＝ｃ^2

を用いずとも，変換前（＝実際のデータ）の値から

　　ｘ’Ｈｘ＝ｃ^2

により，ｃ^2の最大値を求めればよいことになります．なお，ｃ^2は固有値で重みづけした標準化平方ユークリッド距離に相当します．

　

［２］超直方体の場合

　

　超楕円体に外接する超直方体はｃ^2を求めることによって与えられます．これは体積最小の平行多面体ですが，ここで求める超直方体は超楕円体に外接するものではなく，すべてのデータ点を含む最小体積の超直方体です．

　

　その場合，超直方体の稜線の長さを求めるには座標変換後の値が必要となるので，超楕円体に較べて計算式は面倒になります．とはいっても，変換行列Ｕは直交行列ですから，

　　Ｕ’＝Ｕ~

なる性質があり，

　　ｘ＝ＵＸ　　→　　Ｘ＝Ｕ~ｘ＝Ｕ’ｘ

のように，逆行列を求めなければならないところが，転置行列で済ませられます．

　

　このようにして，超直方体の稜線の長さを，座標変換後の絶対値の最大値

　　｜Ｘi｜　　（i=1〜n）

から求めることができます．

　

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

　

【２】実行例

（超楕円体に外接する超直方体に設定しています）

（隠線処理はしておりません）

　

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

　

【３】プログラム

　

1000 '

1010 ' **** mechanical ellipsoid ****

1020 '　　　　　　　　　　　　2002/11/22 (C) ｻﾄｳ ｲｸﾛｳ

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

1040 CLS 3:WIDTH 80,25:COLOR 6

1050 GOSUB *INITIALIZE

1060 GOSUB *DATAIN

1070 GOSUB *HESSIAN

1080 GOSUB *DISTANCE:' [mechanical]

1090　GOSUB *ELLIPSOID

1190 '

1200 CLS 3

1210 END

1220 '

1230 ' *** 初期設定 ***

1240 '

1250 *INITIALIZE:

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

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

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

1290 SCALE.KIND=1:PLOT=1

1300 PI=3.14159

1310 JX=.709:JY=1

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

1330 '

1340 M=3:N1=51

1350 'M=2:N1=49

1360 PP=.05

1370 '

1380 MMAX=10:NMAX=10

1390 DIM V(MMAX)

1400 DIM JACOBI(MMAX,NMAX)

1410 DIM AZ(MMAX,MMAX),BZ(MMAX,MMAX)

1420 DIM HESSE(MMAX,MMAX)

1430 DIM COV(MMAX,MMAX)

1440 DIM VO(MMAX),VL(MMAX)

1450 DIM VD(MMAX),VS(MMAX)

1460 DIM UNIT$(MMAX)

1470 '

1480 FOR I=1 TO 3

1490　 VO(I)=0:VL(I)=5

1500　 VD(I)=5:VS(I)=1

1510　 UNIT$(I)="m/s"

1520 NEXT I

1530 FOR I=4 TO 6

1540　 VO(I)=0:VL(I)=2

1550　 VD(I)=4:VS(I)=1

1560　 UNIT$(I)="rad/s"

1570 NEXT I

1580 '

1590 DIM X(MMAX),Y(MMAX)

1600 DIM XY(MMAX,MMAX)

1610 DIM X1(MMAX),X2(MMAX)

1620 DIM AVE(MMAX),DEV(MMAX)

1630 '

1640 DIM EDGE(MMAX)

1650 DIM S(MMAX),S0(MMAX)

1660 DIM GXS(MMAX),GYS(MMAX)

1670 DIM GXE(MMAX),GYE(MMAX)

1680 DIM GZ(MMAX)

1690 DIM D1(MMAX),D2(MMAX)

1700 RETURN

1710 '

1720 ' *** データの読み込み ***

1730 '

1740 *DATAIN:

1750 CLS 3

1760 FOR I=1 TO M

1770　 X1(I)=0:X2(I)=0

1780　 FOR J=1 TO M

1790　　 XY(I,J)=0

1800　 NEXT J

1810 NEXT I

1820 '

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

1840 RESTORE *V

1850 FOR I=1 TO N1

1860　 PRINT "*";

1870　 FOR J=1 TO M:READ X(J):NEXT J

1880　 FOR J=1 TO M

1890　　 X=X(J):GOSUB *CALC1

1900　　 FOR K=J TO M

1910　　　 Y=X(K):GOSUB *CALC2

1920　　 NEXT K

1930　 NEXT J

1940 NEXT I

1950 PRINT

1960 FOR I=1 TO M

1970　 AVE(I)=X1(I)/N1

1980　 S=X2(I)-X1(I)*X1(I)/N1

1990　 DEV(I)=SQR(S/(N1-1))

2000　 FOR J=I TO M

2010　　 S=XY(I,J)-X1(I)*X1(J)/N1

2020　　 COV(I,J)=S/(N1-1)

2030　　 COV(J,I)=COV(I,J)

2040　 NEXT J

2050 NEXT I

2060 RETURN

2070 '

2080 ' ** CALCULATE **

2090 '

2100 *CALC1:

2110 X1(J)=X1(J)+X

2120 X2(J)=X2(J)+X*X

2130 RETURN

2140 '

2150 *CALC2:

2160 XY(J,K)=XY(J,K)+X*Y

2170 RETURN

2180 '

2190 ' *** ヘシアンと共分散行列 ***

2200 '

2210 *HESSIAN:

2220 FOR M1=1 TO M

2230　 FOR M2=M1 TO M

2240　　 S=COV(M1,M2)

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

2260　 NEXT M2

2270 NEXT M1

2280 NZ=M

2290 '

2300 GOSUB *INVERSE

2310 '

2320 FOR I=1 TO M

2330　 FOR J=1 TO M:HESSE(I,J)=BZ(I,J):NEXT J

2340 NEXT I

2350 RETURN

2360 '

2370 ' ** 逆行列 ( AZ==>BZ ) **

2380 '

2390 *INVERSE:' [GAUSS-JORDAN]

2400 FOR I=1 TO NZ

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

2420　 BZ(I,I)=1

2430 NEXT I

2440 FOR I=1 TO NZ

2450　 IZ=I

2460　 IF AZ(I,I)=0 THEN 2470 ELSE 2530

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

2480　 IF AZ(IZ,I)=0 THEN 2470

2490　 FOR J=1 TO NZ

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

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

2520　 NEXT J

2530　 AII=AZ(I,I)

2540　 FOR J=1 TO NZ

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

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

2570　 NEXT J

2580　 FOR K=1 TO NZ

2590　　 AKI=AZ(K,I):IF K=I THEN 2640

2600　　 FOR J=1 TO NZ

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

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

2630　　 NEXT J

2640　 NEXT K

2650 NEXT I

2660 RZ=0

2670 RETURN

2680 '

2690 ' *** 標準化平方距離 ***

2700 '

2710 *DISTANCE:

2720 RESTORE *V

2730 MAX.D2=0

2740 FOR IS=1 TO N1

2750　 D2=0

2760　 FOR JS=1 TO M:READ X(JS):NEXT JS

2770　 FOR JS=1 TO M

2780　　 FOR KS=1 TO M

2790　　　 D2=D2+(X(JS)-AVE(JS))*(X(KS)-AVE(KS))*HESSE(JS,KS)

2800　　 NEXT KS

2810　 NEXT JS

2820　 IF MAX.D2<D2 THEN MAX.D2=D2

2830 NEXT IS

2840 RETURN

2850 '

2860 ' *** 標準化平方距離 ***

2870 '

2880 *DISTANCE.2:

2890 RESTORE *V

2900 MAX.D2=0

2910 FOR IS=1 TO N1

2920　 D2=0

2930　 FOR JS=1 TO M:READ X(JS):NEXT JS

2940　 FOR JS=1 TO M

2950　　 FOR KS=1 TO M:D1(KS)=0:NEXT KS

2960　　 FOR KS=1 TO M

2970　　　 D1(JS)=D1(JS)+BZ(KS,JS)*(X(KS)-AVE(KS))

2980　　 NEXT KS

2990　　 D2=D2+D1(JS)*D1(JS)/AZ(JS,JS)

3000　 NEXT JS

3010　 FOR JS=1 TO M

3020　　 IF MAX.D2<ABS(D2) THEN MAX.D2=ABS(D2)

3030　 NEXT JS

3040 NEXT IS

3050 RETURN

3060 '

3070 ' *** 稜線 ***

3080 '

3090 *DISTANCE.RIDGE:

3100 RESTORE *V

3110 FOR JS=1 TO M:MAX.D2(JS)=0:NEXT JS

3120 FOR IS=1 TO N1

3130　 FOR JS=1 TO M:D1(JS)=0:NEXT JS

3140　 FOR JS=1 TO M:READ X(JS):NEXT JS

3150　 FOR JS=1 TO M

3160　　 FOR KS=1 TO M

3170　　　 D1(JS)=D1(JS)+BZ(KS,JS)*(X(KS)-AVE(KS))

3180　　 NEXT KS

3190　 NEXT JS

3200　 FOR JS=1 TO M

3210　　 IF D2(JS)<ABS(D1(JS)) THEN D2(JS)=ABS(D1(JS))

3220　 NEXT JS

3230 NEXT IS

3240 RETURN

3250 '

3260 ' *** 楕円体 ***

3270 '

3280 *ELLIPSOID:

3290 IF M=1 THEN RETURN

3300 '

3310 FOR I=1 TO M

3320　 FOR J=1 TO M:AZ(I,J)=COV(I,J):NEXT J

3330 NEXT I

3340 NZ=M:GOSUB *JACOBI

3350 'GOSUB *DISTANCE.2:' [mechanical]

3360 '

3370 FOR I=1 TO M-1

3380　 FOR J=I+1 TO M

3390　　 UO=VO(I):UL=VL(I)

3400　　 VO=VO(J):VL=VL(J)

3410　　 DU=VD(I):SU=VS(I)

3420　　 DV=VD(J):SV=VS(J)

3430 '

3440　　 GOSUB *DRAW.AXIS

3450　　 GOSUB *DRAW.ELLIPSE

3460　　 GOSUB *TEMP

3470　 NEXT J

3480 NEXT I

3490 RETURN

3500 '

3510 ' ** 固有値 ( AZ==>BZ ) **

3520 '

3530 *JACOBI:' [ EIGEN-VECTOR: BZ(J,I) ]

3540 EPSL=.00001

3550 FOR I=1 TO NZ

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

3570　 BZ(I,I)=1

3580 NEXT I

3590 '

3600 CMAX=0

3610 FOR I=1 TO NZ

3620　 FOR J=1 TO NZ

3630　　 IF CMAX<ABS(AZ(I,J)) THEN CMAX=ABS(AZ(I,J))

3640　 NEXT J

3650 NEXT I

3660 EPS=CMAX*EPSL

3670 '

3680 LMAX=100

3690 LOOP=0

3700 SW=0

3710 WHILE LOOP<LMAX AND SW=0

3720 '

3730 FOR L=1 TO NZ-1

3740 FOR I=1 TO NZ-1

3750　 J=I+L:IF J>NZ THEN 3940

3760　 IF ABS(AZ(I,J))<EPS THEN 3930

3770　 DIF=AZ(I,I)-AZ(J,J)

3780　 IF DIF=0 THEN T=PI/4 ELSE T=1/2*ATN(2*AZ(I,J)/DIF)

3790　 C=COS(T):S=SIN(T)

3800　 FOR K=1 TO NZ

3810　　 A1=AZ(I,K):A2=AZ(J,K)

3820　　 AZ(I,K)= C*A1+S*A2

3830　　 AZ(J,K)=-S*A1+C*A2

3840　 NEXT K

3850　 FOR K=1 TO NZ

3860　　 A1=AZ(K,I):A2=AZ(K,J)

3870　　 AZ(K,I)= C*A1+S*A2

3880　　 AZ(K,J)=-S*A1+C*A2

3890　　 B1=BZ(K,I):B2=BZ(K,J)

3900　　 BZ(K,I)= C*B1+S*B2

3910　　 BZ(K,J)=-S*B1+C*B2

3920　 NEXT K

3930 NEXT I

3940 NEXT L

3950 '

3960 SW=1

3970 FOR I=1 TO NZ-1

3980　 FOR J=I+1 TO NZ

3990　　 IF ABS(AZ(I,J))>EPS THEN SW=0

4000　 NEXT J

4010 NEXT I

4020 LOOP=LOOP+1

4030 '

4040 WEND

4050 RZ=0

4060 RETURN

4070 '

4080 ' *** 正射影 ***

4090 '

4100 *DRAW.ELLIPSE:

4110 'DRAW.COLOR=4

4120 WUO=UO-UL:WUL=UO+UL

4130 WVO=VO-VL:WVL=VO+VL

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

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

4160 '

4170 'GOSUB *POLAR

4180 'GOSUB *EQUATOR

4190 'GOSUB *SECTION

4200 GOSUB *CONTOUR

4210　GOSUB *HYPERCUBE

4220 '

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

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

4250 RETURN

4260 '

4270 ' *** 輪郭 ***

4280 '

4290 *CONTOUR:

4300 DIF=COV(I,I)-COV(J,J)

4310 IF DIF=0 THEN TH=PI/4 ELSE TH=1/2*ATN(2*COV(I,J)/DIF)

4320 C=COS(TH):S=SIN(TH)

4330 '

4340 LAMBDA1=COV(I,I)*C*C+2*COV(I,J)*C*S+COV(J,J)*S*S

4350 LAMBDA2=COV(I,I)*S*S-2*COV(I,J)*C*S+COV(J,J)*C*C

4360 T=MAX.D2: ' [mechanical]

4370 GOTO 4410

4380 T=1:' [non-stochastic]

4390 T=-2*LOG(PP): ' [stochastic]

4400 'DF=2:GOSUB *CHI2.PERCENT:T=UUX

4410 '

4420 AAA=SQR(T*LAMBDA1)

4430 BBB=SQR(T*LAMBDA2)

4440 '

4450 G=0

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

4470　 GX=AAA*COS(ANGLE)*C-BBB*SIN(ANGLE)*S+AVE(I)

4480　 GY=AAA*COS(ANGLE)*S+BBB*SIN(ANGLE)*C+AVE(J)

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

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

4510 NEXT ANGLE

4520 RETURN

4530 '

4540 ' *** 超直方体 ***

4550 '

4560 *HYPERCUBE:

4570 'GOSUB *DISTANCE.RIDGE:' [mechanical]

4580 'FOR IS=1 TO M:EDGE(IS)=D2(IS):NEXT IS

4590 T=MAX.D2:' [external]

4600 FOR IS=1 TO M:EDGE(IS)=SQR(T*AZ(IS,IS)):NEXT IS

4610 GOTO 4680

4620 '

4630 T=1:' [non-stochastic]

4640 FOR IS=1 TO M:EDGE(IS)=SQR(T*AZ(IS,IS)):NEXT IS

4650 ' [stochastic]

4660 PP0=PP:PP=1-(1-PP)^(1/M):GOSUB *NORMAL.PERCENT:PP=PP0

4670 FOR IS=1 TO M:EDGE(IS)=EDGE(IS)*UUN:NEXT IS

4680 '

4690 FOR IS=0 TO 2^M-1

4700　 FOR JS=1 TO M:S(JS)=-1:NEXT JS

4710　 FOR JS=1 TO M

4720　　 IF (IS AND 2^(JS-1))=0 THEN S(JS)=1

4730　 NEXT JS

4740　'GOSUB *VERTEX:GXS=GX:GYS=GY

4750　 GOSUB *VERTEX:GXS=GX:GYS=GY

4760　 '

4770　 FOR JS=1 TO M:S0(JS)=S(JS):NEXT JS

4780　 FOR JS=1 TO M

4790　　'FOR KS=1 TO M:S(KS)=S0(KS):NEXT KS

4800　　'S(JS)=-S(JS)

4810　　'GOSUB *VERTEX:GXE=GX:GYE=GY

4820　　'LINE(GXS,-GYS)-(GXE,-GYE),7

4830　　 FOR KS=1 TO M:S(KS)=S0(KS):NEXT KS

4840　　 S(JS)=-S(JS)

4850　　 GOSUB *VERTEX:GXE=GX:GYE=GY

4860　　 LINE(GXS,-GYS)-(GXE,-GYE),7

4870　　'GOSUB *POLYGON.LINE

4880　 NEXT JS

4890 NEXT IS

4900 RETURN

4910 '

4920 ' *** 頂点 ***

4930 '

4940 *VERTEX:

4950 GX=0:GY=0

4960 FOR LS=1 TO M

4970　 GX=GX+BZ(I,LS)*EDGE(LS)*S(LS)

4980　 GY=GY+BZ(J,LS)*EDGE(LS)*S(LS)

4990 NEXT LS

5000 GX=GX+AVE(I):GY=GY+AVE(J)

5010 RETURN

5020 '

5030 ' *** 輪郭抽出 ***

5040 '

5050 *POLYGON.LINE:

5060 EPS=.05

5070 FOR KS=1 TO M

5080　 FOR LS=1 TO M:S(LS)=S0(LS):NEXT LS

5090　 IF KS=JS THEN 5110

5100　 S(KS)=S(KS)*(1-EPS)

5110　 GOSUB *VERTEX:GXS(KS)=GX:GYS(KS)=GY

5120　 '

5130　 FOR LS=1 TO M:S(LS)=S0(LS):NEXT LS

5140　 S(JS)=-S(JS)

5150　 IF KS=JS THEN 5170

5160　 S(KS)=S(KS)*(1-EPS)

5170　 GOSUB *VERTEX:GXE(KS)=GX:GYE(KS)=GY

5180 NEXT KS

5190 '

5200 X=GXS(JS):Y=GYS(JS)

5210 FOR KS=1 TO M

5220　 GZ(KS)=(GYE(KS)-GYS(KS))*(X-GXS(KS))-(GXE(KS)-GXS(KS))*(Y-GYS(KS))

5230 NEXT KS

5240 '

5250 ZMAX=GZ(JS):ZMIN=GZ(JS)

5260 FOR KS=1 TO M

5270　 IF KS=JS THEN 5300

5280　 IF ZMAX<GZ(KS) THEN ZMAX=GZ(KS)

5290　 IF ZMIN>GZ(KS) THEN ZMIN=GZ(KS)

5300 NEXT KS

5310 '

5320 SW=0

5330 IF (ZMAX-GZ(JS))*(ZMIN-GZ(JS))=0 THEN SW=1

5340 GXS=GXS(JS):GYS=GYS(JS)

5350 GXE=GXE(JS):GYE=GYE(JS)

5360 IF SW=1 THEN LINE(GXS,-GYS)-(GXE,-GYE),7

5370 RETURN

5380 '

5390 ' *** 座標軸 ***

5400 '

5410 *DRAW.AXIS:

5420 CLS 3

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

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

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

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

5470 RR=3:GOSUB *DOT.PLOT

5480 GOSUB *TITLE.BACK2

5490 RETURN

5500 '

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

5520 '

5530 *X.SCALE2:

5540 CU=(SX2-SX1)/DU/2

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

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

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

5580 FOR K=0 TO DU*2

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

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

5610 NEXT K

5620 'GOSUB *X.REF2

5630 RETURN

5640 '

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

5660 '

5670 *X.REF2:

5680 REF.COLOR=6

5690 UW=UL-UO

5700 FOR K=0 TO DU*2 STEP SU

5710　 F=UW*(K/DU-1):F$=STR$(F)

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

5730　 FL=LEN(F$)

5740　 FOR L=1 TO FL

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

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

5770　 NEXT L

5780 NEXT K

5790 RETURN

5800 '

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

5820 '

5830 *Y.SCALE2:

5840 CV=(SY2-SY1)/DV/2

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

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

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

5880 FOR K=0 TO DV*2

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

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

5910 NEXT K

5920 'GOSUB *Y.REF2

5930 RETURN

5940 '

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

5960 '

5970 *Y.REF2:

5980 REF.COLOR=6

5990 VW=VL-VO

6000 FOR K=0 TO DV*2 STEP SV

6010　 F=VW*(K/DV-1):F$=STR$(F)

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

6030　 FL=LEN(F$)

6040　 FOR L=1 TO FL

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

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

6070　 NEXT L

6080 NEXT K

6090 RETURN

6100 '

6110 ' *** プロット ***

6120 '

6130 *DOT.PLOT:

6140 WUO=UO-UL:WUL=UO+UL

6150 WVO=VO-VL:WVL=VO+VL

6160 WX=WUL-WUO

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

6180 AX=(SX2-SX1)/WX

6190 WY=WVL-WVO

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

6210 AY=(SY1-SY2)/WY

6220 '

6230 RESTORE *V

6240 FOR K=1 TO N1

6250　 FOR L=1 TO M:READ V(L):NEXT L

6260　 XX=AX*V(I)+BX

6270　 YY=AY*V(J)+BY

6280　 CIRCLE(XX,YY),RR,1

6290　 PAINT(XX,YY),CLR,1

6300　 CIRCLE(XX,YY),RR,CLR

6310 NEXT K

6320 RETURN

6330 '

6340 ' *** タイトル ***

6350 '

6360 *TITLE.BACK2:

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

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

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

6400 L1$="ellipsoid projection of "

6410 L2$="v"+MID$(STR$(I),2)+"-"

6420 L3$="v"+MID$(STR$(J),2)+" plane"

6430 L$=L1$+L2$+L3$

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

6450 L$=UNIT$(J)

6460 IF LEN(L$)=0 THEN 6480

6470 LOCATE 7,LY-1:PRINT "("+L$+")"

6480 L$=UNIT$(I)

6490 IF LEN(L$)=0 THEN 6510

6500 LOCATE LX+3,21:PRINT "("+L$+")"

6510 '

6520 'LOCATE LX+3,LY:PRINT "v components: "

6530 FOR K=1 TO M

6540　 V$="v"+MID$(STR$(K),2)+"="

6550　'LOCATE LX+3,LY+K+1:PRINT V$;:PRINT V(K)

6560 NEXT K

6570 RETURN

6580 '

6590 *TEMP:

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

6610 RETURN

10180 '

10190 *V:

10200 DATA 0.2537,-0.1464,0

10210 DATA -0.75,-1.299,1.5

10220 DATA 0.1036,0.1794,0.5

10230 DATA -0.5,-0.866,1

10240 DATA -1.1036,-1.9114,1.5

10250 DATA -1.7071,-2.9568,2

10260 DATA 0.3536,0.6124,0

10270 DATA -0.25,-0.433,0.5

10280 DATA -0.8536,-1.4784,1

10290 DATA -1.4571,-2.5238,1.5

10300 DATA 1.2071,2.0908,-1

10310 DATA 0.6036,1.0454,-0.5

10320 DATA 0,0,0

10330 DATA -0.6036,-1.0454,0.5

10340 DATA -1.2071,-2.0908,1

10350 DATA -1.8107,-3.1362,1.5

10360 DATA 0.8536,1.4784,-1

10370 DATA 0.25,0.433,-0.5

10380 DATA -0.3536,-0.6124,0

10390 DATA -0.9571,-1.6578,0.5

10400 DATA 0.5,0.866,-1

10410 DATA -0.1501,0.3258,0.5

10420 DATA -0.7537,-0.7196,1

10430 DATA -1.3572,-1.765,1.5

10440 DATA -1.9608,-2.8103,2

10450 DATA 0.0999,0.7588,0

10460 DATA -0.5037,-0.2866,0.5

10470 DATA -1.1072,-1.332,1

10480 DATA -1.7108,-2.3773,1.5

10490 DATA 0.3499,1.1918,-0.5

10500 DATA -0.2537,0.1464,0

10510 DATA -0.8572,-0.8989,0.5

10520 DATA -1.4608,-1.9443,1

10530 DATA 0.5999,1.6248,-1

10540 DATA -0.0037,0.5795,-0.5

10550 DATA -0.6072,-0.4659,0

10560 DATA -1.2108,-1.5113,0.5

10570 DATA -1.0073,-0.5731,1

10580 DATA -1.6109,-1.6185,1.5

10590 DATA -0.1538,0.9053,0

10600 DATA -0.7573,-0.1401,0.5

10610 DATA -1.3609,-1.1855,1

10620 DATA -1.9644,-2.2309,1.5

10630 DATA 0.0962,1.3383,-0.5

10640 DATA -0.5073,0.2929,0

10650 DATA -1.1109,-0.7525,0.5

10660 DATA -1.7144,-1.7979,1

10670 DATA 0.3462,1.7713,-1

10680 DATA -0.2573,0.7259,-0.5

10690 DATA -0.8609,-0.3195,0

10700 DATA -0.761,0.4393,0

10710 '

10720 '*V:

10730 DATA 0.3087,-0.1783

10740 DATA 0.3326,0.1544

10750 DATA 0.1544,-0.0891

10760 DATA -0.0239,-0.3326

10770 DATA -0.2021,-0.5761

10780 DATA -0.3804,-0.8196

10790 DATA -0.5586,-1.0631

10800 DATA 0.3565,0.487

10810 DATA 0.1783,0.2435

10820 DATA 0,0

10830 DATA -0.1783,-0.2435

10840 DATA -0.3565,-0.487

10850 DATA -0.5348,-0.7305

10860 DATA -0.713,-0.974

10870 DATA -0.8913,-1.2175

10880 DATA 0.0239,0.3326

10890 DATA -0.1544,0.0891

10900 DATA -0.3326,-0.1544

10910 DATA -0.5109,-0.3979

10920 DATA -0.6891,-0.6414

10930 DATA -0.8674,-0.8849

10940 DATA -0.1305,0.4217

10950 DATA -0.3087,0.1783

10960 DATA -0.487,-0.0652

10970 DATA -0.6652,-0.3087

10980 DATA -0.8435,-0.5522

10990 DATA -1.0217,-0.7957

11000 DATA -0.2849,0.5109

11010 DATA -0.4631,0.2674

11020 DATA -0.6414,0.0239

11030 DATA -0.8196,-0.2196

11040 DATA -0.9979,-0.4631

11050 DATA -1.1761,-0.7066

11060 DATA -0.4392,0.6

11070 DATA -0.6175,0.3565

11080 DATA -0.7957,0.113

11090 DATA -0.974,-0.1305

11100 DATA -1.1522,-0.374

11110 DATA -1.3305,-0.6175

11120 DATA -0.5936,0.6891

11130 DATA -0.7719,0.4456

11140 DATA -0.9501,0.2021

11150 DATA -1.1284,-0.0414

11160 DATA -1.3066,-0.2849

11170 DATA -0.748,0.7783

11180 DATA -0.9262,0.5348

11190 DATA -1.1045,0.2913

11200 DATA -1.2827,0.0478

11210 DATA -1.0806,0.6239

11220 '

11230 ' ** NORMAL DISTRIBUTION **

11240 '

11250 *NORMAL.PERCENT:

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

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

11280 IF PP>.5 THEN UUN=-UUN

11290 RETURN

11300 '

11310 ' ** T-DISTRIBUTION **

11320 '

11330 *T.PERCENT:

11340 GOSUB *NORMAL.PERCENT

11350 UUN=ABS(UUN)

11360 UU2=UUN*UUN

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

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

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

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

11410 IF PP>.5 THEN UUT=-UUT

11420 RETURN

11430 '

11440 ' ** F-DISTRIBUTION **

11450 '

11460 *F.PERCENT:

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

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

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

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

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

11520 GOSUB *NORMAL.PERCENT

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

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

11550 RETURN

11560 '

11570 ' ** CHI-SQUARE DISTRIBUTION **

11580 '

11590 *CHI2.PERCENT:

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

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

11620 GOSUB *NORMAL.PERCENT

11630 UUX=2/(9*DF)

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

11650 UUX=UUX^3*DF

11660 RETURN

11670 '

11680 ' ** NORMAL DISTRIBUTION **

11690 '

11700 *NORMAL.PROB:

11710 AA1=.196854

11720 AA2=.115194

11730 AA3=.000344

11740 AA4=.019527

11750 W=ABS(UUN)

11760 P0=1+W*(AA1+W*(AA2+W*(AA3+W*AA4)))

11770 P0=P0^4

11780 P0=1-.5/P0

11790 P0=.5+(P0-.5)*SGN(UUN)

11800 PP=1-P0

11810 RETURN

11820 '

11830 ' ** T-DISTRIBUTION **

11840 '

11850 *T.PROB:

11860 W1=(1-2/(3+8*DF))

11870 W2=DF*LOG(UUT*UUT/DF+1)

11880 UUN=W1*SQR(W2)

11890 GOSUB *NORMAL.PROB

11900 PP=1-P0

11910 RETURN

11920 '

11930 ' ** F-DISTRIBUTION **

11940 '

11950 *F.PROB:

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

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

11980 FSW=0

11990 IF UUF=0 THEN P0=0:GOTO 12100

12000 IF UUF<1 THEN UUF=1/UUF:SWAP DF1,DF2:FSW=1

12010 AA1=2/(9*DF1)

12020 AA2=2/(9*DF2)

12030 W=UUF^(1/3)

12040 W1=W+AA1-W*AA2-1

12050 W2=AA2*W*W+AA1

12060 UUN=W1/SQR(W2)

12070 IF DF2<=3 THEN UUN=UUN*(1+.08*(UUN^4)/(DF2^3))

12080 GOSUB *NORMAL.PROB

12090 IF FSW=1 THEN P0=1-P0:SWAP DF1,DF2:FSW=0

12100 PP=1-P0

12110 RETURN

12120 '

12130 ' ** CHI-SQUARE DISTRIBUTION **

12140 '

12150 *CHI2.PROB:

12160 IF DF=2 THEN PP=EXP(-UUX/2):RETURN

12170 AA1=2/(9*DF)

12180 AA2=UUX/DF

12190 W=AA2^(1/3)

12200 W1=W+AA1-1

12210 UUN=W1/SQR(AA1)

12220 GOSUB *NORMAL.PROB

12230 PP=1-P0

12240 RETURN

　

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