無事、この関数でできました…。ほっ。

とりあえず、ちゃんと動いていそうです。

１番と４番目の変数の間に0.7の相関を仮定しました。

平均も大体１になってる。

test.r<-rmvexp(1000,c(1,1,1,1),matrix(c(1,0,0,0.7, 0,1,0,0, 0,0,1,0, 0.7,0,0,1),ncol=4))

apply(test.r,2,mean)

pairs.panels(test.r)

f:id:senedesmus:20170217161649p:plain

ちなみに、lcmixのインストールはこちら。CRANになくて、win var.もないみたい。

install.packages("lcmix", repos="http://R-Forge.R-project.org")

多変量○○乱数発生器

擬似乱数

対数正規

一様

指数

多変量じゃないけど、参考になる

相関のある擬似乱数を”どんな分布であれ”作る方法

擬似乱数

という記事が、ここに載っていた。

統計学の基礎知識がちゃんとないので、原理はよくわからないのだけども、ほんとにそうなら便利だなぁ…。

これがあれば、相関のある二つの正規乱数から二つの相関のある指数乱数なんて簡単にできそうです。いっちょやってみたいと思います。やったらご報告します。

報告しますって言うか…。やってみたら、あっさりできた。

要するに、相関のある擬似正規乱数（例えば４つの変数で）をmvrnormか何かで作って、それをdexpするだけでした（汗）

例えば、４変数で正規乱数を１番目と４番目の変数に仮定したとき、

mvrnorm(n=1000,mu=c(1,1,1,1),Sigma=matrix(c(1,0,0,0.7, 0,1,0,0, 0,0,1,0, 0.7,0,0,1),ncol=4))

f:id:senedesmus:20170217144958p:plain

で、これをdexp()すると…？

f:id:senedesmus:20170217145551p:plain

あれ…？やっぱ、そうは問屋がおろさないのか。

しかし、正規分布している変数（もしくは一様分布とか）を指数分布に変換する関数はありそうな気がします、よね。しかしそれで複数の変数の相関まで維持されるかどうかは疑問。むーん！難しい。もうちょっと統計学の基礎知識があったらいかったのになぁ。

多項回帰

多項の割合の回帰（あるサンプルにa,b,c…がどれくらいの割合で含まれているか？を別の値から推定する）をすることになった。

北大の久保先生は「出来れば割合の回帰はDirichlet分布とかを使わないで（つまり割り算をやめて）、階層モデルでやったほうがいい」とおっしゃっていて、http://hosho.ees.hokudai.ac.jp/~kubo/ce/EcoSj2013.html

でも、割り算を避けられないケースとして

・計器の出力が割り算（何かの相対値）になっているケース

などをあげておられるのだけれども、私のデータはまさにそれにあたるので、やらないといけない（涙）

古くから「割合」の推定の研究はいろいろあるみたいなんだけども、とにかく簡便にやるならこれじゃないか？と思い、

http://dirichletreg.r-forge.r-project.org/

勉強しています。

ステップ１）DR_data()：データの成型

【使い方】

以下のような感じ

AL<-DR_data(ArcticLake[,1:3])

【（大事そうな）オプション】

Y：matrixかdata.frame。割合のデータ
trafo：論理値か実数。変数変換は、極端な0や1の値をさけるのに使える。詳しくはDetailsへ。Tureだったら”変換実施”、Falseなら”変換抑制”。0とか1とかの値が入ってるデータでsuppressすると、エラーの原因になりがち。もしこの値が実数なら、それを閾値に使う、つまり、Yのうちy<trafoやy>(1-trafo)のときに変換を実施するよ。デフォルト値はsqrt(.Machine$double.eps)、つまり、そのPCで取り扱える最小値の二乗（まぁ0と1をよけるってことか）
base：再パラメーター化モデルで使われるbaseコンポーネントデフォルト値は１
norm_tol：数値的な正確性によっては、Yの行の合計が1でないこともある。デフォルト値はsqrt(.Machine$double.eps)で、１とだいたい等しいことをテストするために使うトレランス（許容範囲）をあらわす。

【詳細】

Yの行毎の合計が１でないときには、normalizationが実施される（それぞれの行の値をその行の合計で割る。warningメッセージが出ます）ある
ある一行にNAが入っていたら、その行全部をNAにする。ベータ分布は[0,1]の区間の値を持つひとつのベクトルとして提供されるものだから。
trafo：ベータ分布が仮定されていることでメッセージが出力されるのでこの仮定にはチェックが必要。
もしtrafo=trueで実行されたら、Smithson and Verkuilen(2006)で提案された変換が実施される。Y*=[y(n-1)+1/2]/n　ここでnはYの観察値の数。（この変換はbetaregにも利用されてる。）次元の数が多ければベータ分布はディリクレ分布になるのでY*=[y(n-1)+1/d]/nに。

(つづく)