由买买提看人间百态

topics

全部话题 - 话题: rmvnorm
(共0页)
x***2
发帖数: 946
1
来自主题: Statistics版 - 请教一个R的问题
我在function里面生成multivariate norm变量时出错。
test<-function(rep,n1,n2,mu1,mu2,x,v1,v2,vx,rho){
...
vx1y<-rho*sqrt(v1*vx)
vx2y<-rho*sqrt(v2*vx)
xy1<-rmvnorm(n=n1,mean=c(mu1,x),sigma=matrix(c(v1,vx1y,vx1y,vx),2,2))
xy2<-rmvnorm(n=n2,mean=c(mu2,x),sigma=matrix(c(v2,vx2y,vx2y,vx),2,2))
在这里就报错,错误信息
Error in rmvnorm(n = n1, mean = c(mu1, x), sigma = matrix(c(v1, vx1y, :
mean and sigma have non-conforming size
如果我改成数字,function就没问题。
# xy1<-rmvnorm(n=n1,mean=c(0,0),sigma=matrix(c(1,0... 阅读全帖
p********a
发帖数: 5352
2
☆─────────────────────────────────────☆
PurpleShell (PurpleShell) 于 (Sun Jul 27 19:45:01 2008) 提到:
mean vector 知道,variance 矩阵知道,但不是正定的,如何构造具有多元正态分布
的随机向量?
☆─────────────────────────────────────☆
lisaxy (lisa) 于 (Mon Jul 28 01:12:46 2008) 提到:
variance 矩阵知道,但不是正定的?
variance 矩阵应该都是正定的吧。
☆─────────────────────────────────────☆
howmoney (多少钱) 于 (Mon Jul 28 10:13:12 2008) 提到:
In R:
library(mvtnorm)
then use rmvnorm(n, mean=c(), sigma=matrix())

☆─────────────────────────────────────☆
z******n
发帖数: 397
3
我想我的看法不大受重视。所以构造了一个数值例子。为了使得结果能够重复,我固定
了随机数种子
set.seed(2)
library("mvtnorm")
n<-100
rho<-.9
bet<-c(.1,.1)
sigma<-matrix(c(1, rho, rho, 1), ncol=2)
x<-rmvnorm(n, sigma=sigma)
e<-rnorm(n,sd=.8)
y<-x%*%bet+e
data<-data.frame(y, x)
colnames(data)<-c("y", "x1", "x2")
mdl0<-lm(y~1, data=data)
mdl1<-lm(y~x1,data=data)
mdl2<-lm(y~x2,data=data)
mdl<-lm(y~x1+x2, data=data)
> anova(mdl0, mdl1, test="Chisq")[2, "Pr(>Chi)"]
[1] 0.03725746
> anova(mdl0, mdl2, test="Chisq")[2, "Pr(>Chi)"]
[1] 0.03311402
> anov... 阅读全帖
z******n
发帖数: 397
4
这种情况我在4楼有讨论,结论仍然是各种情况都可能发生。下面是对应于你提到的情
况的一个数值例子
set.seed(29)
library("mvtnorm")

n<-100
rho<--.9
bet<-c(.1,.1)
sigma<-matrix(c(1, rho, rho, 1), ncol=2)
x<-rmvnorm(n, sigma=sigma)
e<-rnorm(n,sd=.5)
y<-x%*%bet+e

data<-data.frame(y, x)
colnames(data)<-c("y", "x1", "x2")

pv.x1<-summary(mdl)$coefficients["x1", "Pr(>|t|)"]
pv.x2<-summary(mdl)$coefficients["x2", "Pr(>|t|)"]
pv.jnt<-anova(mdl0, mdl, test="Chisq")[2, "Pr(>Chi)"]
> c(pv.x1, pv.x2, pv.jnt)
[1] 0.03195767 ... 阅读全帖
(共0页)