« 2015年01月 | メイン

2015年04月 アーカイブ

2015年04月05日

[R] ウエイト付共分散計算

cov.wt(x, wt = rep(1/nrow(x), nrow(x)), cor = FALSE, center = TRUE, method = c("unbiased", "ML"))
method = "ML"

とすると、(n-1)で割った場合でなくnで割った値を算出.

z <- rbind(c(0,1), c(0,1), c(0,1), c(0,1), c(0,1))
z2 <- rbind(z, -z)
z2
      [,1] [,2]
 [1,]    0    1
 [2,]    0    1
 [3,]    0    1
 [4,]    0    1
 [5,]    0    1
 [6,]    0   -1
 [7,]    0   -1
 [8,]    0   -1
 [9,]    0   -1
[10,]    0   -1
w2 <- c(rep(1,5), rep(1,5))
cov.wt(z2, w2, method = "ML")

# --- これと同じ共分散マトリックス
z3 <- c(0,1)
z4 <- rbind(z3, -z)
z4
   [,1] [,2]
z3    0    1
      0   -1
      0   -1
      0   -1
      0   -1
      0   -1
w4 <- c(5, rep(1,5))
cov.wt(z4, w4, method = "ML")


確率を行として入れる場合

z5 <- rbind(c(0,0), c(0,0), c(0,1), c(0,1), c(0,1))
z6 <- rbind(z5, -z)
z6
      [,1] [,2]
 [1,]    0    0
 [2,]    0    0
 [3,]    0    1
 [4,]    0    1
 [5,]    0    1
 [6,]    0   -1
 [7,]    0   -1
 [8,]    0   -1
 [9,]    0   -1
[10,]    0   -1
w6 <- c(rep(1,5), rep(1,5))
cov.wt(z6, w6, method = "ML", center = F)
$cov
     [,1] [,2]
[1,]    0  0.0
[2,]    0  0.8

$center
[1] 0

$n.obs
[1] 10

$wt
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

# --- これと同じ共分散マトリックス
z7 <- c(0,sqrt(0.6))
z8 <- rbind(z7, -z)
z8
   [,1]       [,2]
z7    0  0.7745967
      0 -1.0000000
      0 -1.0000000
      0 -1.0000000
      0 -1.0000000
      0 -1.0000000
w8 <- c(5, rep(1,5))
cov.wt(z8, w8, method = "ML", center = F)
$cov
     [,1] [,2]
[1,]    0  0.0
[2,]    0  0.8

$center
[1] 0

$n.obs
[1] 6

$wt
[1] 0.5 0.1 0.1 0.1 0.1 0.1

各行が0or1のフラグだと考え、その確率を行として加えたい場合は、上記のようにする.
分散は平均周りの二次のモーメントであるため1行目から5行目までの確率0.6になるように

sqrt(0.6)

をセットし、そのときのウエイトを5とすることで、同じ共分散マトリックスが得られる.

ただし、

center = F

の場合に限る.

なぜかというと、「center=T」では、各列の平均を計算し平均で各列を基準化するのだが、z6とz8で平均が同じにならないため.

[R] 多次元正規分布

いまさらですが、多次元正規分布.

ライブラリmvtnormで、多次元正規分布の密度関数などが提供されています.

> dmvnorm(x=rbind(c(0,0),c(0,1)), mean=c(1,1))
[1] 0.05854983 0.09653235

マトリックスをxとして入力すると、各行の密度を返してくれます.

About 2015年04月

2015年04月にブログ「[R言語] Rのお部屋::あーるのおへや[別館]」に投稿されたすべてのエントリーです。過去のものから新しいものへ順番に並んでいます。

前のアーカイブは2015年01月です。

他にも多くのエントリーがあります。メインページアーカイブページも見てください。

Powered by
Movable Type 3.37