p
で積分したくないパラメーターが q
だとすると，
どうしても exp(pq)
というような部分がでてくるぢゃありませんか．
せ き ぶ ん で き ん．
えーい，
どうしろってゆーわけ?!
integrate( function(c) ( dpois(growth, lambda = c * mu.nonconst) * dgamma(c, shape = shape, scale = scale) ), 0, Inf, stop.on.error = FALSE )てなかんぢで， 数値積分を含むカタチに書き換えてみる …… 計算時間は 4 倍以上を要するようになってしまった． どひゃー
dnbinom()
版の 2 倍程度の計算時間，
というところまでは回復した．
ふう．
しかしあいかわらず単純化されたモデルばかりが選ばれる．
気象の影響より「個体差」のほうが重要，
ということなのか?
うう，
それは困るんだけどなぁ．
λi
からそれぞれ独立に得られた Poisson 乱数の和は，
λi
の和をパラメーターにもつ Poisson 分布にしたがう」
というのは正しい定理なんだけど，
どうもこれを樹木の直径成長にこじつけるのは，
どこかに間違いが混入してるような気がする．
少なくとも尤度方程式は一致しない，
というあたりまえのことにいまごろ気づいた．
> x <- list(abc = 1, abd = 2, xyz = 4) > x$ab NULL > x$xy [1] 4なぜかしら「前方一致」 方式になっとる ……
glmmML()
でいーんぢゃないの，
とまた計算してみると
……
やはり呪われた start.sigma
依存性が再発してるな，
と確認できた．
glmmML()
の場合，
データが少ないときにちょっとばかり複雑なモデルを
推定させようとすると，
こういうふうに推定がむちゃくちゃになる．
たぶん
BFGS
とかいう計算方法に問題あるんだろう．
Nelder-Mead
法つかえるようにしてくれ
……
static
を使えばよろしかろう，
ということでテストコードを書いてみる．
何の問題もなく動いた．
ということで，
混合な確率論的モデルの尤度計算関数もこの計算方式に切り替えてみる．
optim()
つまりここでは Nelder-Mead 法による
パラメーター探索だけど，
C 版 (Numerical Recipe 掲載のものとか)
あるいは C++ 版 (久保自作)
に取り換えたところで使いづらくなるわりには
たいして速くならないだろう．
ということで，
速度性能の追求はこのあたりで妥協せんといかんのでは，
という気がする．
そもそもモデルだってまだまだ変更されうるわけで．
> best.model $aic 1275.83 $par scale.log const dbh.log dbh.log.sq ps.prev ct6 -1.01933 -9.40252 5.11672 -0.74374 -0.20741 0.29254 temp.selector: best=18.5/width=8.0 rain: factor=0.20/time=0.30 ppfd: factor=0.01 temp: factor=0.05…… まだまだ改善の余地がいろいろあるような．
0 < λ << 1
かつ
Prob[x > 1] = O(λ2)?
ふーむ
dvipdfm
で pdf ファイルを生成してアップロード．
で，
こういうのはだいたいにおいてアップロードしてから，
いろいろと不備な点が見つかるもんだ．
1615 ごろまで手直しに時間とられる．
sort
や分類用の名前を生成するためにこういう文字列変換が必要．
なんともはや
my $from = "がぎぐげござじずぜぞだぢづでどばびぶべぼぱぴぷぺぽ"; my $to = "かきくけこさしすせそたちつてとはひふへほはひふへほ"; my $name_sort = Jcode->new($input)->tr($from, $to)->euc;