Nothing
agconst <-
function(y,xmat,nsim=1000){
n=length(y)
if(length(xmat)==n){xmat=matrix(xmat,ncol=1)}
k=dim(xmat)[2]
if(dim(xmat)[1]!=n){print("ERROR: number of rows of xmat must be length of y")}
if(k>3){print("ERROR: xmat can have at most 3 columns")}
if(rankMatrix(xmat)[1]<k-1e-8){print("ERROR: xmat is not full rank")}
ans=new.env()
if(k==1){
x=as.numeric(xmat)
delta=makemdelta(x)
tdist=1:nsim
one=matrix(1:n*0+1,ncol=1)
for(isim in 1:nsim){
ysim=rnorm(n)
p1=coneB(ysim,t(delta),one)$yhat
p2=coneB(ysim,-t(delta),one)$yhat
r1=sum((ysim-p1)^2)
r2=sum((ysim-p2)^2)
r0=sum((ysim-mean(ysim))^2)
tdist[isim]=(r0-min(r1,r2))/r0
}
p1=coneB(y,t(delta),one)$yhat
p2=coneB(y,-t(delta),one)$yhat
r1=sum((y-p1)^2)
r2=sum((y-p2)^2)
r0=sum((y-mean(y))^2)
tstat=(r0-min(r1,r2))/r0
pval=sum(tdist>tstat)/nsim
ans$p1=p1
ans$p2=p2
}else if(k==2){
sse=1:4;theta=matrix(nrow=4,ncol=n);obs=1:4
x1=xmat[,1];x2=xmat[,2]
smat=cbind(x1,x2,y)
basis1=basis2d(smat)
amat1=basis1$amat
smat=cbind(x1,-x2,y)
basis2=basis2d(smat)
amat2=basis2$amat
tdist=1:nsim
for(isim in 1:nsim){
yb=rnorm(n)
ans1=coneA(yb,amat1)
sse1=sum((yb-ans1$thetahat)^2)
ans2=coneA(yb,-amat1)
sse2=sum((yb-ans2$thetahat)^2)
ans3=coneA(yb,amat2)
sse3=sum((yb-ans3$thetahat)^2)
ans4=coneA(yb,-amat2)
sse4=sum((yb-ans4$thetahat)^2)
sse0=sum((yb-mean(yb))^2)
tdist[isim]=(sse0-min(sse1,sse2,sse3,sse4))/sse0
}
ans1=coneA(y,amat1)
theta[1,]=ans1$thetahat
sse[1]=sum((y-ans1$thetahat)^2)
ans2=coneA(y,-amat1)
theta[2,]=ans2$thetahat
sse[2]=sum((y-ans2$thetahat)^2)
ans3=coneA(y,amat2)
theta[3,]=ans3$thetahat
sse[3]=sum((y-ans3$thetahat)^2)
ans4=coneA(y,-amat2)
theta[4,]=ans4$thetahat
sse[4]=sum((y-ans4$thetahat)^2)
sse0=sum((y-mean(y))^2)
tstat=(sse0-min(sse))/sse0
pval=sum(tdist>tstat)/nsim
ans$p1=ans1$thetahat
ans$p2=ans2$thetahat
ans$p3=ans3$thetahat
ans$p4=ans4$thetahat
ans$thetahat=theta[min(obs[sse==min(sse)]),]
ans$pval=pval
}else if(k==3){
sse=1:8;theta=matrix(nrow=8,ncol=n);obs=1:8
smat=xmat
basis1=basis3d(smat)
amat1=basis1$amat
smat=xmat;smat[,3]=-xmat[,3]
basis2=basis3d(smat)
amat2=basis2$amat
smat=xmat;smat[,2]=-xmat[,2]
basis3=basis3d(smat)
amat3=basis3$amat
smat=xmat;smat[,1]=-xmat[,1]
basis4=basis3d(smat)
amat4=basis4$amat
tdist=1:nsim
for(isim in 1:nsim){
yb=rnorm(n)
ans1=coneA(yb,amat1)
sse1=sum((yb-ans1$thetahat)^2)
ans2=coneA(yb,-amat1)
sse2=sum((yb-ans2$thetahat)^2)
ans3=coneA(yb,amat2)
sse3=sum((yb-ans3$thetahat)^2)
ans4=coneA(yb,-amat2)
sse4=sum((yb-ans4$thetahat)^2)
ans5=coneA(yb,amat3)
sse5=sum((yb-ans5$thetahat)^2)
ans6=coneA(yb,-amat3)
sse6=sum((yb-ans6$thetahat)^2)
ans7=coneA(yb,amat4)
sse7=sum((yb-ans7$thetahat)^2)
ans8=coneA(yb,-amat4)
sse8=sum((yb-ans8$thetahat)^2)
sse0=sum((yb-mean(yb))^2)
tdist[isim]=(sse0-min(sse1,sse2,sse3,sse4,sse5,sse6,sse7,sse8))/sse0
}
ans1=coneA(y,amat1)
theta[1,]=ans1$thetahat
sse[1]=sum((y-ans1$thetahat)^2)
ans2=coneA(y,-amat1)
theta[2,]=ans2$thetahat
sse[2]=sum((y-ans2$thetahat)^2)
ans3=coneA(y,amat2)
theta[3,]=ans3$thetahat
sse[3]=sum((y-ans3$thetahat)^2)
ans4=coneA(y,-amat2)
theta[4,]=ans4$thetahat
sse[4]=sum((y-ans4$thetahat)^2)
ans5=coneA(y,amat3)
theta[5,]=ans5$thetahat
sse[5]=sum((y-ans5$thetahat)^2)
ans6=coneA(y,-amat3)
theta[6,]=ans6$thetahat
sse[6]=sum((y-ans6$thetahat)^2)
ans7=coneA(y,amat4)
theta[7,]=ans7$thetahat
sse[7]=sum((y-ans7$thetahat)^2)
ans8=coneA(y,-amat4)
theta[8,]=ans8$thetahat
sse[8]=sum((y-ans8$thetahat)^2)
sse0=sum((y-mean(y))^2)
tstat=(sse0-min(sse))/sse0
pval=sum(tdist>tstat)/nsim
ans$p1=ans1$thetahat
ans$p2=ans2$thetahat
ans$p3=ans3$thetahat
ans$p4=ans4$thetahat
ans$p5=ans5$thetahat
ans$p6=ans6$thetahat
ans$p7=ans7$thetahat
ans$p8=ans8$thetahat
ans$thetahat=theta[min(obs[sse==min(sse)]),]
ans$pval=pval
}
ans
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.