|
|
!'最小2乗法によりn次関数を求める
!' A + X*B + X^2*C + X^3*D + X^4*E = F(X)
!'Σ1*1 ΣX*1 ΣX^2*1 ΣX^3*1 ΣX^4*1 ΣF(X)*1
!'Σ1*X ΣX*X ΣX^2*X ΣX^3*X ΣX^4*X ΣF(X)*X
!'Σ1*X^2 ΣX*X^2 ΣX^2*X^2 ΣX^3*X^2 ΣX^4*X^2 ΣF(X)*X^2
!'Σ1*X^3 ΣX*X^3 ΣX^2*X^3 ΣX^3*X^3 ΣX^4*X^3 ΣF(X)*X^3
!'Σ1*X^4 ΣX*X^4 ΣX^2*X^4 ΣX^3*X^4 ΣX^4*X^4 ΣF(X)*X^4
LET N =10
RANDOMIZE
DIM X(N),Y(N)
LET A =RND
LET B =RND
LET C =RND
LET D =RND
FOR I =1 TO N
LET X(I)=I
LET Y(I)=A+B*X(I)+C*X(I)^2+D*X(I)^3 !'仮データ作成(3次式)
NEXT I
FOR I=1 TO N
PRINT "X=";X(I);"Y=";Y(I);"^Y=";FORECASTAUTO(X(I),N,Y,X)
NEXT I
END
EXTERNAL FUNCTION FORECASTAUTO(X,N,YA(),XA())
LET NLEVEL=7 !'当てはめ曲線の最高次数
DIM A(NLEVEL+1,NLEVEL+1),B(NLEVEL+1),R(NLEVEL+1),W(NLEVEL+1),WW(NLEVEL+1)
LET YN=MEAN(N,YA)
FOR I=1 TO N
LET SY=SY+(YA(I)-YN)^2
NEXT I
FOR I=0 TO NLEVEL
FOR J=0 TO NLEVEL
LET A(I+1,J+1)=A(I+1,J+1)+XA(I+1)^(I+J)
NEXT J
LET B(I+1)=B(I+1)+XA(I+1)^I*YA(I+1)
NEXT I
FOR LEVEL=2 TO NLEVEL !'2次式からNLEVEL次まで
MAT W=ZER
CALL CRAMER(LEVEL,A,B,W)
LET SE=0
FOR I=1 TO N
LET YY=W(LEVEL+1)
FOR J=LEVEL-1 TO 0 STEP -1
LET YY=YY*XA(I)+W(J+1)
NEXT J
LET SE=SE+(YA(I)-YY)^2
NEXT I
IF RR<1-SE/SY THEN !'決定係数より最良の次数を求める
LET RR=1-SE/SY
LET JISU=LEVEL
MAT WW=W !'係数の保存
END IF
NEXT LEVEL
LET YY=WW(JISU+1)
FOR J=JISU-1 TO 0 STEP -1
LET YY=YY*X+WW(J+1)
NEXT J
LET FORECASTAUTO=YY
END FUNCTION
EXTERNAL FUNCTION MEAN(N,X()) !'平均
FOR I=1 TO N
LET S=S+X(I)
NEXT I
LET MEAN=S/N
END FUNCTION
EXTERNAL SUB CRAMER (N, X(,), Y(), D()) !'クラーメル法
DIM A(N, N)
FOR I=1 TO N
FOR J=1 TO N
LET A(I,J)=X(I,J)
NEXT J
NEXT I
LET DD = DET(A)
IF DD = 0 THEN STOP !'ERROR
FOR K = 1 TO N
FOR I = 1 TO N
FOR J = 1 TO N
IF J = K THEN LET A(I, J) = Y(I) ELSE LET A(I, J) = X(I, J)
NEXT J
NEXT I
LET D(K) = DET(A) / DD
NEXT K
END SUB
|
|