# analyze the spam data using the Gelman el al. prior
data=read.table('spam.txt',header=FALSE,as.is=F)
data = as.matrix(data)
N = dim(data)[1]
n = round(N/2)
ntest = N-n
p = dim(data)[2]

# a single try

rs = sample(1:N,N)
X = data[rs,1:(p-1)]
Y = data[rs,p]
x = X[1:n,]
y = Y[1:n]
xtest = X[(n+1):N,]
ytest = Y[(n+1):N]

beta.plasso = plasso_logistic(x,y)
beta.alasso = alasso_logistic(x,y)

xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.alasso
PPS.alasso = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
muhat = exp(xb)/(1+exp(xb))
CRPS.alasso = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.plasso
PPS.plasso = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
muhat = exp(xb)/(1+exp(xb))
CRPS.plasso = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest

cbind(PPS.alasso,PPS.plasso)
cbind(CRPS.alasso,CRPS.plasso)
sum(abs(beta.alasso)==0)
sum(abs(beta.plasso)==0)

cbind(beta.alasso,beta.plasso)

# average over many random partitions

N = dim(data)[1]
n = round(N/2)
ntest = N-n
p = dim(data)[2]

nrep = 10
ppsalasso = ppsplasso = crpsalasso = crpsplasso = nzealasso = nzeplasso = rep(0,nrep)

for (irep in 1:nrep) {

rs = sample(1:N,N)
X = data[rs,1:(p-1)]
Y = data[rs,p]
x = X[1:n,]
y = Y[1:n]
xtest = X[(n+1):N,]
ytest = Y[(n+1):N]

beta.plasso = plasso_logistic(x,y)
beta.alasso = alasso_logistic(x,y)

xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.alasso
ppsalasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
muhat = exp(xb)/(1+exp(xb))
crpsalasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest
nzealasso[irep] = sum(abs(beta.alasso)==0)

xb = cbind(as.matrix(rep(1,ntest)),xtest)%*%beta.plasso
ppsplasso[irep] = (-ytest%*%xb+ sum(log(1+exp(xb))))/ntest
muhat = exp(xb)/(1+exp(xb))
crpsplasso[irep] = (sum((muhat[ytest==0])^2)+sum((1-muhat[ytest==1])^2))/ntest
nzeplasso[irep] = sum(abs(beta.plasso)==0)

}
cbind(mean(ppsalasso),mean(ppsplasso))
cbind(mean(crpsalasso),mean(crpsplasso))
cbind(mean(nzealasso),mean(nzeplasso))

par(mfrow=c(1,3))
aLasso = ppsalasso; pLasso = ppsplasso; 
PPS = cbind(aLasso,pLasso)
boxplot(PPS)
title('PPS')

aLasso = crpsalasso; pLasso = crpsplasso; 
CRPS = cbind(aLasso,pLasso)
boxplot(CRPS)
title('CRPS')

aLasso = nzealasso; pLasso = nzeplasso; 
NZE = cbind(aLasso,pLasso)
boxplot(NZE)
title('NZE')



