R の最適化関数を用いて SEM のパラメタを推定する方法
前川眞一
20250610,23,24,25,26,0708, 2025-07-08
この文書は R の Markdown で作成しています。灰色に囲まれた部分が R のコードで、 白枠の部分がその出力です。コードをコピーして R に貼り付けると同じ結果が 出力されるはずです。ただし、文書の最初から順に実行しないと、未定義のオブジェクト やパッケージが出てくる可能性があります。
コピーをするためには灰色部分の左上のボタンをクリックして下さい。
主なパッケージは、数値計算用の optimx です。 それに加えて、記号行列を取り扱うため lazy.symbolic パッケージと 数値微分を行うために lazy.mat パッケージもも必要となります。 これらのパッケージはここに置いてあります。
http://mayekawa.in.coocan.jp/Rpackages/Rpackages.html
# はじめに パッケージの宣言
options(width=200)
library(optimx)
library(lazy.symbolic)
library(lazy.mat)
modify_math_operators()（ここでいくつかのメッセージが出ますが省略してあります。）
1 RAM モデル
以下のパス図を例として解説する。
図 1 多重指標モデル
このパス図が表しているのは以下の方程式と同値である。 \[\begin{align}
Y_1 &= \lambda_1 \; \xi + E_1 \\
Y_2 &= \lambda_2 \; \xi + E_2 \\
Y_3 &= \lambda_3 \; \xi + E_3 \\
Y_4 &= \kappa_4 \; \eta + E_4 \\
Y_5 &= \kappa_5 \; \eta + E_5 \\
Y_6 &= \kappa_6 \; \eta + E_6 \\
\eta &= \gamma \; \xi + \zeta
\end{align}\]
この中で外生変数 exogenous variable は \(\{\xi, E_j,j=1,2,\ldots,6, \zeta\},\) であり 内生変数 endogenous variables は \(\{Y_j, j=1,2,\ldots,6, \eta\},\) である。
いま、それぞれの変数を縦に並べたベクトルを以下の様に定義する。 \[\begin{eqnarray} \boldsymbol y= \mbox{$\left(\begin{array}{c}Y_1 \\ Y_2 \\ \vdots \\ Y_6\end{array}\right)$},\hspace{2ex}\hspace{2ex} \boldsymbol e= \mbox{$\left(\begin{array}{c}E_1 \\ E_2 \\ \vdots \\ E_6\end{array}\right)$} \end{eqnarray}\] また、\(\boldsymbol v\) を観測変数と、誤差以外の分析の対象である潜在変数を縦に並べたベクトル、 \(\boldsymbol u\) を誤差を含む外性変数を縦に並べたベクトルとする。 \[\begin{align} \boldsymbol v= \mbox{$\left(\begin{array}{c} \boldsymbol y\\ \xi \\ \eta \end{array}\right)$},\hspace{2ex}\hspace{2ex}\boldsymbol u= \mbox{$\left(\begin{array}{c} \boldsymbol e\\\xi \\ \zeta \end{array}\right)$} \tag{vandu} \end{align}\] また、外生変数の分散共分散行列を \({{\cal D}}\left(\boldsymbol u\right)\) と置く。 このパス図の場合、\({{\cal D}}\left(\boldsymbol u\right)\) は \(\{\mbox{${\bf\it var}\left(E_j\right)$}, j=1,2,\ldots,6, \; \mbox{${\bf\it var}\left(\xi\right)$}, \mbox{${\bf\it var}\left(\zeta\right)$} \}\) を対角要素として持つ対角行列となる。外生変数間に相関を仮定する場合には \({{\cal D}}\left(\boldsymbol u\right)\) の対応する要素にその共分散に名前を付けたものを配置する。
そして、\(\boldsymbol v\) と \(\boldsymbol u\) の間の関係を以下の様に記す方法を RAM モデルと呼ぶ。 \[\begin{align} \boldsymbol v&= \boldsymbol B\boldsymbol v+ \boldsymbol u\tag{RAM1} \end{align}\] ただし、\(\boldsymbol B\) 行列の要素 \(b_{ij}\) は、 パス図に示される \(\boldsymbol v\) に含まれる変数 \(j\) から変数 \(i\) へのパス係数であり、先のパス図の場合には以下の様になる。 \[\begin{align} \boldsymbol B&= \begin{bmatrix} 0&0&0&0&0&0 & \lambda_1 & 0 \\ 0&0&0&0&0&0 & \lambda_2 & 0 \\ 0&0&0&0&0&0 & \lambda_3 & 0 \\ 0&0&0&0&0&0 & 0 & \kappa_4 \\ 0&0&0&0&0&0 & 0 & \kappa_5 \\ 0&0&0&0&0&0 & 0 & \kappa_6 \\ 0&0&0&0&0&0 & 0 & 0 \\ 0&0&0&0&0&0 & \gamma & 0 \end{bmatrix} \end{align}\] (RAM1) 式より、\(\boldsymbol v\) を移項することにより \[\begin{align} \boldsymbol v&= (\boldsymbol I-\boldsymbol B)^{-1} \boldsymbol u\tag{RAM2} \end{align}\] を得るため観測される の分散共分散行列構造は \[\begin{align} {{\cal D}}\left(\boldsymbol v\right) &= (\boldsymbol I-\boldsymbol B)^{-1} {{\cal D}}\left(\boldsymbol u\right) (\boldsymbol I-\boldsymbol B)^{-1}\mathstrut' \tag{RAM2} \end{align}\] として与えられる。
なお、上記の RAM モデルのパラメタは \(\boldsymbol B\) 行列の名前が付いた要素と \({{\cal D}}\left(\boldsymbol u\right)\) の名前が付いた要素である。
以下に、図 1 のパス図のパラメタを optimx や nlfb 等の R の最適化関数 を用いて推定する方法に関して説明する。
2 R の最適化関数を利用したモデルのパラメタの推定
2.1 データとして用いる観測変数間の分散共分散行列 \(\boldsymbol S\) を読み込む。
また、最尤解で必要となる \(\log|\boldsymbol S|\) を計算している。
# observed dispersion matrix from Toyoda Fig.11-4
# (Data from Hodge & Treiman (1968))
nvar=6
S=diag(nvar)
S[upper.tri(S,diag=TRUE)]=
c(
1.000,
0.304, 1.000,
0.305, 0.344, 1.000,
0.100, 0.156, 0.158, 1.000,
0.284, 0.192, 0.324, 0.360, 1.000,
0.176, 0.136, 0.226, 0.210, 0.265, 1.000
)
# sname=c("inc", "occ", "edu", "chu", "mem", "fri")
sname=paste("Y",1:nvar,sep="")
dimnames(S)=list(sname,sname)
S=S+t(S)
S=S-diag(diag(S)/2,ncol=nvar,nrow=nvar)
logdetS=log(det(S))
printm(S, logdetS, fmt="7.3")## S logdetS
## Y1 Y2 Y3 Y4 Y5 Y6
## Y1 1.000 0.304 0.305 0.100 0.284 0.176 -0.698
## Y2 0.304 1.000 0.344 0.156 0.192 0.136
## Y3 0.305 0.344 1.000 0.158 0.324 0.226
## Y4 0.100 0.156 0.158 1.000 0.360 0.210
## Y5 0.284 0.192 0.324 0.360 1.000 0.265
## Y6 0.176 0.136 0.226 0.210 0.265 1.000# save the copy
Sorg=S2.2 パラメタの数値をベクトル param として与えたとき、\({{\cal D}}\left(\boldsymbol v\right)\) を計算してくれる関数の準備。
R の最適化のための関数、例えば、optimx や nlfb は、パラメタの値を その要素として持つベクトルを第１引数として受け取り、第２引数として、 そのパラメタベクトルを第１引数として与えたときに最適化されるべき関数の値 （例えば最小２乗基準や尤度関数）を返す関数を受け取るように設計されている。
したがって、パラメタの値をその要素として持つベクトル、以下、パラメタベクトル、 を与えたときに、 モデルに基づいて観測される変数間の分散共分散行列を計算する関数が必要となるが、 それは、記号的に計算された \({{\cal D}}\left(\boldsymbol v\right)=(\boldsymbol I-\boldsymbol B)^{-1}{{\cal D}}\left(\boldsymbol u\right)(\boldsymbol I-\boldsymbol B)^{-1}\mathstrut'\) の要素として現れる記号的パラメタに、 対応するパラメタベクトルの数値を代入する関数となる。
以下に示す関数は、 パラメタベクトル param がその要素の記号名を names として持つ場合に、 \({{\cal D}}\left(\boldsymbol v\right)\) の記号的な式を与えたときに その中に現れるパラメタの要素に param で 指定されている数値を代入して、数値的な \({{\cal D}}\left(\boldsymbol v\right)\) を返す関数を生成する。
generate_Dv <- function( param, Dv, print=0 ){
# function factory to generate a function which creates numerical Dv matrix
# Shin-ichi Mayekawa
# 20250620,21
#
# param and Dv must be defined prior to call.
# top of the function
code0=c("f <- function( param, nvar ){"
," # Calculate Disp(v) from param"," ")
# definition of Dv
code=c(" # Get the elements of B matrix from param.")
code=c(code
, sapply( names(param)
, function(x) paste(" ", x,"=param[","'",x,"']", sep="")) )
code1=c(" "," # Define the Dv matrix.")
code1=c(code1," Dvn=diag(nvar)")
for( i in 1:length(Dv) ){
code1=c(code1,paste(" Dvn[",i,"]=",Dv[i], sep=""))
}
code2=c(" return(Dvn)","}")
# generate all
fullcode=c(code0,code,code1,code2)
# printm(code, .lj.=1)
eval(parse(text=fullcode))
if( print ){
cat("\nThe following function was generated.\n\n")
print(f)
if( print >= 2 ){
cat(" with\n")
# printm(code,code1,code2)
eval(parse(text=c(code,code1)))
printm(t(param))
}
}
return( f )
} # end of generate_Dv2.3 パス図を元に、\(\boldsymbol B\) 行列と \({{\cal D}}\left(\boldsymbol u\right)\) 行列を定義する。
まず、パス図が示唆する方程式を RAM モデルとして書くための準備を行う。 そのためには、\(\boldsymbol v\) と \(\boldsymbol u\) を定義し、係数行列 \(\boldsymbol B\) と \(\boldsymbol u\) の分散共分散行列 \({{\cal D}}\left(\boldsymbol u\right)\) を定めなければならない。
\(\boldsymbol v\) はパス図に現れる誤差以外の変数を表す。 変数の順序は、基本的には順不同であるが、ここでは、まず、 観測される変数 \(Y_j, j=1,2,\ldots,p\) を \(\boldsymbol v\) の要素とし、 その後に、誤差以外の潜在変数を追加する。
\(\boldsymbol u\) は \(\boldsymbol v\) と同じ長さのベクトルで、対応する \(\boldsymbol v\) の要素が外生変数の 場合にはその外生変数そのものが、内生変数の場合にはそれに加わる誤差が入る。
なお、\(\boldsymbol v\) と \(\boldsymbol u\) は以降で定義される行列の行や列の名前として、 また、係数行列が正しくパス図を表現しているかを確かめるために利用される。
options(width=200)
# variables
nvar=6
# variables of interest
v=c(paste("Y",1:nvar,sep=""),"xi", "eta" )
# Exogenous
u=c(paste("E",1:nvar,sep=""),"xi", "zeta")
nv=length(v)
locy=1:nvar
printm(nvar,nv,locy,v,u)## nvar nv locy v u
## 6 8 1 1 1 Y1 1 E1
## 2 2 2 Y2 2 E2
## 3 3 3 Y3 3 E3
## 4 4 4 Y4 4 E4
## 5 5 5 Y5 5 E5
## 6 6 6 Y6 6 E6
## 7 xi 7 xi
## 8 eta 8 zeta次に、\({{\cal D}}\left(\boldsymbol u\right)\) に含まれる分散共分散行列の要素のうち、 自由パラメタとなる部分に適当に名前を付ける。 分散の名前を v から始まるように指定すると、最尤法ではその対数がパラメタとなる。
また、分散共分散行列の要素のうち既知の部分は該当する数値を入れておく。 以下の例では、大部分が 0 である。
また、外生変数である \(\boldsymbol\xi\) の分散を 1 に固定している。
options(width=200)
# dispersion of u
Du=Diag(c(paste("ve",1:nvar, sep=""),"vxi","vzeta"))
dimnames(Du)=list(u,u)
# fix the variance of xi to 1.
Du["xi","xi"]=1
printm(v,u,Du)## v u Du
## E1 E2 E3 E4 E5 E6 xi zeta
## 1 Y1 1 E1 E1 ve1 0 0 0 0 0 0 0
## 2 Y2 2 E2 E2 0 ve2 0 0 0 0 0 0
## 3 Y3 3 E3 E3 0 0 ve3 0 0 0 0 0
## 4 Y4 4 E4 E4 0 0 0 ve4 0 0 0 0
## 5 Y5 5 E5 E5 0 0 0 0 ve5 0 0 0
## 6 Y6 6 E6 E6 0 0 0 0 0 ve6 0 0
## 7 xi 7 xi xi 0 0 0 0 0 0 1 0
## 8 eta 8 zeta zeta 0 0 0 0 0 0 0 vzeta最後に \(\boldsymbol B\) に含まれる係数のうち自由パラメタとなる部分に適当に名前を付ける。 名前の付け方は自由であるが v で始まるものは避ける。
\(\boldsymbol B\) 行列の \(i,j\) 要素には \(v_j\) から \(v_i\) へのパスが存在するときに 何らかの名称もしくは数値が付与される。 ここでは \(v_j\) から \(v_i\) へのパスの 係数を \(ij\) とを付加した形で命名している。
以下の例では \(\boldsymbol\eta\) から \(Y_4\) へのパス係数を 1 に固定している。
options(width=200)
#
# path matrix
# B[i,j] is the name of the coefficient of the path from j to i.
B=matrix(0,nv,nv)
dimnames(B)=list(v,v)
B[8,7]="G"
B[1:3,7]=paste(paste("L",1:3,sep=""),"",sep="")
B[4:6,8]=paste(paste("K",1:3,sep=""),"",sep="")
# fix kappa4 to 1
B["Y4","eta"]=1
printm(B)## B
## Y1 Y2 Y3 Y4 Y5 Y6 xi eta
## Y1 0 0 0 0 0 0 L1 0
## Y2 0 0 0 0 0 0 L2 0
## Y3 0 0 0 0 0 0 L3 0
## Y4 0 0 0 0 0 0 0 1
## Y5 0 0 0 0 0 0 0 K2
## Y6 0 0 0 0 0 0 0 K3
## xi 0 0 0 0 0 0 0 0
## eta 0 0 0 0 0 0 G 0パス図に関する入力は以上で終わりである。
2.4 \(\boldsymbol B\) 行列と \({{\cal D}}\left(\boldsymbol u\right)\) 行列から記号的な \({{\cal D}}\left(\boldsymbol v\right)\) を算出する。
\(\boldsymbol B\) 行列が定義されると、 \((\boldsymbol I-\boldsymbol B)^{-1}\) の逆行列を記号的に求めて最終的に \({{\cal D}}\left(\boldsymbol v\right)\) を 記号的に定義する。
また、\(\boldsymbol B\boldsymbol v+\boldsymbol u\) や \((\boldsymbol I-\boldsymbol B)^{-1}\boldsymbol u\) が出力され、パス図との対応を確かめる ことが出来る。
# sort S, u, v, B so that B becomes lower triangular
options(width=120)
# calculate inv(I-B)
invImB=matSweep(diag(nv)-B) |> Expand(all=1)
printm(B,invImB,Du)## B invImB Du
## Y1 Y2 Y3 Y4 Y5 Y6 xi eta Y1 Y2 Y3 Y4 Y5 Y6 xi eta E1 E2 E3 E4 E5 E6 xi zeta
## Y1 0 0 0 0 0 0 L1 0 Y1 1 0 0 0 0 0 L1 0 E1 ve1 0 0 0 0 0 0 0
## Y2 0 0 0 0 0 0 L2 0 Y2 0 1 0 0 0 0 L2 0 E2 0 ve2 0 0 0 0 0 0
## Y3 0 0 0 0 0 0 L3 0 Y3 0 0 1 0 0 0 L3 0 E3 0 0 ve3 0 0 0 0 0
## Y4 0 0 0 0 0 0 0 1 Y4 0 0 0 1 0 0 G 1 E4 0 0 0 ve4 0 0 0 0
## Y5 0 0 0 0 0 0 0 K2 Y5 0 0 0 0 1 0 K2*G K2 E5 0 0 0 0 ve5 0 0 0
## Y6 0 0 0 0 0 0 0 K3 Y6 0 0 0 0 0 1 K3*G K3 E6 0 0 0 0 0 ve6 0 0
## xi 0 0 0 0 0 0 0 0 xi 0 0 0 0 0 0 1 0 xi 0 0 0 0 0 0 1 0
## eta 0 0 0 0 0 0 G 0 eta 0 0 0 0 0 0 G 1 zeta 0 0 0 0 0 0 0 vzetaprintm(B%*%v+u, invImB%*%u)## B %*% v + u invImB %*% u
## [,1] [,1]
## Y1 L1*xi+E1 Y1 E1+L1*xi
## Y2 L2*xi+E2 Y2 E2+L2*xi
## Y3 L3*xi+E3 Y3 E3+L3*xi
## Y4 eta+E4 Y4 E4+G*xi+zeta
## Y5 K2*eta+E5 Y5 E5+K2*G*xi+K2*zeta
## Y6 K3*eta+E6 Y6 E6+K3*G*xi+K3*zeta
## xi xi xi xi
## eta G*xi+zeta eta G*xi+zeta# calculate D(v)=inv(I-B) D(u) inv(I-B)'
Dv=invImB%*%Du%*%t(invImB)
printm(Dv, .lj.=1)## Dv
## Y1 Y2 Y3 Y4
## Y1 ve1+L1*L1 L1*L2 L1*L3 L1*G
## Y2 L2*L1 ve2+L2*L2 L2*L3 L2*G
## Y3 L3*L1 L3*L2 ve3+L3*L3 L3*G
## Y4 G*L1 G*L2 G*L3 ve4+G*G+vzeta
## Y5 K2*G*L1 K2*G*L2 K2*G*L3 K2*G*G+K2*vzeta
## Y6 K3*G*L1 K3*G*L2 K3*G*L3 K3*G*G+K3*vzeta
## xi L1 L2 L3 G
## eta G*L1 G*L2 G*L3 G*G+vzeta
## Y5 Y6 xi eta
## Y1 L1*K2*G L1*K3*G L1 L1*G
## Y2 L2*K2*G L2*K3*G L2 L2*G
## Y3 L3*K2*G L3*K3*G L3 L3*G
## Y4 G*K2*G+vzeta*K2 G*K3*G+vzeta*K3 G G*G+vzeta
## Y5 ve5+K2*G*K2*G+K2*vzeta*K2 K2*G*K3*G+K2*vzeta*K3 K2*G K2*G*G+K2*vzeta
## Y6 K3*G*K2*G+K3*vzeta*K2 ve6+K3*G*K3*G+K3*vzeta*K3 K3*G K3*G*G+K3*vzeta
## xi K2*G K3*G 1 G
## eta G*K2*G+vzeta*K2 G*K3*G+vzeta*K3 G G*G+vzetaこれで、必要となる \({{\cal D}}\left(\boldsymbol v\right)\) が定義できた。
2.5 パラメタベクトルの保管場所を定義する。
次に初期値となる数値的パラメタベクトル praram を定めるが、 param の各要素がどのパラメタの値であるかは、要素の名称 names(param) として 与えられていなければならない。
ここでは param の要素の名前のみを生成し、数値は仮の値を入れておく。 実際の初期値は、関数の最適化の直前に代入する。
options(width=200)
# Parameter vector
bname=unique( B[!sapply(B, is.number, pat=2)] )
dname=unique( Du[!sapply(Du, is.number, pat=2)] )
paramname=c(bname,dname)
nparam=length(paramname)
param=1:nparam
names(param)=paramname
printm(param)## param
## L1 1
## L2 2
## L3 3
## G 4
## K2 5
## K3 6
## ve1 7
## ve2 8
## ve3 9
## ve4 10
## ve5 11
## ve6 12
## vzeta 132.6 R の最適化関数を用いるための準備
まず、先ほど作成した generate_Dv を用いてパラメタベクトルを与えたときに \({{\cal D}}\left(\boldsymbol v\right)\) を数値的に計算してくれる関数 get_Dv を作成する。 ただし、必要となるのは \({{\cal D}}\left(\boldsymbol y\right)\) の部分のみであるため、引数 Dv として Dv[locy,locy] を渡している。
options(width=200)
get_Dy <- generate_Dv( param, Dv[locy,locy], print=2 )##
## The following function was generated.
##
## function (param, nvar)
## {
## L1 = param["L1"]
## L2 = param["L2"]
## L3 = param["L3"]
## G = param["G"]
## K2 = param["K2"]
## K3 = param["K3"]
## ve1 = param["ve1"]
## ve2 = param["ve2"]
## ve3 = param["ve3"]
## ve4 = param["ve4"]
## ve5 = param["ve5"]
## ve6 = param["ve6"]
## vzeta = param["vzeta"]
## Dvn = diag(nvar)
## Dvn[1] = ve1 + L1 * L1
## Dvn[2] = L2 * L1
## Dvn[3] = L3 * L1
## Dvn[4] = G * L1
## Dvn[5] = K2 * G * L1
## Dvn[6] = K3 * G * L1
## Dvn[7] = L1 * L2
## Dvn[8] = ve2 + L2 * L2
## Dvn[9] = L3 * L2
## Dvn[10] = G * L2
## Dvn[11] = K2 * G * L2
## Dvn[12] = K3 * G * L2
## Dvn[13] = L1 * L3
## Dvn[14] = L2 * L3
## Dvn[15] = ve3 + L3 * L3
## Dvn[16] = G * L3
## Dvn[17] = K2 * G * L3
## Dvn[18] = K3 * G * L3
## Dvn[19] = L1 * G
## Dvn[20] = L2 * G
## Dvn[21] = L3 * G
## Dvn[22] = ve4 + G * G + vzeta
## Dvn[23] = K2 * G * G + K2 * vzeta
## Dvn[24] = K3 * G * G + K3 * vzeta
## Dvn[25] = L1 * K2 * G
## Dvn[26] = L2 * K2 * G
## Dvn[27] = L3 * K2 * G
## Dvn[28] = G * K2 * G + vzeta * K2
## Dvn[29] = ve5 + K2 * G * K2 * G + K2 * vzeta * K2
## Dvn[30] = K3 * G * K2 * G + K3 * vzeta * K2
## Dvn[31] = L1 * K3 * G
## Dvn[32] = L2 * K3 * G
## Dvn[33] = L3 * K3 * G
## Dvn[34] = G * K3 * G + vzeta * K3
## Dvn[35] = K2 * G * K3 * G + K2 * vzeta * K3
## Dvn[36] = ve6 + K3 * G * K3 * G + K3 * vzeta * K3
## return(Dvn)
## }
## <environment: 0x000001a267524b58>
## with
## t(param)
## L1 L2 L3 G K2 K3 ve1 ve2 ve3 ve4 ve5 ve6 vzeta
## [1,] 1 2 3 4 5 6 7 8 9 10 11 12 13そして、それを用いて、最適化基準を計算する関数を定義する。 最適化基準としては、最小２乗解の場合には残差の平方和を分散共分散行列の 下三角部分の要素の数で除したものの平方根を、最尤解の場合には負の対数尤度比 を用いる。また、最尤解の場合、\({{\cal D}}\left(\boldsymbol u\right)\) の分散成分の対数をパラメタと して用いるため、\({{\cal D}}\left(\boldsymbol v\right)\) の計算に際してそれを元に戻している。
# rmse
rmse <- function( param ){
# rmse of LS PA ssq(vec(S-Dv)
# Shin-ichi Mayekawa
# 20220617
#
# Note that S and nvar are from the parent environment
# where residSDv2 is defined.
# get Dv from param
Dv=get_Dy( param, nvar )
# rmse
resid=vech(S-Dv)
return( sqrt(sum(resid^2)/length(resid)) )
} # end of rmse
# mllh
mllh <- function( param ){
# minus twice the log likelihood (mllh) of PA
# Shin-ichi Mayekawa
# 20220618
# log variance: 20250622
#
# Note that S and logdetS and nvar are from
# the parent environment where mllh is defined.
# log variance
locv=which( substr(names(param),1,1) == "v" )
p=param
p[locv]=exp(p[locv])
# get Dv from param
Dv=get_Dy( p, nvar )
mllh=tr(S%*%solve(Dv))+log(det(Dv)) - logdetS - nvar
return( mllh )
} # end of mllhpara
d1 <- function( param, func ){
# numerical first derivative of function func.
grad=JacobianMat( param, func )
return( grad )
} # end of d1
d2 <- function( param, func ){
# numerical second derivative of func
res=JacobianMat( param, d1, func=func )
# cat(".d2an.")
return( res )
} # end of d22.7 R の最適化関数を利用したパラメタの推定
これで準備が完了したので、後は、R の何らかの最適化関数を用いて パラメタを推定すれば良い。
パラメタの初期値としては何らかの乱数を用い手も良いが、 ここではパラメタの値を全て 1 としている。 ただ、そおれらの値を param に代入する際に
param[] <- 数値ベクトル
とすると、先に定義した要素の名前が保持される。
パラメタの初期値は大切である。本来ならばもっと良い初期値を用いるべきである。
最適化のための関数としては optimx を用い、方法として nlminb と BFGS を指定した。
最小２乗解：
options(width=200)
# initial value of param: Names must be kept.
set.seed(1701)
param[]=rnorm(length(param),2,0.1)
param[]=rep(1,length(param))
rmse0=rmse(param)
# estimate parameters by optimx
method=c("Nelder-Mead","nlminb","BFGS")
resoptLS=optimx( param, rmse, method=method, control=list(trace=0) )
rmse1=resoptLS$value
names(rmse1)=method
paramLS=t(resoptLS[,1:nparam, drop=0])
printm(param, paramLS, fmt="12.5")## param paramLS
## Nelder-Mead nlminb BFGS
## L1 1.00000 L1 0.35497 0.53002 0.53002
## L2 1.00000 L2 0.31829 0.51464 0.51464
## L3 1.00000 L3 0.51157 0.64477 0.64477
## G 1.00000 G 0.40707 0.29437 0.29436
## K2 1.00000 K2 2.02102 1.59626 1.59635
## K3 1.00000 K3 1.52564 0.95208 0.95206
## ve1 1.00000 ve1 0.98408 0.71908 0.71907
## ve2 1.00000 ve2 1.23464 0.73514 0.73514
## ve3 1.00000 ve3 0.89196 0.58427 0.58427
## ve4 1.00000 ve4 1.32971 0.79575 0.79575
## ve5 1.00000 ve5 0.98389 0.47956 0.47952
## ve6 1.00000 ve6 1.18198 0.81486 0.81487
## vzeta 1.00000 vzeta -0.12520 0.11759 0.11759printm(rmse0,t(rmse1), fmt="12.5")## rmse0 t(rmse1)
## Nelder-Mead nlminb BFGS
## 1.21851 [1,] 0.17203 0.02855 0.02855# gradient
gradLS=apply(paramLS,2, d1, rmse)
rownames(gradLS)=paramname
printm(gradLS, fmt="16.7")## gradLS
## Nelder-Mead nlminb BFGS
## L1 0.0018149 -0.0000002 -0.0000096
## L2 0.0383018 0.0000000 -0.0000068
## L3 0.0589139 0.0000000 -0.0000083
## G 0.1764601 0.0000002 -0.0000165
## K2 0.0096616 0.0000000 0.0000086
## K3 0.0138896 -0.0000000 -0.0000104
## ve1 0.0304738 0.0000004 -0.0000029
## ve2 0.0929919 -0.0000001 -0.0000015
## ve3 0.0425339 0.0000002 -0.0000021
## ve4 0.1024779 0.0000001 -0.0000077
## ve5 0.0413325 0.0000000 0.0000025
## ve6 0.0764692 -0.0000002 -0.0000024
## vzeta 0.1115107 0.0000001 0.0000013初期値で計算された rmse の値は 1.218509 であり、最終的に得られた最小２乗解では その値が 0.1720306, 0.0285524, 0.0285524 にまで減少している。 数値的に計算したグレイディエントからみて、 Nelder-Mead 以外は収束しているように見える。
パラメタの２次微分の計算
モデルが識別されない場合には最適化される関数のパラメタベクトルに関する ２次微分行列が正則でなくなる。
HessLS=d2(paramLS[,2], rmse)
dimnames(HessLS)=list(paramname,paramname)
printm(HessLS, fmt="6.3")## HessLS
## L1 L2 L3 G K2 K3 ve1 ve2 ve3 ve4 ve5 ve6 vzeta
## L1 3.653 0.403 0.631 1.116 0.105 0.059 1.768 0.000 0.000 0.000 0.000 -0.000 -0.000
## L2 0.403 3.573 0.533 1.264 0.143 0.075 0.000 1.717 0.000 -0.000 0.000 -0.000 0.000
## L3 0.631 0.533 4.327 1.335 0.138 0.066 0.000 0.000 2.151 -0.000 0.000 0.000 0.000
## G 1.116 1.264 1.335 15.282 2.994 1.463 0.000 0.000 0.000 0.982 2.502 0.890 13.824
## K2 0.105 0.143 0.138 2.994 0.981 0.121 0.000 0.000 0.000 0.000 1.088 -0.000 3.823
## K3 0.059 0.075 0.066 1.463 0.121 0.638 -0.000 0.000 0.000 0.000 -0.000 0.649 1.834
## ve1 1.768 0.000 0.000 0.000 0.000 -0.000 1.668 0.000 0.000 0.000 0.000 0.000 0.000
## ve2 0.000 1.717 0.000 0.000 0.000 0.000 0.000 1.668 0.000 0.000 0.000 -0.000 0.000
## ve3 0.000 0.000 2.151 0.000 0.000 0.000 0.000 0.000 1.668 0.000 0.000 0.000 0.000
## ve4 0.000 -0.000 -0.000 0.982 0.000 0.000 0.000 0.000 0.000 1.668 0.000 0.000 1.668
## ve5 0.000 0.000 0.000 2.502 1.088 -0.000 0.000 0.000 0.000 0.000 1.668 0.000 4.250
## ve6 -0.000 -0.000 0.000 0.890 -0.000 0.649 0.000 -0.000 0.000 0.000 0.000 1.668 1.512
## vzeta -0.000 0.000 0.000 13.824 3.823 1.834 0.000 0.000 0.000 1.668 4.250 1.512 23.480printm(nparam,qr(HessLS)$rank)## nparam qr(HessLS)$rank
## 13 13printm(eigen(HessLS)$values, fmt="10.5")## eigen(HessLS)$values
## 1 35.65133
## 2 7.10228
## 3 4.75086
## 4 4.35740
## 5 3.79691
## 6 2.04374
## 7 1.73728
## 8 0.73683
## 9 0.67180
## 10 0.53571
## 11 0.33479
## 12 0.16810
## 13 0.05281この場合、Hessian のランクも落ちていない。
最尤解：
初期値は最小２乗解を用いる。最適化のための関数としては optimx を用い、方法として Nelder-Mead, nlminb ならびに BFGS を指定した。 また、\({{\cal D}}\left(\boldsymbol u\right)\) のうち分散成分はその対数をパラメタとし、 負にならないようにしている。
options(width=200)
# initial from LS
param=unlist(resoptLS["nlminb",1:nparam])
# log(variance) as the paramrter
locv=which( substr(names(param),1,1) == "v" )
param[locv]=log(pmax(0.01,param[locv]))
mllh0=mllh(param)
# estimate parameters by optimx
options(warn=-1)
method=c("Nelder-Mead","nlminb","BFGS")
resoptML=optimx( param, mllh, method=method
, control=list(trace=0, maxit=500
, rel.tol=1e-15, step.min=1/16, x.tol=1e-15) )
options(warn=1)
resoptML[abs(resoptML) > 1e9]=1e9
paramML=t(resoptML[,1:nparam, drop=0])
paramML0=paramML
# back to original
locv=which( substr(names(param),1,1) == "v" )
paramML[locv,]=exp(paramML[locv,])
printm(param,paramML, fmt="12.5")## param paramML
## Nelder-Mead nlminb BFGS
## L1 0.53002 L1 0.53332 0.52626 0.52626
## L2 0.51464 L2 0.52093 0.52515 0.52515
## L3 0.64477 L3 0.64111 0.63513 0.63513
## G 0.29437 G 0.30153 0.31033 0.31033
## K2 1.59626 K2 1.57842 1.52412 1.52412
## K3 0.95208 K3 0.90969 0.85116 0.85116
## ve1 -0.32979 ve1 0.73069 0.72306 0.72306
## ve2 -0.30769 ve2 0.72826 0.72421 0.72422
## ve3 -0.53739 ve3 0.58909 0.59660 0.59661
## ve4 -0.22847 ve4 0.79383 0.77437 0.77437
## ve5 -0.73488 ve5 0.46750 0.47589 0.47588
## ve6 -0.20474 ve6 0.82594 0.83654 0.83654
## vzeta -2.14053 vzeta 0.12322 0.12932 0.12932printm(mllh0,resoptML$value, fmt="12.5" )## mllh0 resoptML$value
## 0.03424 1 0.03362
## 2 0.03295
## 3 0.03295# gradient
gradML=apply(paramML0,2, d1, mllh)
rownames(gradML)=paramname
printm(gradML, fmt="16.7")## gradML
## Nelder-Mead nlminb BFGS
## L1 0.0168867 -0.0000000 0.0000006
## L2 -0.0052748 -0.0000000 0.0000008
## L3 0.0011784 -0.0000000 0.0000001
## G 0.0035928 0.0000000 0.0000002
## K2 0.0008277 -0.0000000 0.0000001
## K3 0.0110041 0.0000000 0.0000003
## ve1 0.0101167 0.0000001 0.0000009
## ve2 0.0007346 0.0000000 0.0000005
## ve3 -0.0033929 0.0000000 0.0000007
## ve4 0.0141633 0.0000000 -0.0000001
## ve5 -0.0035172 -0.0000000 -0.0000003
## ve6 -0.0070184 -0.0000000 -0.0000000
## vzeta 0.0025392 -0.0000000 0.0000001数値的に計算したグレイディエントから見て Nelder-Mead 以外は収束している。
パラメタの２次微分の計算
HessML=d2(paramML0[,2], mllh)
dimnames(HessML)=list(paramname,paramname)
printm(HessML, fmt="6.3")## HessML
## L1 L2 L3 G K2 K3 ve1 ve2 ve3 ve4 ve5 ve6 vzeta
## L1 1.668 -0.147 -0.110 -0.205 -0.058 -0.016 0.449 -0.050 -0.157 -0.016 -0.012 -0.005 -0.028
## L2 -0.147 1.667 -0.177 -0.071 -0.044 -0.014 -0.052 0.447 -0.115 0.001 -0.045 -0.006 -0.045
## L3 -0.110 -0.177 1.925 -0.301 -0.076 -0.028 -0.131 -0.095 0.578 -0.012 -0.028 -0.004 -0.039
## G -0.205 -0.071 -0.301 4.466 0.935 0.278 -0.094 -0.163 -0.163 0.045 0.265 0.063 0.341
## K2 -0.058 -0.044 -0.076 0.935 0.512 -0.023 -0.004 -0.011 -0.014 -0.041 0.308 -0.043 0.208
## K3 -0.016 -0.014 -0.028 0.278 -0.023 0.346 -0.002 -0.002 0.000 -0.012 -0.074 0.147 0.075
## ve1 0.449 -0.052 -0.131 -0.094 -0.004 -0.002 0.745 0.007 0.048 0.007 -0.002 -0.001 0.004
## ve2 -0.050 0.447 -0.095 -0.163 -0.011 -0.002 0.007 0.746 0.025 0.000 0.014 0.000 0.017
## ve3 -0.157 -0.115 0.578 -0.163 -0.014 0.000 0.048 0.025 0.576 0.008 0.006 -0.004 0.011
## ve4 -0.016 0.001 -0.012 0.045 -0.041 -0.012 0.007 0.000 0.008 0.812 0.028 0.003 0.046
## ve5 -0.012 -0.045 -0.028 0.265 0.308 -0.074 -0.002 0.014 0.006 0.028 0.389 0.038 0.130
## ve6 -0.005 -0.006 -0.004 0.063 -0.043 0.147 -0.001 0.000 -0.004 0.003 0.038 0.872 0.016
## vzeta -0.028 -0.045 -0.039 0.341 0.208 0.075 0.004 0.017 0.011 0.046 0.130 0.016 0.172printm(nparam,qr(HessML)$rank)## nparam qr(HessML)$rank
## 13 13printm(eigen(HessML)$values, fmt="10.5")## eigen(HessML)$values
## 1 4.81865
## 2 2.26652
## 3 2.00834
## 4 1.50490
## 5 0.91682
## 6 0.81953
## 7 0.71179
## 8 0.60637
## 9 0.53218
## 10 0.30364
## 11 0.27867
## 12 0.09312
## 13 0.03599この場合、Hessian のランクも落ちていない。
2.8 識別されないモデル。
以下に \(\boldsymbol\eta\) から \(Y_4\) へのパス係数を自由なパラメタとした例を示す。 変更する部分は \(\boldsymbol B\) 行列の定義部分とそれに続くパラメタの定義部分、 ならびに、パラメタを与えたときに \({{\cal D}}\left(\boldsymbol v\right)\) を計算する関数の生成部分である。
モデルが識別されない場合、 最適化される関数のパラメタベクトルによる２次微分の行列、 すなわち Hessian 行列が正則ではなくなる。
options(width=200)
# same Du will be used.
# path matrix
# B[i,j] is the name of the coefficient of the path from j to i.
B=matrix(0,nv,nv)
dimnames(B)=list(v,v)
B[8,7]="G"
B[1:3,7]=paste(paste("L",1:3,sep=""),"",sep="")
B[4:6,8]=paste(paste("K",1:3,sep=""),"",sep="")
# do NOT fix kappa4 to 1
#### B["Y4","eta"]=1
printm(B,Du)## B Du
## Y1 Y2 Y3 Y4 Y5 Y6 xi eta E1 E2 E3 E4 E5 E6 xi zeta
## Y1 0 0 0 0 0 0 L1 0 E1 ve1 0 0 0 0 0 0 0
## Y2 0 0 0 0 0 0 L2 0 E2 0 ve2 0 0 0 0 0 0
## Y3 0 0 0 0 0 0 L3 0 E3 0 0 ve3 0 0 0 0 0
## Y4 0 0 0 0 0 0 0 K1 E4 0 0 0 ve4 0 0 0 0
## Y5 0 0 0 0 0 0 0 K2 E5 0 0 0 0 ve5 0 0 0
## Y6 0 0 0 0 0 0 0 K3 E6 0 0 0 0 0 ve6 0 0
## xi 0 0 0 0 0 0 0 0 xi 0 0 0 0 0 0 1 0
## eta 0 0 0 0 0 0 G 0 zeta 0 0 0 0 0 0 0 vzeta# calculate inv(I-B)
invImB=matSweep(diag(nv)-B) |> Expand(all=1)
printm(B,invImB,Du)## B invImB Du
## Y1 Y2 Y3 Y4 Y5 Y6 xi eta Y1 Y2 Y3 Y4 Y5 Y6 xi eta E1 E2 E3 E4 E5 E6 xi zeta
## Y1 0 0 0 0 0 0 L1 0 Y1 1 0 0 0 0 0 L1 0 E1 ve1 0 0 0 0 0 0 0
## Y2 0 0 0 0 0 0 L2 0 Y2 0 1 0 0 0 0 L2 0 E2 0 ve2 0 0 0 0 0 0
## Y3 0 0 0 0 0 0 L3 0 Y3 0 0 1 0 0 0 L3 0 E3 0 0 ve3 0 0 0 0 0
## Y4 0 0 0 0 0 0 0 K1 Y4 0 0 0 1 0 0 K1*G K1 E4 0 0 0 ve4 0 0 0 0
## Y5 0 0 0 0 0 0 0 K2 Y5 0 0 0 0 1 0 K2*G K2 E5 0 0 0 0 ve5 0 0 0
## Y6 0 0 0 0 0 0 0 K3 Y6 0 0 0 0 0 1 K3*G K3 E6 0 0 0 0 0 ve6 0 0
## xi 0 0 0 0 0 0 0 0 xi 0 0 0 0 0 0 1 0 xi 0 0 0 0 0 0 1 0
## eta 0 0 0 0 0 0 G 0 eta 0 0 0 0 0 0 G 1 zeta 0 0 0 0 0 0 0 vzetaprintm(B%*%v+u, invImB%*%u)## B %*% v + u invImB %*% u
## [,1] [,1]
## Y1 L1*xi+E1 Y1 E1+L1*xi
## Y2 L2*xi+E2 Y2 E2+L2*xi
## Y3 L3*xi+E3 Y3 E3+L3*xi
## Y4 K1*eta+E4 Y4 E4+K1*G*xi+K1*zeta
## Y5 K2*eta+E5 Y5 E5+K2*G*xi+K2*zeta
## Y6 K3*eta+E6 Y6 E6+K3*G*xi+K3*zeta
## xi xi xi xi
## eta G*xi+zeta eta G*xi+zeta# calculate D(v)=inv(I-B) D(u) inv(I-B)'
Dv=invImB%*%Du%*%t(invImB)
printm(Dv, .lj.=1)## Dv
## Y1 Y2 Y3 Y4 Y5 Y6 xi
## Y1 ve1+L1*L1 L1*L2 L1*L3 L1*K1*G L1*K2*G L1*K3*G L1
## Y2 L2*L1 ve2+L2*L2 L2*L3 L2*K1*G L2*K2*G L2*K3*G L2
## Y3 L3*L1 L3*L2 ve3+L3*L3 L3*K1*G L3*K2*G L3*K3*G L3
## Y4 K1*G*L1 K1*G*L2 K1*G*L3 ve4+K1*G*K1*G+K1*vzeta*K1 K1*G*K2*G+K1*vzeta*K2 K1*G*K3*G+K1*vzeta*K3 K1*G
## Y5 K2*G*L1 K2*G*L2 K2*G*L3 K2*G*K1*G+K2*vzeta*K1 ve5+K2*G*K2*G+K2*vzeta*K2 K2*G*K3*G+K2*vzeta*K3 K2*G
## Y6 K3*G*L1 K3*G*L2 K3*G*L3 K3*G*K1*G+K3*vzeta*K1 K3*G*K2*G+K3*vzeta*K2 ve6+K3*G*K3*G+K3*vzeta*K3 K3*G
## xi L1 L2 L3 K1*G K2*G K3*G 1
## eta G*L1 G*L2 G*L3 G*K1*G+vzeta*K1 G*K2*G+vzeta*K2 G*K3*G+vzeta*K3 G
## eta
## Y1 L1*G
## Y2 L2*G
## Y3 L3*G
## Y4 K1*G*G+K1*vzeta
## Y5 K2*G*G+K2*vzeta
## Y6 K3*G*G+K3*vzeta
## xi G
## eta G*G+vzeta# Parameter vector
bname=unique( B[!sapply(B, is.number, pat=2)] )
dname=unique( Du[!sapply(Du, is.number, pat=2)] )
paramname=c(bname,dname)
nparam=length(paramname)
param=1:nparam
names(param)=paramname
printm(param)## param
## L1 1
## L2 2
## L3 3
## G 4
## K1 5
## K2 6
## K3 7
## ve1 8
## ve2 9
## ve3 10
## ve4 11
## ve5 12
## ve6 13
## vzeta 14# new dispersion function
get_Dy <- generate_Dv( param, Dv[locy,locy], print=2 )##
## The following function was generated.
##
## function (param, nvar)
## {
## L1 = param["L1"]
## L2 = param["L2"]
## L3 = param["L3"]
## G = param["G"]
## K1 = param["K1"]
## K2 = param["K2"]
## K3 = param["K3"]
## ve1 = param["ve1"]
## ve2 = param["ve2"]
## ve3 = param["ve3"]
## ve4 = param["ve4"]
## ve5 = param["ve5"]
## ve6 = param["ve6"]
## vzeta = param["vzeta"]
## Dvn = diag(nvar)
## Dvn[1] = ve1 + L1 * L1
## Dvn[2] = L2 * L1
## Dvn[3] = L3 * L1
## Dvn[4] = K1 * G * L1
## Dvn[5] = K2 * G * L1
## Dvn[6] = K3 * G * L1
## Dvn[7] = L1 * L2
## Dvn[8] = ve2 + L2 * L2
## Dvn[9] = L3 * L2
## Dvn[10] = K1 * G * L2
## Dvn[11] = K2 * G * L2
## Dvn[12] = K3 * G * L2
## Dvn[13] = L1 * L3
## Dvn[14] = L2 * L3
## Dvn[15] = ve3 + L3 * L3
## Dvn[16] = K1 * G * L3
## Dvn[17] = K2 * G * L3
## Dvn[18] = K3 * G * L3
## Dvn[19] = L1 * K1 * G
## Dvn[20] = L2 * K1 * G
## Dvn[21] = L3 * K1 * G
## Dvn[22] = ve4 + K1 * G * K1 * G + K1 * vzeta * K1
## Dvn[23] = K2 * G * K1 * G + K2 * vzeta * K1
## Dvn[24] = K3 * G * K1 * G + K3 * vzeta * K1
## Dvn[25] = L1 * K2 * G
## Dvn[26] = L2 * K2 * G
## Dvn[27] = L3 * K2 * G
## Dvn[28] = K1 * G * K2 * G + K1 * vzeta * K2
## Dvn[29] = ve5 + K2 * G * K2 * G + K2 * vzeta * K2
## Dvn[30] = K3 * G * K2 * G + K3 * vzeta * K2
## Dvn[31] = L1 * K3 * G
## Dvn[32] = L2 * K3 * G
## Dvn[33] = L3 * K3 * G
## Dvn[34] = K1 * G * K3 * G + K1 * vzeta * K3
## Dvn[35] = K2 * G * K3 * G + K2 * vzeta * K3
## Dvn[36] = ve6 + K3 * G * K3 * G + K3 * vzeta * K3
## return(Dvn)
## }
## <environment: 0x000001a26c143ed8>
## with
## t(param)
## L1 L2 L3 G K1 K2 K3 ve1 ve2 ve3 ve4 ve5 ve6 vzeta
## [1,] 1 2 3 4 5 6 7 8 9 10 11 12 13 14# initial value of param: Names must be kept.
set.seed(1701)
param[]=rnorm(length(param),2,0.1)
param[]=rep(1,length(param))
rmse0=rmse(param)
# estimate parameters by optimx
method=c("Nelder-Mead","nlminb","BFGS")
resoptLS=optimx( param, rmse, method=method, control=list(trace=0) )
rmse1=resoptLS$value
names(rmse1)=method
paramLS=t(resoptLS[,1:nparam, drop=0])
printm(param, paramLS, fmt="12.5")## param paramLS
## Nelder-Mead nlminb BFGS
## L1 1.00000 L1 0.24757 0.53002 0.53002
## L2 1.00000 L2 0.25328 0.51464 0.51464
## L3 1.00000 L3 0.19713 0.64477 0.64478
## G 1.00000 G 0.45558 0.65600 0.86614
## K1 1.00000 K1 0.37059 0.44874 0.33986
## K2 1.00000 K2 2.20509 0.71630 0.54251
## K3 1.00000 K3 -0.09728 0.42724 0.32357
## ve1 1.00000 ve1 1.16473 0.71907 0.71908
## ve2 1.00000 ve2 1.31027 0.73514 0.73514
## ve3 1.00000 ve3 1.07500 0.58427 0.58426
## ve4 1.00000 ve4 1.36845 0.79575 0.79575
## ve5 1.00000 ve5 2.31216 0.47956 0.47955
## ve6 1.00000 ve6 1.15334 0.81486 0.81486
## vzeta 1.00000 vzeta -0.47659 0.58397 1.01808printm(rmse0,t(rmse1), fmt="12.5")## rmse0 t(rmse1)
## Nelder-Mead nlminb BFGS
## 1.21851 [1,] 0.23208 0.02855 0.02855# gradient
gradLS=apply(paramLS,2, d1, rmse)
rownames(gradLS)=paramname
printm(gradLS, fmt="16.7")## gradLS
## Nelder-Mead nlminb BFGS
## L1 -0.0075303 0.0000000 -0.0000026
## L2 0.0250519 -0.0000001 -0.0000062
## L3 -0.0472360 -0.0000000 -0.0000055
## G -0.0771146 0.0000000 -0.0000129
## K1 0.0496193 0.0000001 -0.0000155
## K2 0.0081204 0.0000000 -0.0000145
## K3 0.0188118 0.0000000 -0.0000088
## ve1 0.0463769 -0.0000001 0.0000029
## ve2 0.0768260 0.0000001 -0.0000006
## ve3 0.0233627 -0.0000000 -0.0000069
## ve4 0.0680187 0.0000000 -0.0000039
## ve5 0.0008170 -0.0000000 -0.0000070
## ve6 0.0309407 -0.0000001 -0.0000016
## vzeta -0.0730150 0.0000001 -0.0000024# Hessian matrix
HessLS=d2(paramLS[,2], rmse)
dimnames(HessLS)=list(paramname,paramname)
printm(HessLS, fmt="6.3")## HessLS
## L1 L2 L3 G K1 K2 K3 ve1 ve2 ve3 ve4 ve5 ve6 vzeta
## L1 3.653 0.403 0.631 0.501 0.232 0.234 0.132 1.768 -0.000 0.000 -0.000 0.000 0.000 0.000
## L2 0.403 3.573 0.533 0.567 0.161 0.319 0.167 0.000 1.717 0.000 -0.000 0.000 0.000 0.000
## L3 0.631 0.533 4.327 0.599 0.242 0.308 0.148 -0.000 -0.000 2.151 0.000 0.000 0.000 0.000
## G 0.501 0.567 0.599 3.077 1.578 2.994 1.463 0.000 0.000 0.000 0.441 1.123 0.399 1.249
## K1 0.232 0.161 0.242 1.578 3.266 0.494 0.303 0.000 0.000 0.000 1.518 0.000 0.000 0.782
## K2 0.234 0.319 0.308 2.994 0.494 4.870 0.602 -0.000 0.000 -0.000 0.000 2.423 0.000 1.716
## K3 0.132 0.167 0.148 1.463 0.303 0.602 3.169 0.000 -0.000 0.000 0.000 -0.000 1.445 0.823
## ve1 1.768 0.000 -0.000 0.000 0.000 -0.000 0.000 1.668 -0.000 0.000 0.000 0.000 -0.000 0.000
## ve2 -0.000 1.717 -0.000 0.000 0.000 0.000 -0.000 -0.000 1.668 -0.000 0.000 0.000 -0.000 0.000
## ve3 0.000 0.000 2.151 0.000 0.000 -0.000 0.000 0.000 -0.000 1.668 0.000 0.000 0.000 0.000
## ve4 -0.000 -0.000 0.000 0.441 1.518 0.000 0.000 0.000 0.000 0.000 1.668 0.000 0.000 0.336
## ve5 0.000 0.000 0.000 1.123 0.000 2.423 -0.000 0.000 0.000 0.000 0.000 1.668 0.000 0.856
## ve6 0.000 0.000 0.000 0.399 0.000 0.000 1.445 -0.000 -0.000 0.000 0.000 0.000 1.668 0.304
## vzeta 0.000 0.000 0.000 1.249 0.782 1.716 0.823 0.000 0.000 0.000 0.336 0.856 0.304 0.952printm(nparam,qr(HessLS)$rank)## nparam qr(HessLS)$rank
## 14 14printm(eigen(HessLS)$values, fmt="10.5")## eigen(HessLS)$values
## 1 9.84680
## 2 5.85434
## 3 4.58037
## 4 4.36309
## 5 4.32672
## 6 3.88038
## 7 1.15785
## 8 0.76414
## 9 0.67756
## 10 0.53571
## 11 0.39389
## 12 0.32344
## 13 0.18944
## 14 0.00001QR 分解だとランクは落ちていないが、固有値で見ると怪しい。
3 RAM モデルによる PLS モデルの分析
ここまでは RAM による多重指標モデルの分析を行ってきたが、 そこでは全ての観測される変数は潜在変数の指標 (indicator) の役割を持つ 内生変数であった。
しかし、観測変数を外生変数として取り扱う共分散構造分析も存在する。 ここでは、図 1 の \(\xi\) から \(Y_1,Y_2,Y_3\) へのパスを反転させ、\(\xi\) に誤差を 加えたモデル（PLS モデル）の分析を行う。
図 2 PLS モデル
この場合、\(\boldsymbol v=\{ Y_j, j=1,2,\ldots,6, \; \xi, \eta \}\) であり、 \(\boldsymbol u=\{ Y_j, j=1,2,3, \; E_j, j=4,5,6, \; \zeta_1, \zeta_2 \}\) である。 ただし、\(\xi\) に加わる誤差を \(\zeta_1\)、\(\eta\) に加わる誤差を \(\zeta_2\) としたが、\(\mbox{${\bf\it var}\left(\zeta_1\right)$}=0\) と置くので結局 \(\xi\) に誤差は加わらない。
（\(\xi\) に誤差を加えるとパラメタが識別されない）
また、外生変数となった \(Y_j, j=1,2,3\) には誤差項を加えず、 変数間の共分散を仮定した。
固定パラメタとしては \(\lambda_1 = \kappa_1 = 1\) ならびに \(\mbox{${\bf\it var}\left(\zeta_1\right)$}=0\) を用いた。
options(width=200)
# variables
nvar=6
# variables of interest
v=c(paste("Y",1:nvar,sep=""),"xi", "eta" )
# Exogenous
u=c(paste("Y",1:3,sep=""),paste("E",4:6,sep=""),"zeta1", "zeta2")
nv=length(v)
locy=1:nvar
printm(nvar,nv,locy,v,u)## nvar nv locy v u
## 6 8 1 1 1 Y1 1 Y1
## 2 2 2 Y2 2 Y2
## 3 3 3 Y3 3 Y3
## 4 4 4 Y4 4 E4
## 5 5 5 Y5 5 E5
## 6 6 6 Y6 6 E6
## 7 xi 7 zeta1
## 8 eta 8 zeta2次に、\({{\cal D}}\left(\boldsymbol u\right)\) に含まれる分散共分散行列の要素のうち、 自由パラメタとなる部分に適当に名前を付ける。 分散の名前を v から始まるように指定すると、最尤法ではその対数がパラメタとなる。
また、分散共分散行列の要素のうち既知の部分は該当する数値を入れておく。 以下の例では、\(\zeta_1\) を実質的に無視するためのその分散を 0 としている。
options(width=200)
# dispersion of u
Du=Diag(c(paste("vy",1:3, sep=""),paste("ve",4:6, sep=""),"vzeta1","vzeta2"))
dimnames(Du)=list(u,u)
Du[1,2]=Du[2,1]="c21"; Du[1,3]=Du[3,1]="c31"; Du[2,3]=Du[3,2]="c32"
# fixed value
Du["zeta1","zeta1"]=0
printm(v,u,Du)## v u Du
## Y1 Y2 Y3 E4 E5 E6 zeta1 zeta2
## 1 Y1 1 Y1 Y1 vy1 c21 c31 0 0 0 0 0
## 2 Y2 2 Y2 Y2 c21 vy2 c32 0 0 0 0 0
## 3 Y3 3 Y3 Y3 c31 c32 vy3 0 0 0 0 0
## 4 Y4 4 E4 E4 0 0 0 ve4 0 0 0 0
## 5 Y5 5 E5 E5 0 0 0 0 ve5 0 0 0
## 6 Y6 6 E6 E6 0 0 0 0 0 ve6 0 0
## 7 xi 7 zeta1 zeta1 0 0 0 0 0 0 0 0
## 8 eta 8 zeta2 zeta2 0 0 0 0 0 0 0 vzeta2最後に \(\boldsymbol B\) に含まれる係数のうち自由パラメタとなる部分に適当に名前を付ける。 名前の付け方は自由であるが v で始まるものは避ける。
\(\boldsymbol B\) 行列の \(i,j\) 要素には \(v_j\) から \(v_i\) へのパスが存在するときに 何らかの名称もしくは数値が付与される。 ここでは \(v_j\) から \(v_i\) へのパスの 係数を \(ij\) とを付加した形で命名している。
以下の例では \(Y_1\) から \(\xi\) へのパス係数と \(\heta\) から \(Y_4\) へのパス係数を 1 に固定している。
options(width=200)
#
# path matrix
# B[i,j] is the name of the coefficient of the path from j to i.
B=matrix(0,nv,nv)
dimnames(B)=list(v,v)
B[8,7]="G"
B[7,1:3]=paste("L",1:3,sep="")
B[4:6,8]=paste(paste("K",1:3,sep=""),"",sep="")
# fix
B["xi","Y1"]=1
B["Y4","eta"]=1
printm(B)## B
## Y1 Y2 Y3 Y4 Y5 Y6 xi eta
## Y1 0 0 0 0 0 0 0 0
## Y2 0 0 0 0 0 0 0 0
## Y3 0 0 0 0 0 0 0 0
## Y4 0 0 0 0 0 0 0 1
## Y5 0 0 0 0 0 0 0 K2
## Y6 0 0 0 0 0 0 0 K3
## xi 1 L2 L3 0 0 0 0 0
## eta 0 0 0 0 0 0 G 0パス図に関する入力は以上で終わりである。
3.1 \(\boldsymbol B\) 行列と \({{\cal D}}\left(\boldsymbol u\right)\) 行列から記号的な \({{\cal D}}\left(\boldsymbol v\right)\) を算出する。
\(\boldsymbol B\) 行列が定義されると、 \((\boldsymbol I-\boldsymbol B)^{-1}\) の逆行列を記号的に求めて最終的に \({{\cal D}}\left(\boldsymbol v\right)\) を 記号的に定義する。
また、\(\boldsymbol B\boldsymbol v+\boldsymbol u\) や \((\boldsymbol I-\boldsymbol B)^{-1}\boldsymbol u\) が出力され、パス図との対応を確かめる ことが出来る。
options(width=120)
# calculate inv(I-B)
invImB=matSweep(diag(nv)-B) |> Expand(all=1)
printm(B,invImB,Du)## B invImB
## Y1 Y2 Y3 Y4 Y5 Y6 xi eta Y1 Y2 Y3 Y4 Y5 Y6 xi eta
## Y1 0 0 0 0 0 0 0 0 Y1 1 0 0 0 0 0 0 0
## Y2 0 0 0 0 0 0 0 0 Y2 0 1 0 0 0 0 0 0
## Y3 0 0 0 0 0 0 0 0 Y3 0 0 1 0 0 0 0 0
## Y4 0 0 0 0 0 0 0 1 Y4 G G*L2 G*L3 1 0 0 G 1
## Y5 0 0 0 0 0 0 0 K2 Y5 K2*G K2*G*L2 K2*G*L3 0 1 0 K2*G K2
## Y6 0 0 0 0 0 0 0 K3 Y6 K3*G K3*G*L2 K3*G*L3 0 0 1 K3*G K3
## xi 1 L2 L3 0 0 0 0 0 xi 1 L2 L3 0 0 0 1 0
## eta 0 0 0 0 0 0 G 0 eta G G*L2 G*L3 0 0 0 G 1
##
## Du
## Y1 Y2 Y3 E4 E5 E6 zeta1 zeta2
## Y1 vy1 c21 c31 0 0 0 0 0
## Y2 c21 vy2 c32 0 0 0 0 0
## Y3 c31 c32 vy3 0 0 0 0 0
## E4 0 0 0 ve4 0 0 0 0
## E5 0 0 0 0 ve5 0 0 0
## E6 0 0 0 0 0 ve6 0 0
## zeta1 0 0 0 0 0 0 0 0
## zeta2 0 0 0 0 0 0 0 vzeta2printm(B%*%v+u, invImB%*%u)## B %*% v + u invImB %*% u
## [,1] [,1]
## Y1 Y1 Y1 Y1
## Y2 Y2 Y2 Y2
## Y3 Y3 Y3 Y3
## Y4 eta+E4 Y4 G*Y1+G*L2*Y2+G*L3*Y3+E4+G*zeta1+zeta2
## Y5 K2*eta+E5 Y5 K2*G*Y1+K2*G*L2*Y2+K2*G*L3*Y3+E5+K2*G*zeta1+K2*zeta2
## Y6 K3*eta+E6 Y6 K3*G*Y1+K3*G*L2*Y2+K3*G*L3*Y3+E6+K3*G*zeta1+K3*zeta2
## xi Y1+L2*Y2+L3*Y3+zeta1 xi Y1+L2*Y2+L3*Y3+zeta1
## eta G*xi+zeta2 eta G*Y1+G*L2*Y2+G*L3*Y3+G*zeta1+zeta2# calculate D(v)=inv(I-B) D(u) inv(I-B)'
Dv=invImB%*%Du%*%t(invImB)
printm(Dv, .lj.=1)## Dv
## Y1
## Y1 vy1
## Y2 c21
## Y3 c31
## Y4 (G*vy1+G*L2*c21+G*L3*c31)
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)
## xi (vy1+L2*c21+L3*c31)
## eta (G*vy1+G*L2*c21+G*L3*c31)
## Y2
## Y1 c21
## Y2 vy2
## Y3 c32
## Y4 (G*c21+G*L2*vy2+G*L3*c32)
## Y5 (K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)
## Y6 (K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)
## xi (c21+L2*vy2+L3*c32)
## eta (G*c21+G*L2*vy2+G*L3*c32)
## Y3
## Y1 c31
## Y2 c32
## Y3 vy3
## Y4 (G*c31+G*L2*c32+G*L3*vy3)
## Y5 (K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)
## Y6 (K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)
## xi (c31+L2*c32+L3*vy3)
## eta (G*c31+G*L2*c32+G*L3*vy3)
## Y4
## Y1 vy1*G+c21*G*L2+c31*G*L3
## Y2 c21*G+vy2*G*L2+c32*G*L3
## Y3 c31*G+c32*G*L2+vy3*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+ve4+vzeta2
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*G*L3+K2*vzeta2
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*G*L3+K3*vzeta2
## xi (vy1+L2*c21+L3*c31)*G+(c21+L2*vy2+L3*c32)*G*L2+(c31+L2*c32+L3*vy3)*G*L3
## eta (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+vzeta2
## Y5
## Y1 vy1*K2*G+c21*K2*G*L2+c31*K2*G*L3
## Y2 c21*K2*G+vy2*K2*G*L2+c32*K2*G*L3
## Y3 c31*K2*G+c32*K2*G*L2+vy3*K2*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*K2*G+(G*c21+G*L2*vy2+G*L3*c32)*K2*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K2*G*L3+vzeta2*K2
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*K2*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*K2*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*K2*G*L3+ve5+K2*vzeta2*K2
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*K2*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*K2*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*K2*G*L3+K3*vzeta2*K2
## xi (vy1+L2*c21+L3*c31)*K2*G+(c21+L2*vy2+L3*c32)*K2*G*L2+(c31+L2*c32+L3*vy3)*K2*G*L3
## eta (G*vy1+G*L2*c21+G*L3*c31)*K2*G+(G*c21+G*L2*vy2+G*L3*c32)*K2*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K2*G*L3+vzeta2*K2
## Y6
## Y1 vy1*K3*G+c21*K3*G*L2+c31*K3*G*L3
## Y2 c21*K3*G+vy2*K3*G*L2+c32*K3*G*L3
## Y3 c31*K3*G+c32*K3*G*L2+vy3*K3*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*K3*G+(G*c21+G*L2*vy2+G*L3*c32)*K3*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K3*G*L3+vzeta2*K3
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*K3*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*K3*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*K3*G*L3+K2*vzeta2*K3
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*K3*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*K3*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*K3*G*L3+ve6+K3*vzeta2*K3
## xi (vy1+L2*c21+L3*c31)*K3*G+(c21+L2*vy2+L3*c32)*K3*G*L2+(c31+L2*c32+L3*vy3)*K3*G*L3
## eta (G*vy1+G*L2*c21+G*L3*c31)*K3*G+(G*c21+G*L2*vy2+G*L3*c32)*K3*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K3*G*L3+vzeta2*K3
## xi
## Y1 vy1+c21*L2+c31*L3
## Y2 c21+vy2*L2+c32*L3
## Y3 c31+c32*L2+vy3*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)+(G*c21+G*L2*vy2+G*L3*c32)*L2+(G*c31+G*L2*c32+G*L3*vy3)*L3
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*L3
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*L3
## xi (vy1+L2*c21+L3*c31)+(c21+L2*vy2+L3*c32)*L2+(c31+L2*c32+L3*vy3)*L3
## eta (G*vy1+G*L2*c21+G*L3*c31)+(G*c21+G*L2*vy2+G*L3*c32)*L2+(G*c31+G*L2*c32+G*L3*vy3)*L3
## eta
## Y1 vy1*G+c21*G*L2+c31*G*L3
## Y2 c21*G+vy2*G*L2+c32*G*L3
## Y3 c31*G+c32*G*L2+vy3*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+vzeta2
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*G*L3+K2*vzeta2
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*G*L3+K3*vzeta2
## xi (vy1+L2*c21+L3*c31)*G+(c21+L2*vy2+L3*c32)*G*L2+(c31+L2*c32+L3*vy3)*G*L3
## eta (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+vzeta2これで、必要となる \({{\cal D}}\left(\boldsymbol v\right)\) が定義できた。
3.2 パラメタベクトルの保管場所を定義する。
次に初期値となる数値的パラメタベクトル praram を定めるが、 param の各要素がどのパラメタの値であるかは、要素の名称 names(param) として 与えられていなければならない。
ここでは param の要素の名前のみを生成し、数値は仮の値を入れておく。 実際の初期値は、関数の最適化の直前に代入する。
options(width=200)
# Parameter vector
bname=unique( B[!sapply(B, is.number, pat=2)] )
dname=unique( Du[!sapply(Du, is.number, pat=2)] )
paramname=c(bname,dname)
nparam=length(paramname)
param=1:nparam
names(param)=paramname
printm(param)## param
## L2 1
## L3 2
## G 3
## K2 4
## K3 5
## vy1 6
## c21 7
## c31 8
## vy2 9
## c32 10
## vy3 11
## ve4 12
## ve5 13
## ve6 14
## vzeta2 153.3 R の最適化関数を用いるための準備
まず、先ほど作成した generate_Dv を用いてパラメタベクトルを与えたときに \({{\cal D}}\left(\boldsymbol v\right)\) を数値的に計算してくれる関数 get_Dv を作成する。 ただし、必要となるのは \({{\cal D}}\left(\boldsymbol y\right)\) の部分のみであるため、引数 Dv として Dv[locy,locy] を渡している。
options(width=200)
get_Dy <- generate_Dv( param, Dv[locy,locy], print=2 )##
## The following function was generated.
##
## function (param, nvar)
## {
## L2 = param["L2"]
## L3 = param["L3"]
## G = param["G"]
## K2 = param["K2"]
## K3 = param["K3"]
## vy1 = param["vy1"]
## c21 = param["c21"]
## c31 = param["c31"]
## vy2 = param["vy2"]
## c32 = param["c32"]
## vy3 = param["vy3"]
## ve4 = param["ve4"]
## ve5 = param["ve5"]
## ve6 = param["ve6"]
## vzeta2 = param["vzeta2"]
## Dvn = diag(nvar)
## Dvn[1] = vy1
## Dvn[2] = c21
## Dvn[3] = c31
## Dvn[4] = (G * vy1 + G * L2 * c21 + G * L3 * c31)
## Dvn[5] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31)
## Dvn[6] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31)
## Dvn[7] = c21
## Dvn[8] = vy2
## Dvn[9] = c32
## Dvn[10] = (G * c21 + G * L2 * vy2 + G * L3 * c32)
## Dvn[11] = (K2 * G * c21 + K2 * G * L2 * vy2 + K2 * G * L3 *
## c32)
## Dvn[12] = (K3 * G * c21 + K3 * G * L2 * vy2 + K3 * G * L3 *
## c32)
## Dvn[13] = c31
## Dvn[14] = c32
## Dvn[15] = vy3
## Dvn[16] = (G * c31 + G * L2 * c32 + G * L3 * vy3)
## Dvn[17] = (K2 * G * c31 + K2 * G * L2 * c32 + K2 * G * L3 *
## vy3)
## Dvn[18] = (K3 * G * c31 + K3 * G * L2 * c32 + K3 * G * L3 *
## vy3)
## Dvn[19] = vy1 * G + c21 * G * L2 + c31 * G * L3
## Dvn[20] = c21 * G + vy2 * G * L2 + c32 * G * L3
## Dvn[21] = c31 * G + c32 * G * L2 + vy3 * G * L3
## Dvn[22] = (G * vy1 + G * L2 * c21 + G * L3 * c31) * G + (G *
## c21 + G * L2 * vy2 + G * L3 * c32) * G * L2 + (G * c31 +
## G * L2 * c32 + G * L3 * vy3) * G * L3 + ve4 + vzeta2
## Dvn[23] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31) * G + (K2 * G * c21 + K2 * G * L2 * vy2 + K2 * G *
## L3 * c32) * G * L2 + (K2 * G * c31 + K2 * G * L2 * c32 +
## K2 * G * L3 * vy3) * G * L3 + K2 * vzeta2
## Dvn[24] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31) * G + (K3 * G * c21 + K3 * G * L2 * vy2 + K3 * G *
## L3 * c32) * G * L2 + (K3 * G * c31 + K3 * G * L2 * c32 +
## K3 * G * L3 * vy3) * G * L3 + K3 * vzeta2
## Dvn[25] = vy1 * K2 * G + c21 * K2 * G * L2 + c31 * K2 * G *
## L3
## Dvn[26] = c21 * K2 * G + vy2 * K2 * G * L2 + c32 * K2 * G *
## L3
## Dvn[27] = c31 * K2 * G + c32 * K2 * G * L2 + vy3 * K2 * G *
## L3
## Dvn[28] = (G * vy1 + G * L2 * c21 + G * L3 * c31) * K2 *
## G + (G * c21 + G * L2 * vy2 + G * L3 * c32) * K2 * G *
## L2 + (G * c31 + G * L2 * c32 + G * L3 * vy3) * K2 * G *
## L3 + vzeta2 * K2
## Dvn[29] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31) * K2 * G + (K2 * G * c21 + K2 * G * L2 * vy2 + K2 *
## G * L3 * c32) * K2 * G * L2 + (K2 * G * c31 + K2 * G *
## L2 * c32 + K2 * G * L3 * vy3) * K2 * G * L3 + ve5 + K2 *
## vzeta2 * K2
## Dvn[30] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31) * K2 * G + (K3 * G * c21 + K3 * G * L2 * vy2 + K3 *
## G * L3 * c32) * K2 * G * L2 + (K3 * G * c31 + K3 * G *
## L2 * c32 + K3 * G * L3 * vy3) * K2 * G * L3 + K3 * vzeta2 *
## K2
## Dvn[31] = vy1 * K3 * G + c21 * K3 * G * L2 + c31 * K3 * G *
## L3
## Dvn[32] = c21 * K3 * G + vy2 * K3 * G * L2 + c32 * K3 * G *
## L3
## Dvn[33] = c31 * K3 * G + c32 * K3 * G * L2 + vy3 * K3 * G *
## L3
## Dvn[34] = (G * vy1 + G * L2 * c21 + G * L3 * c31) * K3 *
## G + (G * c21 + G * L2 * vy2 + G * L3 * c32) * K3 * G *
## L2 + (G * c31 + G * L2 * c32 + G * L3 * vy3) * K3 * G *
## L3 + vzeta2 * K3
## Dvn[35] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31) * K3 * G + (K2 * G * c21 + K2 * G * L2 * vy2 + K2 *
## G * L3 * c32) * K3 * G * L2 + (K2 * G * c31 + K2 * G *
## L2 * c32 + K2 * G * L3 * vy3) * K3 * G * L3 + K2 * vzeta2 *
## K3
## Dvn[36] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31) * K3 * G + (K3 * G * c21 + K3 * G * L2 * vy2 + K3 *
## G * L3 * c32) * K3 * G * L2 + (K3 * G * c31 + K3 * G *
## L2 * c32 + K3 * G * L3 * vy3) * K3 * G * L3 + ve6 + K3 *
## vzeta2 * K3
## return(Dvn)
## }
## <environment: 0x000001a26b473240>
## with
## t(param)
## L2 L3 G K2 K3 vy1 c21 c31 vy2 c32 vy3 ve4 ve5 ve6 vzeta2
## [1,] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15以下は、前半と同様である。 rmse や mllh 等はそのまま利用できる。
3.4 R の最適化関数を利用したパラメタの推定
最小２乗解：
options(width=200)
# initial value of param: Names must be kept.
set.seed(1701)
param[]=rnorm(length(param),2,0.1)
param[]=rep(1,length(param))
rmse0=rmse(param)
# estimate parameters by optimx
method=c("Nelder-Mead","nlminb","BFGS")
resoptLS=optimx( param, rmse, method=method, control=list(trace=0) )
rmse1=resoptLS$value
names(rmse1)=method
paramLS=t(resoptLS[,1:nparam, drop=0])
printm(param, paramLS, fmt="12.5")## param paramLS
## Nelder-Mead nlminb BFGS
## L2 1.00000 L2 1.93360 0.49792 -7.37608
## L3 1.00000 L3 1.57796 1.51585 -17.81022
## G 1.00000 G 0.04458 0.09960 -0.01168
## K2 1.00000 K2 0.68317 1.61343 1.52806
## K3 1.00000 K3 1.16026 0.95925 0.92754
## vy1 1.00000 vy1 1.20813 1.00000 0.99634
## c21 1.00000 c21 0.42566 0.30400 0.33137
## c31 1.00000 c31 0.98373 0.30500 0.36971
## vy2 1.00000 vy2 1.48504 1.00000 0.99725
## c32 1.00000 c32 1.41361 0.34400 0.32477
## vy3 1.00000 vy3 0.93768 1.00000 0.97940
## ve4 1.00000 ve4 0.75202 0.79853 0.78522
## ve5 1.00000 ve5 0.77359 0.47555 0.49995
## ve6 1.00000 ve6 1.08901 0.81462 0.81553
## vzeta2 1.00000 vzeta2 0.22605 0.14896 0.15556printm(rmse0,t(rmse1), fmt="12.5")## rmse0 t(rmse1)
## Nelder-Mead nlminb BFGS
## 5.58830 [1,] 0.32822 0.02439 0.04533# gradient
gradLS=apply(paramLS,2, d1, rmse)
rownames(gradLS)=paramname
printm(gradLS, fmt="16.7")## gradLS
## Nelder-Mead nlminb BFGS
## L2 0.0040884 -0.0000000 0.0002428
## L3 0.0036617 -0.0000000 0.0000857
## G 0.3405288 -0.0000005 0.0192483
## K2 -0.0250110 0.0000001 -0.0008642
## K3 0.0510780 -0.0000001 -0.0008658
## vy1 0.0300373 0.0000002 -0.0001182
## c21 0.0192560 0.0000001 0.0006340
## c31 0.0985363 -0.0000001 0.0005167
## vy2 0.0740646 0.0000004 0.0017945
## c32 0.1588041 0.0000000 -0.0007018
## vy3 -0.0085448 -0.0000001 -0.0018361
## ve4 0.0032979 0.0000004 -0.0001772
## ve5 -0.0145167 0.0000003 0.0011068
## ve6 0.0657866 -0.0000001 0.0001777
## vzeta2 0.0794586 0.0000008 -0.0012452初期値で計算された rmse の値は 5.5883043 であり、最終的に得られた最小２乗解では その値が 0.3282239, 0.0243877, 0.0453311 にまで減少している。 数値的に計算したグレイディエントからみて、 nlminb がうまく収束しているように見える。
パラメタの２次微分の計算
HessLS=d2(paramLS[,2], rmse)
dimnames(HessLS)=list(paramname,paramname)
printm(HessLS, fmt="6.3")## HessLS
## L2 L3 G K2 K3 vy1 c21 c31 vy2 c32 vy3 ve4 ve5 ve6 vzeta2
## L2 0.126 0.098 2.985 0.200 0.098 0.034 0.108 0.093 0.045 0.159 0.063 0.051 0.133 0.047 0.746
## L3 0.098 0.151 3.769 0.271 0.131 0.038 0.055 0.162 0.018 0.106 0.158 0.077 0.201 0.071 1.122
## G 2.985 3.769 105.080 7.362 3.607 1.717 2.167 4.805 0.653 3.086 3.339 2.059 5.359 1.894 29.936
## K2 0.200 0.271 7.362 1.139 0.141 0.091 0.113 0.267 0.034 0.167 0.197 0.000 1.269 0.000 4.543
## K3 0.098 0.131 3.607 0.141 0.739 0.047 0.057 0.134 0.017 0.082 0.095 0.000 0.000 0.755 2.169
## vy1 0.034 0.038 1.717 0.091 0.047 2.043 0.046 0.141 0.001 0.004 0.006 0.019 0.050 0.018 0.282
## c21 0.108 0.055 2.167 0.113 0.057 0.046 2.065 0.075 0.044 0.137 0.006 0.019 0.050 0.018 0.281
## c31 0.093 0.162 4.805 0.267 0.134 0.141 0.075 2.267 0.002 0.056 0.152 0.059 0.153 0.054 0.854
## vy2 0.045 0.018 0.653 0.034 0.017 0.001 0.044 0.002 1.974 0.067 0.002 0.005 0.013 0.004 0.070
## c32 0.159 0.106 3.086 0.167 0.082 0.004 0.137 0.056 0.067 2.182 0.076 0.029 0.076 0.027 0.425
## vy3 0.063 0.158 3.339 0.197 0.095 0.006 0.006 0.152 0.002 0.076 2.169 0.045 0.116 0.041 0.647
## ve4 0.051 0.077 2.059 0.000 0.000 0.019 0.019 0.059 0.005 0.029 0.045 1.953 0.000 0.000 1.953
## ve5 0.133 0.201 5.359 1.269 0.000 0.050 0.050 0.153 0.013 0.076 0.116 0.000 1.953 0.000 5.083
## ve6 0.047 0.071 1.894 0.000 0.755 0.018 0.018 0.054 0.004 0.027 0.041 0.000 0.000 1.953 1.797
## vzeta2 0.746 1.122 29.936 4.543 2.169 0.282 0.281 0.854 0.070 0.425 0.647 1.953 5.083 1.797 28.394printm(nparam,qr(HessLS)$rank)## nparam qr(HessLS)$rank
## 15 15printm(eigen(HessLS)$values, fmt="10.5")## eigen(HessLS)$values
## 1 117.32801
## 2 19.11815
## 3 2.38630
## 4 2.23517
## 5 2.13070
## 6 2.03110
## 7 1.95260
## 8 1.95259
## 9 1.95258
## 10 1.93580
## 11 0.85741
## 12 0.19384
## 13 0.06168
## 14 0.03763
## 15 0.01119この場合、Hessian のランクも落ちていない。
最尤解：
初期値は最小２乗解を用いる。
options(width=200)
# initial from LS
param=unlist(resoptLS["nlminb",1:nparam])
# log(variance) as the paramrter
locv=which( substr(names(param),1,1) == "v" )
param[locv]=log(pmax(0.01,param[locv]))
mllh0=mllh(param)
# estimate parameters by optimx
options(warn=-1)
method=c("Nelder-Mead","nlminb","BFGS")
resoptML=optimx( param, mllh, method=method
, control=list(trace=0, maxit=500
, rel.tol=1e-15, step.min=1/16, x.tol=1e-15) )
options(warn=1)
resoptML[abs(resoptML) > 1e9]=1e9
paramML=t(resoptML[,1:nparam, drop=0])
paramML0=paramML
# back to original
locv=which( substr(names(param),1,1) == "v" )
paramML[locv,]=exp(paramML[locv,])
printm(param,paramML, fmt="12.5")## param paramML
## Nelder-Mead nlminb BFGS
## L2 0.49792 L2 0.48071 0.41898 0.41898
## L3 1.51585 L3 1.60640 1.43691 1.43691
## G 0.09960 G 0.09606 0.10811 0.10811
## K2 1.61343 K2 1.63611 1.57902 1.57902
## K3 0.95925 K3 0.93027 0.86238 0.86238
## vy1 0.00000 vy1 1.00152 1.00000 1.00000
## c21 0.30400 c21 0.30413 0.30400 0.30400
## c31 0.30500 c31 0.30409 0.30500 0.30500
## vy2 0.00000 vy2 1.00195 1.00000 1.00000
## c32 0.34400 c32 0.34320 0.34400 0.34400
## vy3 -0.00000 vy3 0.99974 1.00000 1.00000
## ve4 -0.22498 ve4 0.79258 0.78312 0.78312
## ve5 -0.74328 ve5 0.46453 0.45926 0.45926
## ve6 -0.20504 ve6 0.82785 0.83871 0.83871
## vzeta2 -1.90407 vzeta2 0.14954 0.16094 0.16094printm(mllh0,resoptML$value, fmt="12.5" )## mllh0 resoptML$value
## 0.02477 1 0.02434
## 2 0.02363
## 3 0.02363# gradient
gradML=apply(paramML0,2, d1, mllh)
rownames(gradML)=paramname
printm(gradML, fmt="16.7")## gradML
## Nelder-Mead nlminb BFGS
## L2 0.0009378 -0.0000000 0.0000001
## L3 0.0013419 -0.0000000 0.0000002
## G -0.0235709 0.0000000 0.0000043
## K2 0.0034398 -0.0000000 0.0000003
## K3 0.0088628 0.0000000 0.0000001
## vy1 0.0025199 0.0000000 -0.0000007
## c21 -0.0008779 -0.0000000 0.0000008
## c31 -0.0027841 0.0000000 -0.0000012
## vy2 0.0032346 0.0000000 -0.0000004
## c32 -0.0032035 0.0000000 -0.0000022
## vy3 0.0011631 0.0000000 0.0000003
## ve4 0.0042440 0.0000000 -0.0000001
## ve5 0.0020831 -0.0000000 -0.0000001
## ve6 -0.0056830 -0.0000000 0.0000002
## vzeta2 0.0021307 0.0000000 -0.0000001数値的に計算したグレイディエントから見て Nelder-Mead 以外は収束している。
パラメタの２次微分の計算
HessML=d2(paramML0[,2], mllh)
dimnames(HessML)=list(paramname,paramname)
printm(HessML, fmt="6.3")## HessML
## L2 L3 G K2 K3 vy1 c21 c31 vy2 c32 vy3 ve4 ve5 ve6 vzeta2
## L2 0.080 0.028 0.899 0.046 0.011 0.002 -0.002 0.000 0.001 0.000 -0.002 0.004 -0.005 0.002 0.001
## L3 0.028 0.080 1.394 0.068 0.013 -0.002 0.001 -0.001 0.002 -0.002 0.001 -0.005 -0.001 0.005 0.000
## G 0.899 1.394 32.707 1.577 0.362 0.000 0.000 0.001 0.000 0.000 -0.001 -0.142 0.040 0.103 -0.001
## K2 0.046 0.068 1.577 0.505 -0.027 -0.001 -0.001 -0.001 -0.001 0.000 -0.001 -0.038 0.303 -0.040 0.276
## K3 0.011 0.013 0.362 -0.027 0.333 0.001 0.002 -0.001 0.000 0.001 0.001 -0.011 -0.075 0.139 0.096
## vy1 0.002 -0.002 0.000 -0.001 0.001 1.346 -0.606 -0.611 0.068 0.138 0.070 0.001 0.000 0.002 -0.001
## c21 -0.002 0.001 0.000 -0.001 0.002 -0.606 2.906 -0.630 -0.625 -0.457 0.173 0.001 -0.001 0.001 -0.002
## c31 0.000 -0.001 0.001 -0.001 -0.001 -0.611 -0.630 2.911 0.173 -0.451 -0.631 0.000 0.000 0.000 0.000
## vy2 0.001 0.002 0.000 -0.001 0.000 0.068 -0.625 0.173 1.426 -0.789 0.110 0.000 0.000 -0.002 0.000
## c32 0.000 -0.002 0.000 0.000 0.001 0.138 -0.457 -0.451 -0.789 3.071 -0.790 -0.003 0.001 0.001 -0.001
## vy3 -0.002 0.001 -0.001 -0.001 0.001 0.070 0.173 -0.631 0.110 -0.790 1.428 0.000 0.000 0.000 0.000
## ve4 0.004 -0.005 -0.142 -0.038 -0.011 0.001 0.001 0.000 0.000 -0.003 0.000 0.825 0.025 0.002 0.057
## ve5 -0.005 -0.001 0.040 0.303 -0.075 0.000 -0.001 0.000 0.000 0.001 0.000 0.025 0.369 0.043 0.171
## ve6 0.002 0.005 0.103 -0.040 0.139 0.002 0.001 0.000 -0.002 0.001 0.000 0.002 0.043 0.876 0.020
## vzeta2 0.001 0.000 -0.001 0.276 0.096 -0.001 -0.002 0.000 0.000 -0.001 0.000 0.057 0.171 0.020 0.303printm(nparam,qr(HessML)$rank)## nparam qr(HessML)$rank
## 15 15printm(eigen(HessML)$values, fmt="10.5")## eigen(HessML)$values
## 1 32.87338
## 2 3.80560
## 3 3.79633
## 4 2.87853
## 5 1.05810
## 6 1.00252
## 7 0.91449
## 8 0.88780
## 9 0.82225
## 10 0.54720
## 11 0.36747
## 12 0.10156
## 13 0.05770
## 14 0.03549
## 15 0.01758Hessian はランク落ちしていない。
3.5 識別されないモデル。
以下に \(\xi\) に誤差 \(\zeta_1\) を加えた場合の例を示す。 変更する部分は Du 行列の \(\zeta_1\) の分散の定義部分と、 それに続くパラメタの定義部分、 ならびに、パラメタを与えたときに \({{\cal D}}\left(\boldsymbol v\right)\) を計算する関数の生成部分である。
モデルが識別されない場合、 最適化される関数のパラメタベクトルによる２次微分の行列、 すなわち Hessian 行列が正則ではなくなる。
options(width=200)
# dispersion of u
Du=Diag(c(paste("vy",1:3, sep=""),paste("ve",4:6, sep=""),"vzeta1","vzeta2"))
dimnames(Du)=list(u,u)
Du[1,2]=Du[2,1]="c21"; Du[1,3]=Du[3,1]="c31"; Du[2,3]=Du[3,2]="c32"
# do NOT fix
### Du["zeta1","zeta1"]=0
printm(v,u,Du)## v u Du
## Y1 Y2 Y3 E4 E5 E6 zeta1 zeta2
## 1 Y1 1 Y1 Y1 vy1 c21 c31 0 0 0 0 0
## 2 Y2 2 Y2 Y2 c21 vy2 c32 0 0 0 0 0
## 3 Y3 3 Y3 Y3 c31 c32 vy3 0 0 0 0 0
## 4 Y4 4 E4 E4 0 0 0 ve4 0 0 0 0
## 5 Y5 5 E5 E5 0 0 0 0 ve5 0 0 0
## 6 Y6 6 E6 E6 0 0 0 0 0 ve6 0 0
## 7 xi 7 zeta1 zeta1 0 0 0 0 0 0 vzeta1 0
## 8 eta 8 zeta2 zeta2 0 0 0 0 0 0 0 vzeta2# same B matrix will be used.
# calculate D(v)=inv(I-B) D(u) inv(I-B)'
Dv=invImB%*%Du%*%t(invImB)
printm(Dv, .lj.=1)## Dv
## Y1
## Y1 vy1
## Y2 c21
## Y3 c31
## Y4 (G*vy1+G*L2*c21+G*L3*c31)
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)
## xi (vy1+L2*c21+L3*c31)
## eta (G*vy1+G*L2*c21+G*L3*c31)
## Y2
## Y1 c21
## Y2 vy2
## Y3 c32
## Y4 (G*c21+G*L2*vy2+G*L3*c32)
## Y5 (K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)
## Y6 (K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)
## xi (c21+L2*vy2+L3*c32)
## eta (G*c21+G*L2*vy2+G*L3*c32)
## Y3
## Y1 c31
## Y2 c32
## Y3 vy3
## Y4 (G*c31+G*L2*c32+G*L3*vy3)
## Y5 (K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)
## Y6 (K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)
## xi (c31+L2*c32+L3*vy3)
## eta (G*c31+G*L2*c32+G*L3*vy3)
## Y4
## Y1 vy1*G+c21*G*L2+c31*G*L3
## Y2 c21*G+vy2*G*L2+c32*G*L3
## Y3 c31*G+c32*G*L2+vy3*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+ve4+G*vzeta1*G+vzeta2
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*G*L3+K2*G*vzeta1*G+K2*vzeta2
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*G*L3+K3*G*vzeta1*G+K3*vzeta2
## xi (vy1+L2*c21+L3*c31)*G+(c21+L2*vy2+L3*c32)*G*L2+(c31+L2*c32+L3*vy3)*G*L3+vzeta1*G
## eta (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+G*vzeta1*G+vzeta2
## Y5
## Y1 vy1*K2*G+c21*K2*G*L2+c31*K2*G*L3
## Y2 c21*K2*G+vy2*K2*G*L2+c32*K2*G*L3
## Y3 c31*K2*G+c32*K2*G*L2+vy3*K2*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*K2*G+(G*c21+G*L2*vy2+G*L3*c32)*K2*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K2*G*L3+G*vzeta1*K2*G+vzeta2*K2
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*K2*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*K2*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*K2*G*L3+ve5+K2*G*vzeta1*K2*G+K2*vzeta2*K2
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*K2*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*K2*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*K2*G*L3+K3*G*vzeta1*K2*G+K3*vzeta2*K2
## xi (vy1+L2*c21+L3*c31)*K2*G+(c21+L2*vy2+L3*c32)*K2*G*L2+(c31+L2*c32+L3*vy3)*K2*G*L3+vzeta1*K2*G
## eta (G*vy1+G*L2*c21+G*L3*c31)*K2*G+(G*c21+G*L2*vy2+G*L3*c32)*K2*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K2*G*L3+G*vzeta1*K2*G+vzeta2*K2
## Y6
## Y1 vy1*K3*G+c21*K3*G*L2+c31*K3*G*L3
## Y2 c21*K3*G+vy2*K3*G*L2+c32*K3*G*L3
## Y3 c31*K3*G+c32*K3*G*L2+vy3*K3*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*K3*G+(G*c21+G*L2*vy2+G*L3*c32)*K3*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K3*G*L3+G*vzeta1*K3*G+vzeta2*K3
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*K3*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*K3*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*K3*G*L3+K2*G*vzeta1*K3*G+K2*vzeta2*K3
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*K3*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*K3*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*K3*G*L3+ve6+K3*G*vzeta1*K3*G+K3*vzeta2*K3
## xi (vy1+L2*c21+L3*c31)*K3*G+(c21+L2*vy2+L3*c32)*K3*G*L2+(c31+L2*c32+L3*vy3)*K3*G*L3+vzeta1*K3*G
## eta (G*vy1+G*L2*c21+G*L3*c31)*K3*G+(G*c21+G*L2*vy2+G*L3*c32)*K3*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*K3*G*L3+G*vzeta1*K3*G+vzeta2*K3
## xi
## Y1 vy1+c21*L2+c31*L3
## Y2 c21+vy2*L2+c32*L3
## Y3 c31+c32*L2+vy3*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)+(G*c21+G*L2*vy2+G*L3*c32)*L2+(G*c31+G*L2*c32+G*L3*vy3)*L3+G*vzeta1
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*L3+K2*G*vzeta1
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*L3+K3*G*vzeta1
## xi (vy1+L2*c21+L3*c31)+(c21+L2*vy2+L3*c32)*L2+(c31+L2*c32+L3*vy3)*L3+vzeta1
## eta (G*vy1+G*L2*c21+G*L3*c31)+(G*c21+G*L2*vy2+G*L3*c32)*L2+(G*c31+G*L2*c32+G*L3*vy3)*L3+G*vzeta1
## eta
## Y1 vy1*G+c21*G*L2+c31*G*L3
## Y2 c21*G+vy2*G*L2+c32*G*L3
## Y3 c31*G+c32*G*L2+vy3*G*L3
## Y4 (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+G*vzeta1*G+vzeta2
## Y5 (K2*G*vy1+K2*G*L2*c21+K2*G*L3*c31)*G+(K2*G*c21+K2*G*L2*vy2+K2*G*L3*c32)*G*L2+(K2*G*c31+K2*G*L2*c32+K2*G*L3*vy3)*G*L3+K2*G*vzeta1*G+K2*vzeta2
## Y6 (K3*G*vy1+K3*G*L2*c21+K3*G*L3*c31)*G+(K3*G*c21+K3*G*L2*vy2+K3*G*L3*c32)*G*L2+(K3*G*c31+K3*G*L2*c32+K3*G*L3*vy3)*G*L3+K3*G*vzeta1*G+K3*vzeta2
## xi (vy1+L2*c21+L3*c31)*G+(c21+L2*vy2+L3*c32)*G*L2+(c31+L2*c32+L3*vy3)*G*L3+vzeta1*G
## eta (G*vy1+G*L2*c21+G*L3*c31)*G+(G*c21+G*L2*vy2+G*L3*c32)*G*L2+(G*c31+G*L2*c32+G*L3*vy3)*G*L3+G*vzeta1*G+vzeta2#
# Parameter vector
bname=unique( B[!sapply(B, is.number, pat=2)] )
dname=unique( Du[!sapply(Du, is.number, pat=2)] )
paramname=c(bname,dname)
nparam=length(paramname)
param=1:nparam
names(param)=paramname
printm(param)## param
## L2 1
## L3 2
## G 3
## K2 4
## K3 5
## vy1 6
## c21 7
## c31 8
## vy2 9
## c32 10
## vy3 11
## ve4 12
## ve5 13
## ve6 14
## vzeta1 15
## vzeta2 16# new dispersion function
get_Dy <- generate_Dv( param, Dv[locy,locy], print=2 )##
## The following function was generated.
##
## function (param, nvar)
## {
## L2 = param["L2"]
## L3 = param["L3"]
## G = param["G"]
## K2 = param["K2"]
## K3 = param["K3"]
## vy1 = param["vy1"]
## c21 = param["c21"]
## c31 = param["c31"]
## vy2 = param["vy2"]
## c32 = param["c32"]
## vy3 = param["vy3"]
## ve4 = param["ve4"]
## ve5 = param["ve5"]
## ve6 = param["ve6"]
## vzeta1 = param["vzeta1"]
## vzeta2 = param["vzeta2"]
## Dvn = diag(nvar)
## Dvn[1] = vy1
## Dvn[2] = c21
## Dvn[3] = c31
## Dvn[4] = (G * vy1 + G * L2 * c21 + G * L3 * c31)
## Dvn[5] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31)
## Dvn[6] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31)
## Dvn[7] = c21
## Dvn[8] = vy2
## Dvn[9] = c32
## Dvn[10] = (G * c21 + G * L2 * vy2 + G * L3 * c32)
## Dvn[11] = (K2 * G * c21 + K2 * G * L2 * vy2 + K2 * G * L3 *
## c32)
## Dvn[12] = (K3 * G * c21 + K3 * G * L2 * vy2 + K3 * G * L3 *
## c32)
## Dvn[13] = c31
## Dvn[14] = c32
## Dvn[15] = vy3
## Dvn[16] = (G * c31 + G * L2 * c32 + G * L3 * vy3)
## Dvn[17] = (K2 * G * c31 + K2 * G * L2 * c32 + K2 * G * L3 *
## vy3)
## Dvn[18] = (K3 * G * c31 + K3 * G * L2 * c32 + K3 * G * L3 *
## vy3)
## Dvn[19] = vy1 * G + c21 * G * L2 + c31 * G * L3
## Dvn[20] = c21 * G + vy2 * G * L2 + c32 * G * L3
## Dvn[21] = c31 * G + c32 * G * L2 + vy3 * G * L3
## Dvn[22] = (G * vy1 + G * L2 * c21 + G * L3 * c31) * G + (G *
## c21 + G * L2 * vy2 + G * L3 * c32) * G * L2 + (G * c31 +
## G * L2 * c32 + G * L3 * vy3) * G * L3 + ve4 + G * vzeta1 *
## G + vzeta2
## Dvn[23] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31) * G + (K2 * G * c21 + K2 * G * L2 * vy2 + K2 * G *
## L3 * c32) * G * L2 + (K2 * G * c31 + K2 * G * L2 * c32 +
## K2 * G * L3 * vy3) * G * L3 + K2 * G * vzeta1 * G + K2 *
## vzeta2
## Dvn[24] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31) * G + (K3 * G * c21 + K3 * G * L2 * vy2 + K3 * G *
## L3 * c32) * G * L2 + (K3 * G * c31 + K3 * G * L2 * c32 +
## K3 * G * L3 * vy3) * G * L3 + K3 * G * vzeta1 * G + K3 *
## vzeta2
## Dvn[25] = vy1 * K2 * G + c21 * K2 * G * L2 + c31 * K2 * G *
## L3
## Dvn[26] = c21 * K2 * G + vy2 * K2 * G * L2 + c32 * K2 * G *
## L3
## Dvn[27] = c31 * K2 * G + c32 * K2 * G * L2 + vy3 * K2 * G *
## L3
## Dvn[28] = (G * vy1 + G * L2 * c21 + G * L3 * c31) * K2 *
## G + (G * c21 + G * L2 * vy2 + G * L3 * c32) * K2 * G *
## L2 + (G * c31 + G * L2 * c32 + G * L3 * vy3) * K2 * G *
## L3 + G * vzeta1 * K2 * G + vzeta2 * K2
## Dvn[29] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31) * K2 * G + (K2 * G * c21 + K2 * G * L2 * vy2 + K2 *
## G * L3 * c32) * K2 * G * L2 + (K2 * G * c31 + K2 * G *
## L2 * c32 + K2 * G * L3 * vy3) * K2 * G * L3 + ve5 + K2 *
## G * vzeta1 * K2 * G + K2 * vzeta2 * K2
## Dvn[30] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31) * K2 * G + (K3 * G * c21 + K3 * G * L2 * vy2 + K3 *
## G * L3 * c32) * K2 * G * L2 + (K3 * G * c31 + K3 * G *
## L2 * c32 + K3 * G * L3 * vy3) * K2 * G * L3 + K3 * G *
## vzeta1 * K2 * G + K3 * vzeta2 * K2
## Dvn[31] = vy1 * K3 * G + c21 * K3 * G * L2 + c31 * K3 * G *
## L3
## Dvn[32] = c21 * K3 * G + vy2 * K3 * G * L2 + c32 * K3 * G *
## L3
## Dvn[33] = c31 * K3 * G + c32 * K3 * G * L2 + vy3 * K3 * G *
## L3
## Dvn[34] = (G * vy1 + G * L2 * c21 + G * L3 * c31) * K3 *
## G + (G * c21 + G * L2 * vy2 + G * L3 * c32) * K3 * G *
## L2 + (G * c31 + G * L2 * c32 + G * L3 * vy3) * K3 * G *
## L3 + G * vzeta1 * K3 * G + vzeta2 * K3
## Dvn[35] = (K2 * G * vy1 + K2 * G * L2 * c21 + K2 * G * L3 *
## c31) * K3 * G + (K2 * G * c21 + K2 * G * L2 * vy2 + K2 *
## G * L3 * c32) * K3 * G * L2 + (K2 * G * c31 + K2 * G *
## L2 * c32 + K2 * G * L3 * vy3) * K3 * G * L3 + K2 * G *
## vzeta1 * K3 * G + K2 * vzeta2 * K3
## Dvn[36] = (K3 * G * vy1 + K3 * G * L2 * c21 + K3 * G * L3 *
## c31) * K3 * G + (K3 * G * c21 + K3 * G * L2 * vy2 + K3 *
## G * L3 * c32) * K3 * G * L2 + (K3 * G * c31 + K3 * G *
## L2 * c32 + K3 * G * L3 * vy3) * K3 * G * L3 + ve6 + K3 *
## G * vzeta1 * K3 * G + K3 * vzeta2 * K3
## return(Dvn)
## }
## <environment: 0x000001a2629d8bf8>
## with
## t(param)
## L2 L3 G K2 K3 vy1 c21 c31 vy2 c32 vy3 ve4 ve5 ve6 vzeta1 vzeta2
## [1,] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16# initial value of param: Names must be kept.
set.seed(1701)
param[]=rnorm(length(param),2,0.1)
param[]=rep(1,length(param))
rmse0=rmse(param)
# estimate parameters by optimx
method=c("Nelder-Mead","nlminb","BFGS")
resoptLS=optimx( param, rmse, method=method, control=list(trace=0) )
rmse1=resoptLS$value
names(rmse1)=method
paramLS=t(resoptLS[,1:nparam, drop=0])
printm(param, paramLS, fmt="12.5")## param paramLS
## Nelder-Mead nlminb BFGS
## L2 1.00000 L2 1.74563 0.49792 -12.94722
## L3 1.00000 L3 1.34090 1.51585 -21.26167
## G 1.00000 G 0.03549 0.09960 -0.00844
## K2 1.00000 K2 1.02115 1.61343 1.55324
## K3 1.00000 K3 0.06961 0.95925 0.93753
## vy1 1.00000 vy1 1.25816 1.00000 0.99698
## c21 1.00000 c21 0.76894 0.30400 0.33864
## c31 1.00000 c31 1.52992 0.30500 0.35973
## vy2 1.00000 vy2 1.25549 1.00000 0.98670
## c32 1.00000 c32 0.79858 0.34400 0.31935
## vy3 1.00000 vy3 1.41559 1.00000 0.99421
## ve4 1.00000 ve4 0.87637 0.79853 0.79111
## ve5 1.00000 ve5 0.67558 0.47555 0.49697
## ve6 1.00000 ve6 1.17968 0.81462 0.81653
## vzeta1 1.00000 vzeta1 1.28575 0.82985 0.49213
## vzeta2 1.00000 vzeta2 0.56325 0.14073 0.15399printm(rmse0,t(rmse1), fmt="12.5")## rmse0 t(rmse1)
## Nelder-Mead nlminb BFGS
## 6.09505 [1,] 0.36690 0.02439 0.04516# gradient
gradLS=apply(paramLS,2, d1, rmse)
rownames(gradLS)=paramname
printm(gradLS, fmt="16.7")## gradLS
## Nelder-Mead nlminb BFGS
## L2 0.0001059 0.0000000 -0.0001273
## L3 -0.0001586 0.0000000 0.0004462
## G 0.0050943 0.0000002 0.6564093
## K2 0.0554323 0.0000000 -0.0025250
## K3 -0.0389696 -0.0000000 -0.0022957
## vy1 0.0333710 0.0000001 -0.0004620
## c21 0.0600733 -0.0000000 0.0002195
## c31 0.1583034 0.0000001 -0.0004966
## vy2 0.0330982 0.0000002 -0.0010024
## c32 0.0580892 -0.0000001 -0.0015363
## vy3 0.0532754 0.0000001 -0.0010559
## ve4 0.0602357 -0.0000001 -0.0002475
## ve5 0.0374367 0.0000001 0.0003867
## ve6 0.0236902 -0.0000001 -0.0000776
## vzeta1 0.0001608 0.0000000 -0.0000009
## vzeta2 0.1276311 0.0000001 -0.0127591# Hessian matrix
HessLS=d2(paramLS[,2], rmse)
dimnames(HessLS)=list(paramname,paramname)
printm(HessLS, fmt="6.3")## HessLS
## L2 L3 G K2 K3 vy1 c21 c31 vy2 c32 vy3 ve4 ve5 ve6 vzeta1 vzeta2
## L2 0.126 0.098 3.108 0.200 0.098 0.034 0.108 0.093 0.045 0.159 0.063 0.051 0.133 0.047 0.007 0.746
## L3 0.098 0.150 3.955 0.271 0.131 0.038 0.055 0.162 0.018 0.106 0.158 0.077 0.201 0.071 0.011 1.122
## G 3.108 3.955 115.753 8.113 3.966 1.763 2.213 4.946 0.665 3.156 3.446 2.381 6.199 2.191 0.344 34.630
## K2 0.200 0.271 8.113 1.139 0.141 0.091 0.113 0.267 0.034 0.167 0.197 0.000 1.269 0.000 0.045 4.543
## K3 0.098 0.131 3.966 0.141 0.739 0.047 0.057 0.134 0.017 0.082 0.095 0.000 0.000 0.755 0.022 2.169
## vy1 0.034 0.038 1.763 0.091 0.047 2.043 0.046 0.141 0.001 0.004 0.006 0.019 0.050 0.018 0.003 0.282
## c21 0.108 0.055 2.213 0.113 0.057 0.046 2.065 0.075 0.044 0.137 0.006 0.019 0.050 0.018 0.003 0.281
## c31 0.093 0.162 4.946 0.267 0.134 0.141 0.075 2.267 0.002 0.056 0.152 0.059 0.153 0.054 0.008 0.854
## vy2 0.045 0.018 0.665 0.034 0.017 0.001 0.044 0.002 1.974 0.067 0.002 0.005 0.013 0.004 0.001 0.070
## c32 0.159 0.106 3.156 0.167 0.082 0.004 0.137 0.056 0.067 2.182 0.076 0.029 0.076 0.027 0.004 0.425
## vy3 0.063 0.158 3.446 0.197 0.095 0.006 0.006 0.152 0.002 0.076 2.169 0.045 0.116 0.041 0.006 0.647
## ve4 0.051 0.077 2.381 0.000 0.000 0.019 0.019 0.059 0.005 0.029 0.045 1.953 -0.000 0.000 0.019 1.953
## ve5 0.133 0.201 6.199 1.269 0.000 0.050 0.050 0.153 0.013 0.076 0.116 -0.000 1.953 0.000 0.050 5.083
## ve6 0.047 0.071 2.191 0.000 0.755 0.018 0.018 0.054 0.004 0.027 0.041 0.000 0.000 1.953 0.018 1.797
## vzeta1 0.007 0.011 0.344 0.045 0.022 0.003 0.003 0.008 0.001 0.004 0.006 0.019 0.050 0.018 0.003 0.282
## vzeta2 0.746 1.122 34.630 4.543 2.169 0.282 0.281 0.854 0.070 0.425 0.647 1.953 5.083 1.797 0.282 28.394printm(nparam,qr(HessLS)$rank)## nparam qr(HessLS)$rank
## 16 16printm(eigen(HessLS)$values, fmt="10.5")## eigen(HessLS)$values
## 1 129.84598
## 2 17.27893
## 3 2.38634
## 4 2.23517
## 5 2.13067
## 6 2.03126
## 7 1.95259
## 8 1.95258
## 9 1.95255
## 10 1.93295
## 11 0.85737
## 12 0.19388
## 13 0.06175
## 14 0.03761
## 15 0.01116
## 16 -0.00001QR 分解のランクは落ちていないが固有値を見ると怪しい。
knitr::knit_exit()