はじめに
東京大学/株式会社Nospare リサーチャーの栗栖です．
前回の記事では経験尤度法を用いて多様体データのフレシェ平均を推定する方法について紹介しました．今回の記事では，Kurisu and Otsu(2024)に基づいて，Rを利用して実際に経験尤度法を用いた多様体データの中でも円周上のデータ(方向データ)の解析を行う方法について紹介します．
フレシェ平均の信頼領域
前回の記事では一般の多様体データに対するフレシェ平均の推定と信頼領域の構成法について紹介しましたが，この記事では特にデータが円周上に値をとる場合について解説します．
円周$\mathbb{S}=\{\nu= \nu(\theta) = (\cos \theta, \sin \theta)':\theta \in [0,2\pi)\}$は1次元の多様体であり，円周上の点$\nu$における接空間$T_\nu \mathbb{S} \subset \mathbb{R}$から$\mathbb{S}$への指数写像$\mathrm{Exp}_\nu: T_\nu \mathbb{S} \to \mathbb{S}$は以下の形で与えられます：
\mathrm{Exp}_\nu(x) = \cos(\|\tilde{x}\|)\nu + \sin(\|\tilde{x}\|){\tilde{x} \over \|\tilde{x}\|}.
ここで$\tilde{x} \in \mathbb{R}^2$は$x$を適当に変換したベクトルです．この指数写像を用いると，$\mathbb{S}$に値をとる確率変数$X$に対するフレシェ平均$\mu$は以下を満たします．
\begin{align*}
\mathbb{E}\left[\left.{d \over d x} d^2\left(\mathrm{Exp}_\mu(x), X \right)\right|_{x=0}\right] = 0 \in \mathbb{R}.
\end{align*}
上記の関係をフレシェ平均が満たすモーメント条件として経験尤度法でフレシェ平均$\mu$を推定することを考えましょう．$\{X_i\}_{i=1}^n$を円周上に値をとる独立同分布なデータとします．経験尤度が以下で定義されることに注意しておきます(詳細については前回の記事を参照してください)：
\begin{align*}
Q(\nu)&= -\sum_{i=1}^{n}\log (1 + \lambda(\nu)'g(X_i;\nu)).
\end{align*}
$q_{1-\alpha}$を自由度$1$の$\chi^2$分布の$1-\alpha$分位点とすると，上記の結果を用いて$Q(\nu)$のレベル集合 $\{\nu \in \mathbb{S}: Q(\nu) \leq q_{1-\alpha}\}$によりフレシェ平均$\nu$の$100(1-\alpha)$%信頼領域を構成することができます．
Turtleデータ
データの説明
以下ではMardia and Jupp(1999)で紹介されている turtleデータを分析してみましょう．このデータでは76匹のカメが卵から生まれた直後にどの方向へ移動したかが記録されています．特にここで扱うデータを円周上のデータ（方向データ）とみなして，このデータのフレシェ平均(平均的なカメの移動方向)とその95%信頼領域を計算してみましょう．
分析結果
下の図はturtleデータに対して計算した経験尤度(黒線)とフレシェ平均(縦の赤線)，自由度$1$の$\chi^2$分布の$0.95$分位点$q_{0.95}$(横の赤線)をプロットしたものです．縦軸は経験尤度の値，横軸は角度を表しています．
図を見ると，経験尤度の値が$q_{0.95}$より小さくなる角度の領域が2つ(フレシェ平均が含まれる領域$R_1$とそうでない領域$R_2$)ありますが，フレシェ平均が含まれない領域は，フレシェ平均の定義において評価する関数(下図)の極大点を含む領域に対応しています．フレシェ平均は下図の関数を最小化する角度に対応するので，信頼領域を構成する際には$R_2$は除いて考えます．
以上の準備の下でフレシェ平均の95%信頼領域を計算した結果が以下の図になります．特に，灰色の丸は実際のデータに対応し，オレンジの丸はフレシェ平均(平均的なカメの移動方向)，赤線は95%信頼領域に対応します．
Rを利用した分析
最後にRを用いて今回紹介したデータ分析を実行するコードの例を紹介しておきます．
まずデータの読み込みと方向データの作成を行います．
library(MASS)
library(progress)
library(circular)
turtle <- c(8,9,13,13,14,18,22,27,30,34,
38,38,40,44,45,47,48,48,48,48,
50,53,56,57,58,58,61,63,64,64,
64,65,65,68,70,73,78,78,78,83,
83,88,88,88,90,92,92,93,95,96,
98,100,103,106,113,118,138,153,153,155,
204,215,223,226,237,238,243,244,250,251,
257,268,285,319,343,350)
turtle <- turtle*pi/180
n <- length(turtle)
t_direction <- array(0,dim=c(n,2))
for(i in 1:n){
t_direction[i,1] <- cos(turtle[i])
t_direction[i,2] <- sin(turtle[i])
}
次にデータからフレシェ平均を計算します．
#Compute Frechet mean
Fre_cand <- 128
Fre <- rep(0,Fre_cand)
ph_cand <- (seq(1:Fre_cand)-1)/Fre_cand*(2*pi)
for(i in 1:Fre_cand){
m1 <- cos(ph_cand[i])
m2 <- sin(ph_cand[i])
m <- c(m1,m2)
A <- 0
for(k in 1:n){
A <- A + acos(sum(m*t_direction[k,]))^2
}
Fre[i] <- A/n
}
ind_FM <- which(Fre == min(Fre))
FM <- c(cos(ph_cand[ind_FM]),sin(ph_cand[ind_FM]))
次に経験尤度の計算に用いる関数を作成します．
######################
#g function
g_func <- function(x,ph){#x: vector, th: theta, ph: phi
m1 <- cos(ph)
m2 <- sin(ph)
m <- c(m1,m2)
#g1 <- x[1]*cos(ph) + x[2]*sin(ph)
g1 <- -x[1]*sin(ph) + x[2]*cos(ph)
A <- acos(sum(m*x))
B <- sqrt(1-(sum(m*x))^2)
C <- -2*sin(1)*(A/B)
D <- C*g1
return(D)
}
g_func2 <- function(x,ph){#x: vector, th: theta, ph: phi
m1 <- cos(ph)
m2 <- sin(ph)
m <- c(m1,m2)
g1 <- x[1]*cos(ph) + x[2]*sin(ph)
C <- -2*sin(1)
D <- C*g1
return(D)
}
#######################
ph_vec <- (seq(1:128)-1)/128*(2*pi)
ph_N <- length(ph_vec)
EL_dual_opt <- array(0,dim=c(ph_N))
#################################
g_mat <- array(0,dim=c(n,ph_N))
for(i in 1:n){
for(k in 1:ph_N){
g_mat[i,k] <- g_func(t_direction[i,],ph_vec[k])
}
}
#エラーが出ないように調整
g_mat[15,17] <- g_func2(t_direction[15,],ph_vec[17])
g_mat[45,33] <- g_func2(t_direction[45,],ph_vec[33])
g_mat[15,81] <- 10000
g_mat[45,97] <- 10000
次にフレシェ平均の95%信頼領域を計算します．
pb <- progress_bar$new(total = ph_N,
format = "[:bar] :percent 終了まで: :eta",
clear = TRUE)
############################################
for(j in 1:ph_N){
pb$tick()
#Compute dual form of EL
EL_dual <- function(lam1){
g_vec2 <- rep(1,n)
for(i in 1:n){
g_vec2[i] <- -2*log(1+lam1*g_mat[i,j])
}
return(sum(g_vec2))
}
lambda <- optim(0,EL_dual)$par
EL_dual_opt[j] <- -EL_dual(lambda)
Sys.sleep(1/ph_N)
}
########################################
#compute confidence region
sp_point <- function(ph){
s1 <- cos(ph)
s2 <- sin(ph)
s <- c(s1,s2)
return(s)
}
dim <- 1
q_EL <- qchisq(0.95,df=dim)
ind_EL <- which(EL_dual_opt<=q_EL)
ind_CR <- which(Fre[ind_EL]<q_EL)
cr_point <- length(ind_CR)
CR_point <- array(0,dim=c(cr_point,2))
for(i in 1:cr_point){
CR_point[i,1] <- cos(ph_vec[ind_EL[i]])
CR_point[i,2] <- sin(ph_vec[ind_EL[i]])
}
計算した経験尤度をプロットしてその形状とフレシェ平均を可視化してみます．
xlabels <- c(expression(paste(0,pi,"/",128)),
expression(paste(40,pi,"/",128)),
expression(paste(80,pi,"/",128)),
expression(paste(120,pi,"/",128)),
expression(paste(160,pi,"/",128)),
expression(paste(200,pi,"/",128)),
expression(paste(240,pi,"/",128)))
xvalues <- c(1:13-1)*10
xvalues2 <- c(1:7-1)*20
plot(EL_dual_opt,type="l",xlim=c(0,ph_N),ylim=c(0,50),xaxt="n",xlab=expression(theta),ylab="EL")
axis(side=1,at=xvalues,labels=F,tcl=-0.25) #２重目盛りにする
axis(side=1,at=xvalues2,labels=xlabels) #x軸のラベルを付けたい位置に目盛りを付ける
abline(h=q_EL,col="red")
abline(v=ind_FM,col="red")
最後にフレシェ平均と信頼領域をプロットして可視化してみます．
#plot data/mean direction
par(mar = c(4, 4, 4, 2))
plot(t_direction,xlim=c(-1,1),ylim=c(-1,1),pch=19,xlab="",ylab="",cex=1.5,col=rgb(0,0,0,0.3))
par(new=T)
#Frechet mean 95% CR
plot(CR_point[,1],CR_point[,2],col=rgb(0.9,0,0,0.95),#col="red"
xlim=c(-1,1),ylim=c(-1,1),lwd=4,cex=1.5,xlab="",ylab="",type="l",xaxt="n",yaxt="n")
par(new=T)
plot(FM[1],FM[2],col=rgb(1,0.5,0,0.90),#col="darkorange",
xlim=c(-1,1),ylim=c(-1,1),lwd=2,pch=19, cex=2.5,xlab="",ylab="",xaxt="n",yaxt="n")
par(new=F)
まとめ
次回の記事ではRを利用して，実際に経験尤度法を用いた多様体データ分析を行う方法について解説しました．株式会社Nospareには統計学の様々な分野を専門とする研究者が所属しています．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．
参考文献
[1] Kurisu, D. and Otsu, T. (2024) Empirical likelihood for manifolds. R&R at Journal of the Royal Statistical Society Series B.
[2] Mardia, K.V. and Jupp, P.E. (1999) Directional Statistics, Wiley.