ぱらぱらめくる『Functional Data Analysis with R and MATLAB』

Functional Data Analysis with R and MATLAB (Use R!)

作者: James Ramsay,Giles Hooker,Spencer Graves
出版社/メーカー: Springer New York
発売日: 2009/07/02
メディア: ペーパーバック
購入: 1人クリック: 1回
この商品を含むブログを見る

MATLABで作ってRに変換し、R用の本にした、という経緯なのではないかと思われる本。Rに特化して書いてほしかった…と思う点が多い
色々な処理の薀蓄より、「データ」はこれ、「処理」はこれ、「結果」はこれ、「その背景は…」と書いていてくれるとありがたいのだが…細かいことは、パッケージの関数説明を見ればわかるわけで、fdaパッケージの関数ヘルプを見てもわからないのは、「こういうデータにこうしたいとき、どの関数をどう使う？」という点なわけなので
と、少し、内容に不満はあるけれど、それでも、まあ、『ぱらぱらめくる』だけだから我慢、と(そこを自分でメモしておくと勉強になる、とポジティブに考えるのもよいかもしれません)
目次
- 1. 関数解析への導入
- 2. RとMATLABの必須事項
- 3. 関数を構成するための、Basis Systemの作り
- 4. 関数解析用オブジェクトの作り方
- 5. 平滑化：ノイズありデータからの曲線計算
- 6. 関数解析データの記述統計
- 7. 関数解析データのばらつきのとらえ方：関数解析のためのPrincipal component analysis, Canonical component analysis
- 8. レジストレーション：カーブデータのアラインメント
- 9. スカラー観測データの関数解析的線形モデル
- 10. 複数次元データの関数解析的線形モデル
- 11. 関数モデルと力学系
1. 関数解析への導入
- 時系列データは１階微分・２階微分、と微分してこそ見える特徴がある
- 時系列データは対数化することで見える特徴がある
- 時系列データは、時間軸をそろえることで見える特徴がある
- 時系列デーはは、複数のものを使って、時間をパラメタとしたデータ空間曲線とすることで見える特徴がある
- 関数解析の初歩
  - Data representation: 平滑化と補間
  - Data registration: アラインメント
  - Grouping: グループにまとめる・クラスタリング
  - Phase-Plane plots: 微分ｘ微分のプロットで特徴を可視化
- 関数解析データのばらつきの評価
  - 記述統計
  - PCA
  - Canonical Correlation
  - 線形モデル
- 関数解析における微分の利用
2. MATLAB と R
- 相互の異同
3. 関数を構成するための、Basis Systemの作り
- 関数のセット"basis functions"を作って、その線形和のための係数を納めるオブジェクトを作って、それらによって関数を取り扱う
- basis functionsに1,2,..階の微分も登録して関数解析に備える
- Basis systemsのいろいろ
  - フーリエ(周期データ)
  - スプライン(非周期な滑らかなデータ)
  - Constant, Monomial and Other Bases

# install.packages("fda")
library(fda)
rangeval <- c(-1,2)
basisobj <- create.constant.basis(rangeval)
plot(basisobj)
basisobj2 <- create.monomial.basis(rangeval,5)
plot(basisobj2)
basisobj3 <- create.fourier.basis(rangeval,7)
plot(basisobj3)

4. 関数解析用オブジェクトの作り方
- Basis 関数の線形和としてあらわすためのオブジェクト
  - たとえば、サインカーブに乗っているデータ点を、bsplineの関数たちで表させると、結構似た形にできることが以下でわかる

n.pt <- 15 # データ点の数を基底関数の数としてやる(それ以上増やしても『情報がない』から仕方ない)
basis13 <- create.bspline.basis(c(0,1),n.pt)
tvec <- seq(from=0,to=1,length=n.pt)
sinecoef <- sin(2*pi*tvec)
sinefd <- fd(sinecoef,basis13,list("t","","f(t)"))
plot(sinefd)
points(tvec*1,sinecoef)

- Linear differential operator (Lfd class)というのがあって、指定階までの微分シリーズをハンドリングできる
5. 平滑化：ノイズありデータからの曲線計算
- 平滑化はRにもやり方があるけれど、fdaパッケージでは、簡単に一つの関数smooth.basis()でできるようにしてある

n <- 51
t <- sort(runif(n))
x <- sin(4*pi*t)
s <- 0.2 
y = x + rnorm(x)*s
nbasis = 13
basisobj = create.bspline.basis(c(0,1),nbasis)
plot(basisobj)
# 出力に簡易版といろいろな情報が返ってくる版がる
ys = smooth.basis(argvals=t, y=y, fdParobj=basisobj)
Ys = smooth.basis(argvals=t, y=y, fdParobj=basisobj,returnMatrix=TRUE)
xfd = ys$fd
Xfd = Ys$fd
# データと平滑結果をプロットする
plotfit.fd(y, t, xfd)

- Basisの構成関数の数を増やすとフィットはよくなるが、それを自由度に照らして、よいあたりを得るためにgcv(Generalized Cross Validation)なる統計量が返ってくる。これを最小化するようなbasisのサイズを探してみる
  - 最小二乗(の平方根)よりも、gcvの方がよいらしい

ns <- 4:30
gcvs <- rmses <- rep(0,length(ns))
for (i in 1:length(ns)) {
	nbasis <- ns[i]
  basisobj = create.bspline.basis(c(0,1),nbasis)
  ys = smooth.basis(t, y, basisobj)
  xfd = ys$fd
  gcv = ys$gcv
  RMSE = sqrt(mean((eval.fd(t, xfd) - x)^2))
  gcvs[i] <- gcv
  rmses[i] <- RMSE
}
matplot(ns,cbind(gcvs,rmses),type="l")

6. 関数解析データの記述統計
- フィッティングなり平滑化なりをして、それをどう「見せるか」の話
7. 関数解析データのばらつきのとらえ方：関数解析のためのPrincipal component analysis, Canonical component analysis
- pca.fd()とその引数を作るfd.Par()
- 同様にcca.fd()にもfd.Par()の出力を使う
8. レジストレーション：カーブデータのアラインメント
- landmarkreg()とregister.fd()
9. スカラー観測データの関数解析的線形モデル 10. 複数次元データの関数解析的線形モデル
- 曲線らしいデータがいくつもあったときに、それに「同様の曲線表現」を想定しつつ、係数を変えてフィッティングする、しかも、それを曲線の属性との関連を想定してフィッティングする
- ここでは、線形モデルだから、曲線の属性を説明変数として、係数を線形モデルで決める、ということ(のようだ)
- fRegress()で実施
- 線形回帰部分はlm()を使う
- 検定をするが、それはパーミュテーションテスト。Fperm.df()
時間の関数であったり、と何かの関数であることを前提とするが、そうでなければ、ただの線形回帰。それを示す

# 適当なデータ
     ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
     trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
     group <- gl(2,10,20, labels=c("Ctl","Trt"))
     weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)
fRegress.D9 <- fRegress(weight ~ group)
(lm.D9.coef <- coef(lm.D9))
(fRegress.D9.coef <- sapply(fRegress.D9$betaestlist, coef))
all.equal(as.numeric(lm.D9.coef), as.numeric(fRegress.D9.coef))

.D9.coef <- coef(lm.D9))
(Intercept)    groupTrt 
      5.032      -0.371 
> 
> (fRegress.D9.coef <- sapply(fRegress.D9$betaestlist, coef))
    const group.Trt 
    5.032    -0.371

- こちらの記事のRコードの丸写し:こうやれば、「動く」ようなので、これが何をどういう手順でやっているかがわかればよさそう

# 35箇所のうち、1番目を予測対象として、34箇所のデータを使う。その1 vs. 34の分離
# pick out data and 'new data'
dailydat <- daily$precav[,2:35]
dailytemp <- daily$tempav[,2:35]
dailydatNew <- daily$precav[,1]
dailytempNew <- daily$tempav[,1]
# 年間の日平均気温
par(mfcol=c(1,2))
matplot(dailytemp,type="l")
# 年間の日平均降水量
matplot(dailydat,type="l")
par(mfcol=c(1,1))

- 周期性データなので、フーリエ変換する。粗さは適当に65basis 関数で表される程度にする
  - そうして得られた平滑カーブをプロットしておく

# set up response variable
annualprec <- log10(apply(dailydat,2,sum))

# create basis objects for and smooth covariate functions
tempbasis <- create.fourier.basis(c(0,365),65) # ある程度の細かさでフーリエ
tempSmooth <- smooth.basis(day.5,dailytemp,tempbasis)
tempfd <- tempSmooth$fd
plot(tempfd)

# create design matrix object
templist <- vector("list",2)
templist[[1]] <- rep(1,34)
templist[[2]] <- tempfd

# create constant basis (for intercept) and
# fourier basis objects for remaining betas
conbasis <- create.constant.basis(c(0,365))
betabasis <- create.fourier.basis(c(0,365),35)
betalist <- vector("list",2)
betalist[[1]] <- conbasis
betalist[[2]] <- betabasis

# set roughness penalty for betas 
# ハーモニック・アクセレレーションは、関数の微分の重みづけルールの一つでうまく周期性データの平滑性を達成するもの
# 0-365の範囲で、(0,k,0)という重みで平滑化する、という(年間のdailyデータにおける『一種のデフォルト』)
Lcoef <- c(0,(2*pi/365)^2,0)
harmaccelLfd <- vec2Lfd(Lcoef, c(0,365))
# 微分による平滑の情報ができたので、それを使って、
# ペナルティの重み変数lambda ここでは、とても大きい値を使っているが、Rのfdaパッケージのヘルプ記事のExamplesでは1を使っているのだが…これはどう決める？
# とは言え、fdPar()によって、作られている関数たちの元になっているのはbetabasisというオブジェクトで、これは、create.fourier.basis(c(0,365),35)の出力なわけで、betafdParというオブジェクトは、フーリエ使って、微分に関する重みづけルールとしてharmonic accelerationして、適当にlambdaでいじった関数群、という、『ある意味のデフォルト』ということだろう
lambda <- 10^12.5 
betafdPar <- fdPar(betabasis, harmaccelLfd, lambda)
betalist[[2]] <- betafdPar

# regress
# 最後は、そんな定数項用の関数(定数関数)とフーリエの重みに工夫した関数たちとでできたリストbetalistを引数に、フーリエ平滑化した、年間気温データtemplistを説明変数ならぬ説明関数として、annualprecというスカラー値を回帰していますよ、というのが以下の１行
annPrecTemp <- fRegress(annualprec, templist, betalist)

- ここまでで、データセットについて回帰する話は終わり
  - まとめると、従属変数データと、説明変数データ(これは平滑化とかをしておく)のリストと、係数に関するリストとを使って、関数回帰する
  - ただし、説明変数データリストと係数に関するリストとは、相互に対応づけてある
  - 説明変数データリストの方は、「切片項」と「関数項」とになっていたりする
  - 説明変数データリストの「関数項」の方は、平滑化などで作られたfdオブジェクトにしておく
  - 係数に関するリストのうち、「関数項」に対応する方は、平滑化で使ったのと同じ関数基底を使ってある。それに重みづけルール情報がある
  - ここまでわかれば、だいたいわかる
  - 後はfRegress()の出力がどうなっているのかがわかればよい
  - 幸い、このfRegress()の出力を使って、新データセットからの予測について、以下のように、「既存関数」を使わずにマニュアル実行しているので、それを見て、「fRegress()出力」の様子を見てみよう
- ここから先が、新データセットについて、予測する話

# create basis objects for and smooth covariate functions for new data
# 訓練用データ(34地点分)と同じ条件で平滑化
tempSmoothNew <- smooth.basis(day.5,dailytempNew,tempbasis)
tempfdNew <- tempSmoothNew$fd

# create design matrix object for new data
# 平滑化を表している部分と、「切片部分」とのリストを作って、予想用データを、訓練用データで行ったプロセスと併せる
templistNew <- vector("list",2)
templistNew[[1]] <- rep(1,1)
templistNew[[2]] <- tempfdNew

# convert the intercept into an fd object
# onebasisに相当する部分を、訓練用データのときはただのベクトルで扱っていたのだが、(おそらく、予想処理での都合上)、「ただの１」であっても、それなりに関数オブジェクト化しないといけないらしく、そのための処理
onebasis <- create.constant.basis(c(0,365))
templistNew[[1]] <- fd(matrix(templistNew[[1]],1,1), onebasis)

# set up yhat matrix (in our case it's 1x1)
# 予測対象は、具体的な値を持った行列。今の例では、回帰が34地点の年間降水量であって、ある１地点のそれを日別気温から予測するだけなので、予測する値は１つだけ(１ｘ１行列)
yhatmat <- matrix(0,1,1)

# loop through covariates
# この解析では、切片と級数関数セットとがリストに分けて入っているので、それぞれについて、予測対象への寄与分を計算して合算することになる。そのループ
# ここで「ループを回して合算」できるというところが「線形回帰」
p <- length(templistNew)
for(j in 1:p){
# やりやすいように、ループにおける「今」の関数セットのオブジェクトを取り出す
    xfdj       <- templistNew[[j]]
# その関数セットオブジェクトはいろいろな情報をリストという構造化したものになっているが
# そのままだとわかりにくいので、それぞれをオブジェクトとして独立させる
# 基底関数が何か
    xbasis     <- xfdj$basis
# 基底関数の数
    xnbasis    <- xbasis$nbasis
# 台(サポート)範囲(下限値と上限値)
    xrng       <- xbasis$rangeval
# 台をいくつに均等分割するかのその数。基底関数の数の10倍程度をめどにするが、それが少なすぎるときには、501は確保するようにする(501とは上下限を500セグメントに分割する、ということ)
    nfine      <- max(501,10*xnbasis+1)
    tfine      <- seq(xrng[1], xrng[2], len=nfine)
# 分割幅
    deltat     <- tfine[2]-tfine[1]
# ここまでは新データの平滑化結果templistNewについて、その出力xfdiから必要な情報を取り出して使いやすくする過程
# ついで、評価
# 予測用データの説明変数情報を平滑化して、それを関数セットオブジェクトにしてあるxfdiに、「評価するべき台上の点」を与えて、それを行列として返す
# この例の場合は、平滑化した気温のカーブになる
    xmat       <- eval.fd(tfine, xfdj)
# plot(xmat)
# 観察データの回帰の結果オブジェクトannPrecTempを使うのはこの後
    betafdParj <- annPrecTemp$betaestlist[[j]]
    betafdj    <- betafdParj$fd
# betamatは気温カーブ上の値がいくつだったら、それを従属変数にどういう寄与分にするかの係数
# したがってサイズはxmatと同じ
    betamat    <- eval.fd(tfine, betafdj)
    # estimate int(x*beta) via trapezoid rule
# 台全体について、評価点とその値があるので、「台形法」で関数カーブの下面積を計算する
# crossprod(xmat,betamat)は台に沿って「気温」の値と「日ごとに、気温ｘその日の係数」をかけて、すべて足し合わせるという作業が
# ベクトルの内積になるわけだけれど
# それを従属変数の数だけ一息で計算する、と言う工程
# 今の例では、従属変数が１スカラーなので、単純に２ベクトルの内積になっている
# それにdeltatを掛けているのは、「ベクトルの個々の要素」が寄与する部分は、小刻みにした幅に比例するから、という意味
# -0.5 * 第１行と最終行との和、というのは、右端と左端とを「台形」で考えると、足し過ぎているから、引きましょう、ということ
    fitj       <- deltat*(crossprod(xmat,betamat) - 
                      0.5*(outer(xmat[1,],betamat[1,]) +
              outer(xmat[nfine,],betamat[nfine,])))
    yhatmat    <- yhatmat + fitj
}

11. 関数モデルと力学系
- pda.fd()関数がPrincipal Differential Analysisの略