はじめに
東京大学/株式会社Nospare リサーチャーの栗栖です．
前回の記事で，最近の研究成果「大域的フレシェ回帰に対するモデル平均」 (Kurisu and Otsu (2025)) について理論的結果を紹介しました．今回の記事では数値実験の結果を紹介します．また記事の最後ではフレシェ回帰をRで実行する方法についても紹介します．
大域的フレシェ回帰のモデル平均
まず大域的フレシェ回帰のモデル平均の復習から始めます．$(\mathbb{Y},d)$を距離空間とし，独立同分布なデータ$\{(Y_i,X_i)\}_{i=1}^{n}$, $Y_i \in \mathbb{Y}$, $X_i \in \mathbb{R}^p$が利用可能であるとします．また予測子(predictor) の組 $\boldsymbol{X}_m:=(X_1,\dots, X_{k_m}) \in \mathbb{R}^{k_m}$, $1 \leq k_1 < k_2 < \dots < k_M$について，各組に対応するフレシェ回帰関数を$L_\oplus^{(m)}$, $m=1,\dots,M$とします．フレシェ回帰関数の定義については前回の記事を参照してください．
以上の設定のもと，大域的フレシェ回帰関数$L_\oplus^{(m)}$, $m=1,\dots,M$のモデル平均は以下で定義されます．
\hat{m}_{\oplus}(\hat{{\bf w}},x),\ \hat{\bf w} = \text{argmin}_{{\bf w} \in \mathbb{W}}\text{CV}_n({\bf w}).
ここで
\begin{align*}
\hat{m}_\oplus({\bf w},x) &= \text{argmin}_{\eta \in \mathbb{Y}}\sum_{m=1}^Mw_md^2(\hat{L}_\oplus^{(m)}(x),\eta),\\
\mathbb{W} &= \{{\bf w}=(w_1,\dots,w_M)' \in [0,1]^M: \sum_{m=1}^Mw_m=1\},\\
\text{CV}_n({\bf w}) &= {1 \over n}\sum_{i=1}^{n}d^2(Y_i,\hat{m}_{\oplus,-i}({\bf w},X_i)),\\
\end{align*}
です．$\hat{L}_\oplus^{(m)}(x)$や$\hat{m}_{\oplus,-i}({\bf w},x)$の定義については前回の記事を参照してください．
数値実験
以下では数値実験で提案手法のパフォーマンスを確認してみます．ここでは距離空間として，$5 \times 5$の対称正定値行列の空間$\mathcal{S}_5$を考え，コレスキー分解による距離を導入します．$A_1,A_2 \in \mathcal{S}_5$とすると，コレスキー分解により$A_j = (A_j^{1/2})'A_j^{1/2}$, $j=1,2$と表現できます．ここで，$A_j^{1/2}$は上三角かつ正の対角成分を持つ行列です．この分解に対してコレスキー分解距離を以下で定義します．
d_C(A_1,A_2) = \sqrt{\text{trace}((A_1^{1/2} - A_2^{1/2})'(A_1^{1/2} - A_2^{1/2}))}.
ここで$\text{trace}(A)$は行列$A$のトレースです．
データ生成過程
まず予測子として9次元のベクトル$X_i = (X_{i,1},\dots,X_{i,9})'$を考えます．ここで，予測子の各成分は以下のように9次元の多変量正規分布に従う確率ベクトルを変換した関数から生成されているとします：
$Z_i = (Z_{i,1},\dots, Z_{i,9})'$を平均$0$の多変量正規分布とし，各成分の相関を$\text{Cov}(Z_{i,j},Z_{i,k}) = 2^{-|j-k|}$となるように設定します．これらの値を用いて$X_{i,j} = 2\Phi(Z_{i,j})$とします．ここで$\Phi(\cdot)$は標準正規分布の累積分布関数です．
このように設定した予測子は各成分が相関を持ちます．次に応答変数$Y_i$を以下のように生成します．
Y = A'A,\ A = (\mu + \sigma)I_5 + \sigma V.
ここで，$I_5$は$5 \times 5$の単位行列，$V=(v_{j,k})_{1 \leq j,k \leq 5}$は$j < k$なら$v_{j,k}=1$，$j \geq k$なら$v_{j,k} = 0$となるような行列です．また$\mu, \sigma$はそれぞれ予測子$X$を与えたときの分布が正規分布，ガンマ分布に従う確率変数で，以下のように設定します．
\begin{align*}
\mu|X &\sim N\left(3 + 2\sum_{j=1}^5{X_{2j-1} \over 2j-1}, 1\right),\\
\sigma|X &\sim \text{Gamma}\left({1 \over 2}\left(\sum_{j=1}^4{X_{2j} \over 2j}\right)^2, {2 \over 3}{1 \over 1 + \left(\sum_{j=1}^4{X_{2j} \over 2j}\right)}\right).
\end{align*}
結果
サンプルサイズ$n=50$とし，モデル平均を考える予測子の候補として，以下の6つのモデルを考えます．
\text{M}_m: X_i^{(m)} = (X_{i,1},\dots,X_{i,m})',\ m = 1,...,6.
パフォーマンスの評価として，同じモデルから発生させた独立なデータ$\{X_s,Y_s\}_{s=1}^{100}$を用いて, 以下のように$\hat{m}_\oplus(\hat{\bf w},x)$の final prediction error (FPE) (定義については前回の記事を参照してください)を100回繰り返し計算して平均をとった値 ($\text{FPE}$)を比較します．
\text{FPE} = {1 \over 100}\sum_{r=1}^{100}\text{FPE}(r).
ここで，$\text{FPE}(r)$は以下で定義される値です．
\text{FPE}(r) = {1 \over 100}\sum_{r=1}^{100}d_C^2(Y_s, \hat{m}_\oplus(\hat{w},X_s))
下の図の横軸は左から$\text{M}_1$のみの場合，$\text{M}_1$と$\text{M}_2$を組み合わせる場合,...,$\text{M}_1$から$\text{M}_6$のモデルを組み合わせる場合のモデル平均推定量に対応し，縦軸はそれぞれの場合における$\text{FPE}$の値です．図の中でCVに対応する青丸が今回我々が提案した手法の$\text{FPE}$の値で，赤三角や緑十字は我々の提案手法の対抗手法としたその他のモデル平均の方法です．結果を見ると，2つ以上のモデルを組み合わせる場合はどの場合も提案手法の$\text{FPE}$が一番小さくなっていることがわかります．対抗手法の詳細や応答変数が確率分布の場合など，その他の設定での数値実験についてはKurisu and Otsu (2025)を参照してください．
以下では，$\text{M}_1, \text{M}_2, \text{M}_3$の３つのモデルを組み合わせる場合に上記の手続きを実行するためのRのコードの例を紹介します．
まずは上記のモデルからデータを生成する関数を準備します．
#データの生成
library(psych)
library(matlab)
library(ramify)
gendata=function(n,p,m,r) {
#' @param n: サンプルサイズ
#' @param p: 予測子の数
#' @param m: 行列のサイズ
#' @param r: 予測子の相関
x=2*pnorm(matrix(rnorm(n*p),n,p)%*%chol(r^(as.matrix(dist(1:p)))))
M <- array(0,c(m,m,n))
Vmu0=3
Vsigma0=3
Vbeta=2
Vgamma=3
Vnu1=1
Vnu2=2
Vmu=Vbeta*(x[,1] + x[,3]/3 + x[,5]/5 + x[,7]/7 + x[,9]/9)+rnorm(n, Vmu0,sqrt(Vnu1))
GamP1=(Vsigma0+Vgamma*(x[,2]/2 + x[,4]/4 + x[,6]/6 + x[,8]/8))^2/Vnu2
GamP2=Vnu2/(Vsigma0+Vgamma*(x[,2]/2 + x[,4]/4 + x[,6]/6 + x[,8]/8))
I =diag(m)
off_diag = matlab::ones(m)
off_diag[lower.tri(off_diag, diag = TRUE)] <- 0
for (j in 1:n) {
Vsigma=rgamma(1,shape=GamP1[j], scale=GamP2[j])
A = (Vmu[j]+Vsigma)*I + Vsigma*off_diag
aux<-t(A)%*%A
M[,,j]<-aux
}
return(list(x=x, Min=M))
}
次に，重みベクトルの候補を生成する関数を準備しておきます．
# 重みベクトルの候補の生成
weight_generate <- function(d, step) {
x <- seq(0, 1, by = step)
v <- expand.grid(replicate(d, x, simplify = FALSE))
th <- 1 + x[2]/10
v <- v[rowSums(v) <th & rowSums(v)>=1, ]
v <- matrix(unlist(v),ncol=d)
return(v)
}
さらに，与えられた重みベクトル${\bf w}=(w_1,w_2,w_3)'$の下でモデル平均推定量を計算する関数を準備しておきます．
GFRCovCholesky_w <- function(MC, Mest, w, optns = list()){
#' @param MC: 新たに発生させるデータの数
#' @param Mest: 各モデルの下でのフレシェ回帰の値
#' @param w: 重みベクトル
#' @param optns: 応答変数を相関行列とするかどうか，行列間の距離関数の選択
if(is.null(optns$corrOut)){
corrOut = FALSE
} else {
corrOut = optns$corrOut
}
if(is.null(optns$metric)){
metric = 'log_cholesky'
} else {
metric = optns$metric
}
m = nrow(Mest[,,1,1])
m_w_in = length(w)
M_JMA_w_in = array(dim=c(m,m,MC))
if(metric == 'log_cholesky'){
for (j in 1:MC) {
ss = 0
U = 0
E = 0
s = array(0,m_w_in)
for (i in 1:m_w_in) {
M = (Mest[,,i,j] +t(Mest[,,i,j]))/2
LL = chol(M)
L = LL - diag(diag(LL))
D = diag(LL)
s[i] = w[i]
ss = ss + s[i]
U = U + s[i]*L
E = E + s[i]*log(D)
}
SS = U/ss + diag(exp(E/ss))
M_JMA_w_in[,,j] = t(SS)%*%SS
}
} else{
for (j in 1:MC) {
ss = 0
U = 0
s = array(0,m_w_in)
for (i in 1:m_w_in) {
M = (Mest[,,i,j] +t(Mest[,,i,j]))/2
L = chol(M)
s[i] = w[i]
ss = ss + s[i]
U = U + s[i]*L
}
M_JMA_w_in[,,j] = t(U/ss) %*% (U/ss)
}
}
if(corrOut){
for(j in 1:MC){
D=diag(1/sqrt(diag(M_JMA_w[,,j])))
M_JMA_w_in[,,j]=D%*%M_JMA_w_in[,,j]%*%D
M_JMA_w_in[,,j]=as.matrix(Matrix::forceSymmetric(M_JMA_w_in[,,j]))
}
}
JMA_w = list(M_JMA_w=M_JMA_w_in, optns=list(corrOut=corrOut,metric=metric))
return(JMA_w)
}
最後にモデル平均推定量の$\text{FPE}$を計算する関数を準備します．
library(frechet)
library(progress)
#########################
# Generate training data
#########################
n=50 #サンプルサイズ
MC=100 #FPEの計算用にデータと独立にモデルから新たに発生させるデータの数
p=9 #予測子の数
m=5 #行列のサイズ
r=0.5 #予測子の相関
n_model=3 #組み合わせるモデルの数
n_iter = 100 #繰り返し回数
####################################################
# Compute JMA estimator
####################################################
FPE_weight_save1 = array(dim=c(n_iter,(1+n_model)))
w_JMA1 = weight_generate(n_model,0.01)
w_model1 = nrow(w_JMA1)
m_w1 = ncol(w_JMA1)
JMA_w1 = array(0,w_model1)
pb1 <- progress_bar$new(total = n_iter,
format = "[:bar] :percent complete",
clear = TRUE)
for (iter in 1:n_iter){
pb1$tick()
#############################
# Generate training/test data
#############################
train_data=gendata(n,p,m,r)
Mtrain=train_data$Min
test_data=gendata(MC,p,m,r)
Mtest=test_data$Min
#############################################
## Compute Global Frechet regression function
#############################################
#Compute prediction values
Cov_est = array(dim=c(n_model,m,m,MC))
for (im in 1:n_model){
GlCov_est = GloCovReg(x=train_data$x[,1:im],M=Mtrain,xout=test_data$x[,1:im],optns=list(corrOut=FALSE,metric="cholesky"))
for(imm in 1:MC){
Cov_est[im,,,imm]=GlCov_est$Mout[[imm]]
}
}
##########################################
## Compute weight and FPE of MA estimator
##########################################
for (jt in 1:w_model1){
m_JMA_w_CV1 = array(0,n)
for (it in 1:n){
##Compute Frechet regression function for nested models
Cov_est1_CV1=GloCovReg(x=train_data$x[-it,1],M=train_data$Min[,,-it],xout=train_data$x[it,1],optns=list(corrOut=FALSE,metric="cholesky"))
Cov_est2_CV1=GloCovReg(x=train_data$x[-it,1:2],M=train_data$Min[,,-it],xout=matrix(train_data$x[it,1:2],ncol=2),optns=list(corrOut=FALSE,metric="cholesky"))
Cov_est3_CV1=GloCovReg(x=train_data$x[-it,1:3],M=train_data$Min[,,-it],xout=matrix(train_data$x[it,1:3],ncol=3),optns=list(corrOut=FALSE,metric="cholesky"))
Mest_CV1 = array(dim=c(m,m,m_w1,1))
Mest_CV1[,,1,1] = Cov_est1_CV1$Mout[[1]]
Mest_CV1[,,2,1] = Cov_est2_CV1$Mout[[1]]
Mest_CV1[,,3,1] = Cov_est3_CV1$Mout[[1]]
## Compute JMA estimators for a given weight
m_JMA_w1 = GFRCovCholesky_w(1,Mest_CV1,w_JMA1[jt,],optns = list(metric="cholesky"))
JMA_CV1 = m_JMA_w1$M_JMA_w
A = train_data$Min[,,it]
B = JMA_CV1[,,1]
CV_dist = (dist4cov(A=A,B=B,optns=list(metric="cholesky"))$dist)^2
m_JMA_w_CV1[it] = mean.default(CV_dist)
}
## Compute the value of CV criterion for a given weight
JMA_w1[jt] = mean.default(m_JMA_w_CV1)
}
## Compute estimated weight
w_est1 = w_JMA1[which(JMA_w1==min(JMA_w1)),]
## Compute MA estimator for the estimated weight
Mest = array(dim=c(m,m,n_model,MC))
for (im in 1:n_model){
for (i in 1:MC){
Mest[,,im,i] = Cov_est[im,,,i]
}
}
m_JMA_w1 = GFRCovCholesky_w(MC,Mest,w_est1,optns = list(metric="cholesky"))
###################################
## Compute FPE of the MA estimator
###################################
FPE_dist_JMA = array(dim=c(1,MC))
for (i in 1:MC){
A = Mtest[,,i]
B1 = m_JMA_w1$M_JMA_w[,,i]
FPE_dist_JMA[1,i] = (dist4cov(A=A,B=B1,optns=list(metric="cholesky"))$dist)^2
}
FPE_JMA1 = mean.default(FPE_dist_JMA[1,])
FPE_weight_save1[iter,] = c(FPE_JMA1, w_est1)
Sys.sleep(1 / n_iter)
}
FPE_JMA_mean1 = mean.default(FPE_weight_save1[1:n_iter,1])
JMA_weight_mean1 = apply(FPE_weight_save1[1:n_iter,2:4],2, mean.default)
FPE_JMA_mean1
JMA_weight_mean1
まとめ
この記事では大域的フレシェ回帰のモデル平均に関する研究成果 (Kurisu and Otsu (2025)) を紹介しました．株式会社Nospareには今回の記事で紹介したランダムオブジェクト解析に限らず，統計学の様々な分野を専門とする研究者が所属しています．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．
参考文献
[1] Kurisu, D. and Otsu, T. (2025) Model averaging for global Frechet regression. R&R at Journal of Multivaraite Analysis.