|
|
> No.4858[元記事へ]
大熊 正さんへのお返事です。
この掲示板の常連投稿者である山中和義さんの個人サイト(現在は閉鎖)で公開されていたプログラムです。
!重回帰分析 - 重回帰式 Y=a+b*X1+c*X2+ … +d*Xp
!●重回帰式 Y=2.504+0.8161*X1-0.2749*X2+1.055*X3
PUBLIC NUMERIC P,N
LET P=3 !変数の数
LET N=8 !データの数
DATA 8,4,8,18 !X1,X2,X3,Y
DATA 7,7,7,12
DATA 5,8,9,14
DATA 4,3,3,6
DATA 6,8,8,12
DATA 2,5,3,8
DATA 3,6,6,10
DATA 9,9,7,16
DIM X(N,P),Y(N) !データを読み込む
FOR i=1 TO N
FOR j=1 TO P
READ X(i,j) ![X1,X2,X3]
NEXT j
READ Y(i) !Y
NEXT i
DIM C(P+1)
CALL LINEST(Y,X, C) !係数
MAT PRINT C;
PRINT "決定係数の２乗=";RSQ(Y,X) !0.83405
DIM t(P) !X1,X2,X3
DATA 5,8,9
MAT READ t
PRINT FORECAST(t,Y,X) !予測値
PRINT "予測値の標準誤差=";STEYX(Y,X)
!●ロジスティック回帰モデル Y = 1 / ( 1 + EXP(a+b*X) )
LET P=1
LET N=8
DATA 5.42, 0.05 !X,Y
DATA 5.61, 0.18
DATA 5.78, 0.25
DATA 5.95, 0.47
DATA 6.12, 0.82
DATA 6.28, 0.89
DATA 6.43, 0.98
DATA 6.58, 1.00
DIM XX(N,P),YY(N) !データを読み込む
FOR i=1 TO N
READ XX(i,1) ![X1]
READ YY(i) !Y
NEXT i
!線形回帰式に変形する。LOG(1/Y-1)=a+b*X → Y'=a+b*X → Y=1/(1+EXP(Y'))
DIM LogYY(N)
FOR i=1 TO N
IF YY(i)=1 THEN
LET LogYY(i)=-5 !※
ELSE
LET LogYY(i)=LOG(1/YY(i)-1)
END IF
NEXT i
CALL LINEST(LogYY,XX, C) !係数
MAT PRINT C;
PRINT "決定係数の２乗=";RSQ(LogYY,XX)
SET WINDOW 5,7,-0.2,1.2 !グラフを描いてみる
DRAW grid(0.5,0.25)
FOR i=1 TO N !データ
PLOT POINTS: XX(i,1),YY(i)
NEXT i
FOR i=5 TO 7 STEP 0.1 !近似式
LET t(1)=i
PLOT LINES: i,1/(1+EXP( FORECAST(t,LogYY,XX) ));
NEXT i
PLOT LINES
END
!最小２乗法による
!
!残差平方和 S=Σ[i=1,n]{(Y[i]-(a+b*X1+c*X2+ … +d*Xp))^2} を最小にする係数を求める。
!偏微分（∂S/∂a、∂S/∂bなど）して、連立方程式を得る。
! ┌ Σ 1*1 Σ 1*X1 Σ 1*X2 … Σ 1*Xp ┐┌ a ┐=┌ Σ 1*Y ┐
! │ ΣX1*1 ΣX1*X1 ΣX1*X2 … ΣX1*Xp ││ b │ │ ΣX1*Y │
! │ ΣX2*1 ΣX2*X1 ΣX2*X2 … ΣX2*Xp ││ c │ │ ΣX2*Y │
! ： ： ： ：
! └ ΣXp*1 ΣXp*X1 ΣXp*X2 … ΣXp*Xp ┘└ d ┘ └ ΣXp*Y ┘
EXTERNAL SUB LINEST(Y(),X(,), t()) !係数を返す
DIM tW(P+1,N),W(N,P+1) !内積の計算に用いる
FOR i=1 TO N ![1,X1,X2, … ,Xp]
FOR k=1 TO P
LET W(i,k+1)=X(i,k)
NEXT k
LET W(i,1)=1
NEXT i
MAT tW=TRN(W)
DIM A(P+1,P+1),b(P+1) !連立方程式 A*t=b
MAT A=tW*W !左辺 A
MAT b=tW*Y !右辺 b
!!!MAT PRINT A; !debug
!!!MAT PRINT b;
DIM iA(P+1,P+1) !連立方程式を解く
MAT iA=INV(A)
MAT t=iA*b !求める係数
!!!MAT PRINT t; !debug
END SUB
EXTERNAL FUNCTION FORECAST(T(),Y(),X(,)) !予測値を求める
DIM C(P+1)
CALL LINEST(Y,X, C) !係数を求める
LET s=C(1)
FOR i=2 TO P+1 !近似式に代入する
LET s=s+T(i-1)*C(i)
NEXT i
LET FORECAST=s
END FUNCTION
EXTERNAL FUNCTION STEYX(Y(),X(,)) !予測値の標準誤差
DIM T(P)
LET s=0
FOR i=1 TO N
FOR k=1 TO P ![X1,X2,…,Xp]
LET T(k)=X(i,k)
NEXT k
LET s=s+(Y(i)-FORECAST(T,Y,X))^2
NEXT i
LET STEYX=SQR(s/(N-2)) !SQR( Σ(Y[k]-f[k])^2 / (N-2) )
END FUNCTION
EXTERNAL FUNCTION DEVSQ(Y()) !偏差平方和Σ[k=1,N]{(Y[k]-AY)^2}
LET s=0
LET s2=0
FOR i=1 TO N
LET s=s+Y(i)
LET s2=s2+Y(i)^2
NEXT i
LET DEVSQ=s2-s^2/N !ΣY*Y-(ΣY)^2/N
END FUNCTION
EXTERNAL FUNCTION RSQ(Y(),X(,)) !決定係数 R^2（寄与率）
DIM T(P)
LET s=0
FOR i=1 TO N
FOR k=1 TO P ![X1,X2,…,Xp]
LET T(k)=X(i,k)
NEXT k
LET s=s+(Y(i)-FORECAST(T,Y,X))^2
NEXT i
LET RSQ=1-s/DEVSQ(Y) !1 - Σ(Y[k]-f[k])^2 / Σ(Y[k]-AY)^2
END FUNCTION
!別解
EXTERNAL FUNCTION RSQ2(Y(),X(,)) !決定係数R^2
DIM W(N),T(P)
FOR i=1 TO N
FOR k=1 TO P ![X1,X2,…,Xp]
LET T(k)=X(i,k)
NEXT k
LET W(i)=FORECAST(T,Y,X) !予測値
NEXT i
LET RSQ2=DEVSQ(W)/DEVSQ(Y) !Σ(f[k]-Af)^2/Σ(Y[k]-AY)^2 推定値の分散を標本値の分散で割ったもの
END FUNCTION
|
|