|
|
最小化 (最小2乗法により回帰式を求める)
PUBLIC NUMERIC X(30),Y(30),N,H
DIM U(2,2),V(2),W(2)
RANDOMIZE
LET N=20 !'データー数
LET H=1/1024
LET AA=RND !'乱数でパラメータを決める
LET BB=RND
LET EPS=1E-6
FOR I=1 TO N !'テスト用データの作成
LET X(I)=I/5
LET Y(I)=FF(AA,BB,X(I))
NEXT I
LET A=1 !'初期値
LET B=1
DO !'ニュートン法
LET U(1,1)=DIFF(A,2,B,0)
LET U(1,2)=DIFF(A,1,B,1)
LET U(2,1)=DIFF(A,1,B,1)
LET U(2,2)=DIFF(A,0,B,2)
LET V(1)=-DIFF(A,1,B,0)
LET V(2)=-DIFF(A,0,B,1)
MAT U=INV(U)
MAT W=U*V
LET A=A+W(1)
LET B=B+W(2)
IF W(1)^2+W(2)^2<EPS THEN EXIT DO
LET L=L+1
IF L>100 THEN
PRINT "収束しません"
EXIT DO
END IF
LOOP
PRINT "AA=";AA,"BB=";BB
PRINT "A =";A," B =";B
LET P$="###.########"
FOR I=1 TO N
PRINT "X=";
PRINT USING P$:X(I);
PRINT " Y=";
PRINT USING P$:Y(I);
PRINT " ^Y=";
PRINT USING P$:FF(A,B,X(I))
NEXT I
END
EXTERNAL FUNCTION FF(A,B,X) !'回帰式
LET FF=A+B*X
!'LET FF=A*EXP(B*X)
!'LET FF=A+B*LOG(X)
!'LET FF=A+B*SQR(X)
!'LET FF=A*B^X
!'LET FF=A*X^B
!'LET FF=A+B/X
!'LET FF=1/(A+B*X)
!'LET FF=1/(A+B/X)
!'LET FF=1/(A+B*X^(1/3))
!'LET FF=1/(1+A*EXP(-B*X))
!'LET FF=X/(A+B*X^2)
!'LET FF=X/(A+B*LOG(X))
!'LET FF=SQR(X)/(A+B*LOG(X))
!'LET FF=LOG(X)/(A+B*SQR(X))
!'LET FF=EXP(X)/SQR(A+B*X^2)
!'LET FF=X/SQR(A+B*X)
!'LET FF=A^(B^X)
!'LET FF=1/(1+1/EXP(A+B*X))
END FUNCTION
EXTERNAL FUNCTION FUNC(A,B) !'最小化する関数(A=AA,B=BBの時、最小値 0)
FOR I=1 TO N
LET S=S+(Y(I)-FF(A,B,X(I)))^2
NEXT I
LET FUNC=S
END FUNCTION
EXTERNAL FUNCTION DIFF(X,M,Y,N)
IF M>0 THEN
LET DIFF=(-DIFF(X+2*H,M-1,Y,N)+8*DIFF(X+H,M-1,Y,N)-8*DIFF(X-H,M-1,Y,N)+DIFF(X-2*H,M-1,Y,N))/(12*H)
EXIT FUNCTION
END IF
IF N>0 THEN
LET DIFF=(-DIFF(X,M,Y+2*H,N-1)+8*DIFF(X,M,Y+H,N-1)-8*DIFF(X,M,Y-H,N-1)+DIFF(X,M,Y-2*H,N-1))/(12*H)
EXIT FUNCTION
END IF
IF M=0 OR N=0 THEN LET DIFF=FUNC(X,Y)
END FUNCTION
|
|