Post-Selection Inferenceの切断分布を求める処理を書いてみました。

J. Lee and Dennis L. Sun and Yuekai Sun and Jonathan E. Taylor, “Exact post-selection inference, with application to the lasso”, The Annals of Statistics, volume 44,
number=3, pages 907-927 (2016}

# X, yを中心化する関数
cent=function(z){
    if(is.matrix(z)){
        p=ncol(z)
        for(j in 1:p)z[,j]=z[,j]-mean(z[,j])
    }
    else z=z-mean(z)
  return(z)
}

# 線形回帰のLassoの係数を求める関数
lasso=function(X,y,lambda){
## X,yは、中心化されていると仮定 
    p=ncol(X)
    out=NULL
    for(j in 1:p){
        SD=sd(X[,j])
        X[,j]=X[,j]/SD
        beta=sum(X[,j]*y)
        if(beta>lambda)out=c(out,(beta-lambda)/SD)
        else if(beta< -lambda)out=c(out,(beta+lambda)/SD)
        else out=c(out,0)
    }
    return(out)
}

# 多面体で表現される区間(モデルと符号で条件付)
intervals=function(X,y,lambda,M,s,k){
    n=nrow(X)
    p=ncol(X)
    m=length(M)
    P=X[,M]%*%solve(t(X[,M])%*%X[,M])%*%t(X[,M])
    XX=X[,M]%*%solve(t(X[,M])%*%X[,M])
    A=rbind(1/lambda*t(X[,-M])%*%(diag(n)-P),
        -1/lambda*t(X[,-M])%*%(diag(n)-P),
        -diag(s)%*%solve(t(X[,M])%*%X[,M])%*%t(X[,M])
        )
    b=c(rep(1,p-m)-as.vector(t(X[,-M])%*%XX%*%s),
        as.vector(rep(1,p-m)+t(X[,-M])%*%XX%*%s),
        -as.vector(lambda*diag(s)%*%solve(t(X[,M])%*%X[,M])%*%s)
        )
    eta=(solve(t(X)%*%X)%*%(t(X)))[k,]
    cc=eta/as.vector(t(eta)%*%eta)
    z=(diag(n)-cc%*%t(eta))%*%as.matrix(y)
    Ac=as.vector(A%*%cc)
    Az=as.vector(A%*%z)
    nu.max=0
    nu.min=0
    nu.zero=0
    for(j in 1:p){
      if(Ac[j]>0)nu.max=max(nu.max,(b[j]-Az[j])/Ac[j])
      else if(Ac[j]<0)nu.min=min(nu.min,(b[j]-Az[j])/Ac[j])
      else nu.zero=min(nu.zero,b[j]-Az[j])
    }
    return(c(as.vector(t(eta)%*%as.matrix(y)),nu.min, nu.max,nu.zero))
}

# 10進数を2進数に変換する関数
binary=function(i,m){
  out=NULL
  for(j in 1:m){
    out=c(2*(i%%2)-1,out)
    i=i%/%2
  }
  return(out)
}

# モデルと符号に関する多面体の区間を、符号を動かして、合併させている。
bind.intervals=function(u,v){
  p=length(u); q=length(v)
  u=c(u,Inf); v=c(v,Inf)
  u.state=1; v.state=1
  w=NULL
  i=1; j=1
  while(i<=p||j<=q){
    if(u[i]<v[j]){
      if(v.state==1)w=c(w,u[i])
      i=i+1
      u.state=-u.state
    }
    else if(u[i]>v[j]){
      if(u.state==1)w=c(w,v[j])
      j=j+1
      v.state=-v.state
    }
    else {
      if(i!=p+1&&u.state==v.state)w=c(w,v[j])
      i=i+1; j=j+1
      u.state=-u.state; v.state=-v.state
    }
  }
  return(w)
}
bind.intervals(u,v)

##  実行例
# データ生成
n=100
p=5
X=matrix(rnorm(n*p),n,p)
y=X[,1]*2-X[,2]*3+rnorm(n)*0.4
# 中心化して、Lassoの係数を求め、アクティブ集合とそれぞれの符号を求める
X=cent(X)
y=cent(y)
lambda=40
beta=lasso(X,y,lambda)
M=NULL
ss=NULL
for(j in 1:p){
  if(beta[j]>0){M=c(M,j); ss=c(ss,1)}
  else if(beta[j]<0){M=c(M,j); ss=c(ss,-1)}
}
# 多面体の区間を求める
m=length(M)
L=2^m
print(intervals(X,y,lambda,M,ss,1)[1])
print(M)
print(ss)
S=NULL
for(i in 1:L){
  s=binary(i,m)
  u=intervals(X,y,lambda,M,s,1)
  print(s)
  print(u[2:3])
  S=bind.intervals(S,u[2:3])
}
S
pnorm(-1)