草薙の研究ログスマートフォン用の表示で見る草薙の研究ログ英語の先生をやってます。統計GAMLSSとBAMLSS新幹線に乗ってコンビニにいくようなものだけど,GAMLSSとそのベイズ的モデルであるBAMLSSを使って,観測のデータに対して分布をフィットさせる。shape)」のこと。詳しくはwikiや,開発者らのページを。modelling大雑把にいえば,分布ベースの回帰モデルのひとつ。基本的には,(ひとつの)応答変数になんらかの分布を仮定する。正規分布とかガンマとか。ガンマ分布だったら,α(形状)とβ(比率)母数のふたつがあると。応答変数のそれぞれの値は,これらの母数によってランダムに発生していると。GAMLSSでは,これら,それぞれの母数の値自体を予測する。予測変数の値という条件における,応答変数が従うであろう分布の母数を予測して,そこから応答変数ができてるみたいな。この母数に対する予測は,線形な関係であってもいいし,加法モデルというくらいだから,スプラインとかでもいい。基本,4母数までの分布ならできるらしく(μ,σ,ν,τ),Rのgamlssパッケージには莫大な数の分布がサポートされている。変量効果も入れられる。BAMLSSというのは,これとほぼ同じようなことをベイズ的にやるんだそうだ(細かい事前分布の設定とかアルゴリズムはよくわからない)。これが面白いのは,平均値を予測するだけでなくて,散布度や形状を予測することができるようになるところ。実験条件Aは,応答変数の平均値には影響を及ぼさないけど,分散を大きくする効果をもつ,とか。歪んだ分布を仮定して,実験条件Bは,応答変数の形状に影響を及ぼす,とか。いうモデリングができるようになる。「分位点回帰のパラメトリック版を非線形でもできるようにしてる」,みたいな感じ。...それは置いておいて,このパッケージは特に便利なので,まずは単変量の分布の母数をこれらのパッケージで推定する。ある分布に従う単変量の母数の推定MASSパッケージのfitdistrやfitdistrplusパッケージでもoptim関数でも,なんなら場合によってはglm関数でも,同じことをやることができる場合は多いのだけど,gamlssパッケージやbamlssパッケージは,普通に単変量のデータの分布を推定できる。データにワイブル分布を当てはめて母数を推定する例をやってみる。これにはリンク関数としてlogが噛んでいるので,係数をexpする。すると,まあ,1と2くらい。これはただの最尤推定と同じなので,他のパッケージとかとも一致する。楽な関数があって,gamlssMLという最尤推定の関数。m2<-gamlssML(d1,"WEI")summary(m2)これもおなじことになる。切片からの回帰式とか書かなくてもよい。便利なのがあって,histDistという関数。これは,ヒストグラムとフィットした確率密度曲線を描いてくれるだけでなくて,上のサマリーも返してくれる。histDist(d1,"WEI")図はこんな感じ。オプションもあるので,色々いじれるけど,自分はやらないだろうな。BAMLSSでも基本は同じ。少なくともこの状態で,それぞれの事前分布は基本的に無情報っぽい。MCMCサンプラーはJAGSとか色々指定できるみたいだけど,デフォルトではGMCMCというのを使っているそう(これはよくわからない)。library(bamlss)#いろいろなオプションを省略して...m3<-bamlss(list(mu=d1~1,sigma=~1),family="weibull")summary(m3)それぞれの母数の95%信用区間(パーセンタイル法)も出してくれる。以下のようにしたら,MCMCサンプルも可視化できる。par(mfrow=c(1,2))plot(density(exp(unlist(m3$s[,1]))),main="Scale",xlab="Value")plot(density(exp(unlist(m3$s[,5]))),main="Shape",xlab="Value")簡単なモデルえっと,せっかくだから,ただ生データをフィットさせるだけでなくて,GAMLSSやBAMLSSの中でも,もっとも簡単な回帰モデル的なものをやってみよう。length)を応答変数,あやめの種類(3水準)を予測変数としよう。がく片の長さは正規分布に従うとして,種類は,がく片の長さの平均と標準偏差にそれぞれどのような影響をあたえるだろうか。ここであやめの種類はfactor型。データとしてはこんな感じ。virginicaという大きさのようになっている。普通の線形モデルなら,summary(lm(iris[,1]~iris[,5]))anova(lm(iris[,1]~iris[,5]))みたいにする。gamlssでは,こんな感じにする。******************************************************************普通の線形モデルと違って,標準偏差の大きさの違いも同時にモデリングしている。見た感じの結果になっている。これは,こんな感じに可視化するとわかりやすい。term.plotという関数もある。par(mfrow=c(1,2))term.plot(m4,"mu",main="Mu")term.plot(m4,"sigma",main="Sigma")もちろん,診断プロットもある。plot(m4)それぞれの予測された母数から,確率密度曲線を描くとこんな感じ。標準偏差には影響を及ぼさないとか,標準偏差には影響を及ぼすけど,平均値には及ぼさない,といったモデルを作り,適合度を比較したり尤度比検定をしたりすることもできる。これ予測変数が因子型だったら,教師ありの混合分布モデルみたいなもん。blmssパッケージでは,m5<-bamlss(list(mu=iris[,1]~iris[,5],sigma=~iris[,5]),family="gaussian")で,(当然だけど)同じような結果を得ることができる。MCMCサンプルを可視化したり,信用区間を構築してもよい。平均は変わらないけど,標準偏差だけ違う2群のモデル回帰モデルは自由に作れるので,ある因子が標準偏差だけに影響を及ぼしているモデルとかも作れる。こんなデータだとしよう。group<-as.factor(c(rep("A",100),rep("B",100)))score<-c(rnorm(100,50,10),rnorm(100,50,30))d2<-data.frame(group,class)m6<-gamlss(score~1,sigma.formula=~group,data=d2)summary(m6)term.plot(m6,"sigma",main="Sigma")結果はこんな感じ。複雑なモデルえっと,実質的に意味があるかはわからないけど(ないだろう),がく片の長さの分布を正規分布だと仮定して,がく片の幅と花びらの長さから,その予測値という条件下における分布の母数(平均,標準偏差)を推定するとしよう。その予測値と予測される母数の関係は非線形だと。#ここでは罰則付きスプラインの一種m7<-gamlss(Sepal.Length~pb(Sepal.Width)+pb(Petal.Length),sigma.formula=~pb(Sepal.Width)+pb(Petal.Length),data=iris,family="NO")summary(m6)par(mfrow=c(2,2))term.plot(m7,"mu",main="Mu")term.plot(m7,"sigma",main="Sigma")可視化するとこんな感じ。がく片の幅と花びらの長さが長くなれば,がく片の長さの分布における平均はそれぞれ単調に増加しそう(花びらの長さはちょっと指数的に見える)。がく片の幅は,がく片の長さの分布における標準偏差とはあまり関係なさそうだ。花びらの長さは,がく片の長さの分布における標準偏差に非線形の影響を及ぼしているかも(意味ないと思うけど)。これに種類を変量効果と考えて...というようにやっていくと,大分自由なモデリングができそう。ある信号の有無について多段階...2017-01-31一般化パレート分布をデータに当てはめる一般化パレート分布は所得の分布などに使われるそうだ。pnormのように...2017-01-11最尤推定した母数のもとでの確率密度曲線をヒストグラムに描き足すえっと,リクエストがあったのでここに書く。(2)はてなブログをはじめよう!kusanagikさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。