コロナ関連の予測

 投稿者:大熊 正  投稿日:2020年 7月 9日(木)07時32分45秒
  今、エクセルでコロナ患者の増加曲線を作り予測等をしてます、式は、
 Y=A/(1+EXP(-M*(X-T)))  の
ロジスチック曲線で、データとYの差を2乗し、最小2乗を取り、エクセルのソルバーで、
Yの A,T,M を求めて、います。これを10進BASIC で同様に行いたいと思い色々やってますがうまくゆきません。すでにあるBASICのソフトの、組み合わせか投稿のソフトで同様のものがあるか、またはできるかを御教えをお願いします。
 

Re: コロナ関連の予測

 投稿者:nagram  投稿日:2020年 7月 9日(木)21時53分13秒
  > 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 "決定係数の2乗=";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 "決定係数の2乗=";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


!最小2乗法による
!
!残差平方和 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
 

コロナ関連の予測

 投稿者:大熊 正  投稿日:2020年 7月11日(土)14時35分41秒
  元記事 NO 4858  の御礼メールです。

早速の御返事ありがとう御座います。早速これからいろいろやってみます。

敬具
 

戻る