bayerse / esreg Goto Github PK

library(esreg)
integraternorm = function(x) {qnorm(x)}
integratert = function(x) {qt(x, df = df)}

n = 2000
alpha = 0.1
p = 10
nrep = 100


i=29
set.seed(i)
gamma = 2 * rbinom(p, 1, 1/2) - 1
eta = 0.5 * rbinom(p, 1, 1/2) 
qr_norm = qnorm(alpha)
inte = integrate(integraternorm, lower = 0, upper = alpha)
es_norm = inte$value / alpha
qr = qr_norm
es = es_norm
beta_qr = gamma + eta * qr
beta_es = gamma + eta * es

err = rnorm(n, 0, 1)
X = matrix(runif(n*p, 0, 2*sqrt(3)), nrow = n, ncol = p)
Xint = cbind(1,X)
Y =  1+X%*%gamma + X%*%eta * err

# Joint regression approach
jointres = esreg(Y ~ X, alpha = alpha)
jointcoef = jointres$coefficients_e
vcovjoint = vcov(jointres)