[とらにもつ]
とら着替えがいっぱい入ってる．
glm(..., family = "gaussian")
の「AIC をどう評価しているのか?」調査のつづき
……
writeLines(as.character(body((glm.fit))), con = "glm.fit.R")
で glm.fit
の内部を
glm.fit.R
出力してみると，
aic.model <- aic(y, n, mu, weights, dev) + 2 * rankという行がある． この
aic()
なる関数は「上」のほうで，
aic <- family$aicと定義されているので，
family
を調べないと行けない．
そこで，
今回はなぜかしら
family = gaussian
の場合について調べないといけないので，
writeLines(as.character(body(gaussian)), con = "gaussian.R")というふうにファイル gaussian.R に出力してみると，
aic = function(y, n, mu, wt, dev) {
nobs <- length(y)
nobs * (log(dev/nobs * 2 * pi) + 1) + 2 - sum(log(wt))
}
という定義があったので，
このように AIC が評価されているのがわかる．
もとの glm.fit()
ないで aic()
を呼ぶときの residual deviance はやはり family
中の dev.resids() 関数を使っていて，
これは gaussian の場合，
dev.resids = function(y, mu, wt) wt * ((y - mu)^2)このように定義されている．
たしかに「統計数理」の普通の号に掲載されているような数理統計学的な記事はとうてい準備できませんが，ときどき「普通ではない号」もまじることもあるので，これなら生態学関係者が執筆してもケチをつけられるいわれはないだろうと考えた次第です．そもそもこちらがお願いしたわけではなく，「統計数理」の編集担当のかたが依頼されてきたので，数理統計になっていない点について当方が気にする必要はないだろう，と．
いっぽうで，「統計数理」の記事は PDF 化されネット上で誰でも読めるようになっているので (おそらく統数研が存続するかぎり維持されるでしょう)，生態学研究者にとっては日本語で読める「統計モデリングのあれこれ例題のまとめ」みたいなものが存在することは便利だろうと考え，こういう統数研のシステムを「こっちが利用してやろう」ぐらいの気分で取り組んでおります．
source("fig_b.R")
b1 <- function(i) bm[,sprintf("b[%i,1]", i)]
plot(density(b1(1) - b1(3))) # log(NV) - log(EW)
plot(density(b1(2) + b1(3) - 2 * b1(1))) # log(TA) + log(EW) - 2 * log(NV)
options(pager="/usr/bin/less -X")
だと help() なんかでエラーとなる．
姑息な回避わざとしては，
~/.bashrc とかで
export LESS=-XF と設定して
options(pager="/usr/bin/less")
とする．
library(rgl)
調査．
library(rgl)
を使った三次元作図プログラミングは初めてなのだが
……
いろいろ工夫してみる．
library(rgl)
の shade3d と cylinder3d()
をこんなふうに組み合わせて作図している．
shade3d( cylinder3d(center = xzy, radius = radius), color = "#a04040" # alpha = 0.5 などと透過 (semi-transparent) 指定可能 )円筒といっても下の図 (上下転倒した視点から見ている図) でわかるように， ホントの円筒ではなくポリゴンで近似していることがわかる． 円周を分割するポリゴン数は
cylinder3d()
のオプション sides で指定できる．
下の図の場合は， デフォルトである sides = 8．
library(rgl)
によるヒノキデータ作図問題のつづき．
昨日は，
かくのごとくとりあえずできてめでたしめでたし
……
とかよろこんでいたわけだが
……
cylinder3d()
の debug = TRUE を指定して作図させてみる
……
library(rgl) 作図は途中で Ctrl-C で停止しても問題ない
……
で，
わかったのは，
ねじれているところで局所座標系が「反転」している!
ということ
……
図中の緑・赤の線分 (たぶん青もあるけど幹の中に隠れている)
は局所座標系をあらわしている．
cylinder3d()
の引数 e2
で法線ヴェクター (normal vector) を指定できるので，
そこに e2 = cbind(1, 0, rep(0,nrow(xzy)))
てなかんじで無理やり指定すればよい，
というかんじ
……
このあたりの曲線指定みたいなものは
フレネ・セレの公式
(Frenet–Serret formulas)
も参照．
(cylinder3d()
の 関数定義コード)
