non-normalな乱数を発生するパッケージあれこれ

はじめに

Rは基本のパッケージに様々な確率分布の乱数を発生させる関数が備わっているが，研究の都合で正規分布に従わない歪んだデータを発生させる必要がでてきたので，Rで歪んだデータの乱数を発生させる際に使えるパッケージをいくつか試してみた．

fungible

まずはfungibleというパッケージに入っているmonte1という関数．変数の数，サンプルサイズ，母相関行列，歪度のベクトル，尖度のベクトルなど必要な情報を入れると，指定した通りの標準化された多変量の歪んだ乱数を生成してくれる．

library(pacman)
p_load(fungible)
Sigma <- matrix(0.3,nrow=4,ncol=4)
diag(Sigma) <- 1
skew <- c(0, 0.25, -0.3, 1.25)
kurt <- c(0, 0, 1.5, 2.5)
samplesize <- 1000
nvar <- ncol(Sigma)

dat <- monte1(seed=1234,nvar=nvar, nsub=samplesize, cormat=Sigma,
              skewvec=skew, kurtvec=kurt)

psych::describe(dat$data)
psych::pairs.panels(dat$dat, smooth = F, hist.col = "white")

サンプルサイズが1000程度だと歪度も尖度もそこそこずれているがサンプルサイズを5000程度まで増やすとそれなりに正確になるようだ．

seedの指定が必ず必要なので，シミュレーションとかで使う時はseedをずらしていくようにプログラムを組む必要がある． 乱数発生の手法の元となっているのは，Fleishman(1978)とVale & Maurelli(1983)の論文．

https://link.springer.com/article/10.1007/BF02293811
https://link.springer.com/article/10.1007/BF02293687

SimMultiCorrData

次はSimMultiCorrDataというパッケージに入っているnonnormvar1という関数では1変量の歪んだ乱数を生成することができる．引数には，FleishmanかPolynomialのどちらの方法か，平均，分散，歪度，尖度のほか，サンプルサイズや，seedを指定できる．デフォルトでは平均=0で分散=1である．

p_load(SimMultiCorrData)
samplesize <- 1000
dat1 <- nonnormvar1(method="Fleishman",
                   skews=-0.6, skurts=0,n=samplesize)
dat2 <- nonnormvar1(method="Fleishman",
                    skews=0, skurts=1.25,n=samplesize)

この関数で発生させた乱数のデータは$continuous_variableに入っており，要約統計量は$summary_continuousにまとまっており，狙い通りの変数ができたかどうかの確認ができて便利である．

ちなみに引数でmethod="Polynomial"を指定すると次の論文の手順で乱数を生成するようである．

https://www.sciencedirect.com/science/article/abs/pii/S0167947302000725?via%3Dihub

このパッケージは多変量のデータもより細かく指定して作ることができる．色々な分布に従うものを作れるようである．やり方は以下に載っているのでそのうち試す．

https://cran.r-project.org/web/packages/SimCorrMix/vignettes/workflow.html

miceadds

miceaddsというパッケージのfleisman_simという関数も単変量のnon-normalなデータを生成することができる．使い方はシンプル．

p_load(miceadds)
dat1 <- fleishman_sim(samplesize, mean=0, sd=1, skew=0.6,kurt=0)
dat2 <- fleishman_sim(samplesize, mean=0, sd=1, skew=0,kurt=2)

このパッケージにはfleishman_coefという関数があってFleishmanのpower taransformationに必要な係数を得ることができる．Fleishmanによる変換では，もとの変数を$X$として

$$ Y = a + bX + cX^2 + dX^3 $$

次のような変換で，non-normalなデータを生成する．fleishman_coefでは平均，標準偏差，歪度と尖度を指定するとこの変換に必要な4つの係数$a, b, c, d$を算出してくれる．例えば，歪度=0，尖度=0である正規分布を歪度=0.6，尖度=0である変数$Y_1$と，歪度=0,尖度=1.25の変数$Y_2$に変換してみるとする．

x <- rnorm(samplesize)
coef_y1 <-fleishman_coef(mean=0,sd=1,skew=0.6, kurt=0)
# >coef_y1
# a           b           c           d
# -0.11830134  1.06186167  0.11830134 -0.02600121
coef_y2 <-fleishman_coef(mean=0,sd=1,skew=0, kurt=1.25)
#> coef_y2
# a             b             c             d
# 2.446955e-06  8.843270e-01 -2.446955e-06  3.718215e-02
y1 <- coef_y1["a"] + coef_y1["b"]*x + coef_y1["c"]*x^2 +  coef_y1["d"]*x^3
y2 <- coef_y2["a"] + coef_y2["b"]*x + coef_y2["c"]*x^2 +  coef_y2["d"]*x^3

おわりに

紹介したパッケージには色々と便利な関数がありそうなので，必要に応じて調べてここに追記していく予定である．

• ←Previous
• Next→