> On Jul 27, 2016, at 8:56 AM, Anamika Chaudhuri <canam...@gmail.com> wrote:
>
>
> Y<-matrix(1:40,ncol=2)
> Y1<-Y/60 # estimates of p
>
> #print(Y1)
>
> sigma2<-
> matrix(c(var(Y1[,1]),cov(Y1[,1],Y1[,2]),cov(Y1[,1],Y1[,2]),var(Y1[,2])),2,2)
> #print(sigma2)
>
> rho<-sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2])
> #rho
>
> mean(Y1[,1])
> mean(Y1[,2])
>
> #within<-matrix(data=0,nrow=20,ncol=1)
>
> for (rate3 in 1:20){
>
> rate<-Y1[i,]
> #print(rate)
>
>
> rate1<-rate/(1-rate)
>
>
> rate2<-log(rate1)
>
>
> Sigma11<-(1/(rate[1]*(1-rate[1]))^2)*sigma2[1,1]
> Sigma22<-(1/(rate[2]*(1-rate[2]))^2)*sigma2[2,2]
> Sigma12<-(1/((rate[1]*(1-rate[1]))*(rate[2]*(1-rate[2]))))*sigma2[1,2]
>
> Sigma2<-matrix(c(Sigma11,Sigma12,Sigma12,Sigma22),2,2)
> #print(Sigma2)
>
> rate3<-mvrnorm(1000, mu=c(rate2[1],rate2[2]), Sigma2)
>
> #print(rate3)
>
> x<-exp(rate3[,1])/(1+exp(rate3[,1]))
> y<-exp(rate3[,2])/(1+exp(rate3[,2]))
> print(x) # Need to be able to stack the x's to produce one matrix
> }
Not reproducible.
No loading of package for mvnorm (probably MASS) and suspect that you meant the
loop to be over 'i' rather than 'rate3' since Y1[i,] throws an error. If my
guesses are correct then perhaps:
require(MASS)
X <- matrix(NA, nrow=20, ncol=1000) # predimentsion with NA's
Y<-matrix(1:40,ncol=2)
Y1<-Y/60
sigma2<-
matrix(c(var(Y1[,1]),cov(Y1[,1],Y1[,2]),cov(Y1[,1],Y1[,2]),var(Y1[,2])),2,2)
rho<-sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2])
mean(Y1[,1])
mean(Y1[,2])
for (rate3 in 1:20){
rate<-Y1[rate3,]
rate1<-rate/(1-rate)
rate2<-log(rate1)
Sigma11<-(1/(rate[1]*(1-rate[1]))^2)*sigma2[1,1]
Sigma22<-(1/(rate[2]*(1-rate[2]))^2)*sigma2[2,2]
Sigma12<-(1/((rate[1]*(1-rate[1]))*(rate[2]*(1-rate[2]))))*sigma2[1,2]
Sigma2<-matrix(c(Sigma11,Sigma12,Sigma12,Sigma22),2,2)
rate3<-mvrnorm(1000, mu=c(rate2[1],rate2[2]), Sigma2)
x<-exp(rate3[,1])/(1+exp(rate3[,1]))
# this is not really "stacking"
# because your code was overwriting x values each time through the loop.
X[ rate3 , ] <- x # save to single row at a time
y<-exp(rate3[,2])/(1+exp(rate3[,2])) # not sure what you wanted to do with
these
# print(x) # and DO NOT PRINT 1000 item vectors to screen!
}
David Winsemius
Alameda, CA, USA
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
