demo/rvg-mvnorm.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/rvg-mvnorm.R
## 
## Generation of multivariate normal variates
## 2013.10.15 - final test 

library(ipeglim)

## case 1. 
## y ~ b1*x1 + b2*x2 (no intercept)
rm(list=ls())
pars <- c(-1,1)
xm <- c(3,-1)
xsd <- c(1,2)
xr <- diag(2)
xr[2] <- xr[3] <- r <- 0
xsigma <- diag(xsd)%*%xr%*%t(diag(xsd))

tmp <- rvg4yx(N=1e3, Xreg=TRUE, ztrunc=TRUE, kind="mvnorm", pars=pars, mean=xm, sigma=xsigma, center.X=FALSE)
dat <- tmp$yX
X <- tmp$X
head(X)
colMeans(X)
cov(X)
cor(X)
 
fit <- summary(ztpr(formula=y~x1+x2-1, data=dat, dist="poisson", ztrunc=TRUE))
fit

## case 2. 
## y ~ b0*x0 + b1*x1
rm(list=ls())
library(ipeglim)
pars <- c(1)
xm <- 1
xsd <- 2

tmp <- rvg4yx(N=1e2, Xreg=TRUE, ztrunc=TRUE, kind="mvnorm", cf0=-1, pars=pars, mean=xm, sd=xsd, center.X=FALSE)
dat <- tmp$yX
head(dat)
table(dat$y)
suppressWarnings(cor(tmp$X))

fit <- summary(ztpr(formula=y~x1, data=dat, dist="poisson", ztrunc=TRUE))
fit