実際のデータからRの二変量正規分布のパラメータを推定するにはどうすればよいですか？

Question

2つの正規分布XとYからなる二変量正規分布でモデル化したい実際のデータからxyペアのセットを取得しました。使用しなくても再作成できるようにパラメータを計算したいそれはあまりにも高価です（百万行）として元のソースデータ。ディストリビューションXのX

標準偏差分布の

• 平均：私は、次のパラメータを推定する必要があると思う

  hexbinplot(x~y, data=xyPairs, xbins=16) 
  

  ：私は成功したと、このデータをプロットしています現時点では

• 平均分布Y
• 配布の標準偏差シグマ行列

を作成するために使用されるイオンY

のRhoは、その後変量正規がで指定されています

Rでこれを行うには、パッケージはありますか？

多くのパッケージを調べましたが、実際のデータをモデル化する2変量正規分布を作成するのではなく、ランダムデータを使用して2変量をシミュレートするのに役立ちます。

詳細をご希望の場合はお知らせください。

Answer 1

オクラホマので、のは、いくつかの事実から始めましょう：

1. あなたが多変量正規分布を持っている場合は、周辺分布が出て取り残されている変数に関連するすべてのパラメータに依存しません。 here
2. パラメータmuおよびsigma^2の最尤推定量は、サンプルの類似体に対応することがよく知られています。一変量の場合に解析解を得る方法の例については、hereを参照してください。

これらのパラメータを次のように推定できます。まず、私はいくつかのサンプルデータを生成してみましょう：

n <- 10000 
set.seed(123) #for reproducible results 
dat <- MASS::mvrnorm(n=n, 
       mu=c(5, 10), 
       Sigma= matrix(c(1,0.5,0.5,2), byrow=T, ncol=2) 
) 

をここで、私は5と10、それぞれであることをmu1とmu2を選択しました。また、sigma1^2は、1 rho*sigma1*sigma2等しい0.5に等しく、sigma2^2はrho * sigma1 * sigma2 = 0.5以来、私たちはrho = 0.5/sqrt(1*2) = 0.35

は、現在知られている（分析）最大尤度を

推定使用することを持っていることを、私たちはパラメータmu1を推定してみましょう。2.ノートに等しいですデータから最初にmu2を入力します。ここでは、事実1は私が依存関係について心配する必要がないことを保証するので、それぞれの変数のサンプル手段を使用します。つまり、限界分布には同一のパラメータがあるので、それらは2値的には正常であるとは無視できます。単変量の場合のこれらのパラメータのMLEが標本平均であることがわかりました。

> colMeans(dat) 
[1] 5.006143 9.993642 

これは、以前にデータを生成するときに指定した真の値に非常に近いことがわかります。

今、私たちはx1とx2の分散を推定してみましょう：

> apply(dat, 2, var) 
[1] 0.9956085 2.0008649 

はまた、これが真の値にかなり近づきます。このアプローチは今のところうまくいくと思われます。 :)

ここで残っているのはすべてrhoです。分散共分散行列の対角外の項目はrho*sigma1*sigma2 = rho * 1 * sqrt(2)です。これは0.5と定義されています。従って、rho = 0.35。

ここで、サンプルの相関関係を見てみましょう。サンプル相関はすでに共分散を標準化しているため、相関係数を得るために手動でsqrt(2)で除算する必要はありません。

> cor(dat) 
      [,1]  [,2] 
[1,] 1.0000000 0.3481344 
[2,] 0.3481344 1.0000000 

これは、以前に指定された真のパラメータにかなり近くなっています。後者は小さなサンプルに偏っていると主張し、修正を加えることができることに注意してください。議論については、Wikipediaの記事を参照してください。それをやりたければ、最後の言葉にn/(n-1)を掛けるだけです。 n=10000などのサンプルサイズでは、大きな違いはありません。

ここで何をしましたか？私は、これらの量の分析的最尤推定量がどのように見えるかを知っていたので、これらのパラメータを推定するために使用しました。分析的にどのようにソリューションが見えるのかわからない場合はどうしますか？原則として、あなたは尤度関数を知っています。あなたはそのデータを持っています。パラメータの関数として尤度関数を書き留めて、多くの利用可能なオプティマイザの1つを使用して、サンプル尤度を最大にするパラメータの値を見つけることができます。これはダイレクトMLアプローチになります。 hereを参照してください。

だから、試してみましょう。尤度数値

上記の手順を最大化

は、我々が解析的最尤推定量を得ることができたという事実を使用します。すなわち、尤度関数の導関数を取り、それをゼロに等しく設定し、未知量について解くことによって、これらの量に対する閉形式解を見いだした。しかし、コンピュータを使って数値を数値で見つけることもできます。これは、扱いやすい解析ソリューションが見つからない場合に便利です。それを試してみましょう。

まず、関数を最大化するので、組み込み関数optimを使ってみましょう。 optimでは、パラメータのベクトルに最初の開始値を指定し、パラメータベクトルを引数として渡す関数を指定する必要があります。この関数は、最大化または最小化される値を返すものとします。

この機能はサンプルの可能性になります。サイズnのiid-サンプルが与えられると、サンプル尤度は全ての尤度（すなわち、確率密度関数）の積である。大規模な製品の数値最適化は可能ですが、一般的には、製品を合計に換算するために対数をとっています。可能性を取得するには、ちょうど変量正規分布の個々のPDFファイルで、長い期間懸命に見て凝視し、そしてあなたは、サンプルの可能性は、この関数は、引数を介して最大化する

-n*(log(sig1) + log(sig2) + 0.5*log(1-rho^2)) - 
0.5/(1-rho^2)*(sum((x1-mu1)^2)/sig1^2 + 
       sum((x2-mu2)^2)/sig2^2 - 
       2*rho*sum((x1-mu1)*(x2-mu2))/(sig1*sig2)) 

のように書くことができていることがわかります。 optimは1つのパラメータベクトルを供給するために私を必要とし、私はこのためにラッパーを使用して、次のように最大化問題を設定して以来：

​​

次のようにoptimへの私の呼び出しは、その後になります

ここ

eps <- eps <- .Machine$double.eps # get a small value for bounding the paramter space to avoid things such as log(0). 

numML <- optim(rep(0.5,5), wrap, dat=dat, 
     method="L-BFGS-B", 
     lower = c(-Inf, -Inf, eps, eps, -1+eps), 
     upper = c(Inf, Inf, 100, 100, 1-eps), 
     control = list(fnscale=-1)) 

を、 rep(0.5,5)は開始値を提供し、wrapは関数の上にあり、lowerとupperはパラメータの境界であり、引数はで、関数が最大化されていることを確認します。結果として、私は得る：

numML$par 
[1] 5.0061398 9.9936433 0.9977539 1.4144453 0.3481296 

注これらの要素がmu1、mu2、sig1、sig2とrhoに対応します。 sig1とsig2を四角で囲むと、私が最初に提供した分散が再現されることがわかります。だから、うまくいくようです。 :)