library(MASS) source("EEL.R") options(digits=3, warn=-1) #options(error=quote(dump.frames("testdump", TRUE))) # suppress warnings and better output set.seed(10) n=100 p=20 # true coef q=0 # IVs beta=c(3,1.5,0,0,2,rep(0,p-5)) sm1=matrix(0.5,p,p) diag(sm1)=1 #model parameter variance covarianc matrix lamv=seq(0.2,3,by=0.1) thres.v=0.05 cc=rep(0,5) m=function(theta,xy) { x=xy[,2:(p+1)] y=xy[,1] #browser() xy[,-1]*drop(y-x%*%theta) } mstar=function(theta,xy,lambda) { arg=1+m(theta,xy)%*%lambda x=xy[,2:(p+1)] wts=drop(llogp(arg,1/nrow(x))) -t(xy[,-1])%*%diag(wts)%*%x } ### Now try a linear regression x=mvrnorm(n,rep(0,p),sm1) y=x%*%beta+rnorm(n) ## q more IVs, which are in the last q columns xyz=cbind(y,x, matrix(rnorm(n*q),n,q)) theta0=drop(solve(t(x)%*%x,t(x)%*%y)) dd=bic.ee.el(m,xyz,theta0,mstar,scadlam=seq(0.2,3,by=0.05))