## ----coxBids-------------------------------------------------------------
x<-read.csv(header=TRUE,text="value,bid
0.0343,0.0001
0.0767,0.0004
0.1005,0.0005
0.1269,0.0006
0.1481,0.0007
0.2085,0.1780
0.2481,0.2024
0.3035,0.2509
0.4011,0.3748
0.4355,0.3776
0.4566,0.3804
0.4671,0.4180
0.4671,0.4287
0.4934,0.4584
0.5251,0.4773
0.6016,0.5581
0.6201,0.5743
0.6334,0.5610
0.6386,0.5932
0.7312,0.5883
0.7522,0.6555
0.7629,0.6046
0.7894,0.5940
0.8843,0.7017
0.8869,0.7259")
with(x,plot(bid~value,xlab="Valuation $x$",ylab="Bid $b$"))
abline(a=0,b=1,lty="dotted")
abline(a=0,b=1/2)
legend("topleft",c("$b=x$","$\\RNBNE$"),lty=c("dotted","solid"))

## ----libs,cache=FALSE----------------------------------------------------
source("markdownHelp.R")
loadData()

## ----eqBids,fig.width=6.5------------------------------------------------
bids<-read.table("out.csv",col.names=c("bid","x","rho","r"),colClasses='numeric',sep=",")
bids<-within(bids,{rho<-reorder(as.factor(rho),-rho);risk<-reorder(sprintf("$r=%.3g$",r),-r)})
rPal<-brewer.pal(9,"Spectral")
rTheme<-custom.theme(symbol = rPal,fill = rPal,region = rPal,bg = "white", fg = "black",pch = 1:7,lwd=1) ;
rTheme$par.main.text$cex<-1
rTheme$par.main.text$font<-1
rTheme$par.xlab.text$cex<-.8
rTheme$par.ylab.text$cex<-.8
key<-list(title="$\\rho$:",points=FALSE,lines=TRUE,size=1.5,between=1,between.columns=1,corner=c(1,0),columns=2,cex=.8,cex.title=.8,transparent=TRUE)
xyplot(bid ~ x | risk,group=rho,data=bids,layout=c(3,1),xlab="$x$",ylab="$\\gamma(x)$",xlim=c(0,1),type="l",par.settings=rTheme,auto.key=key) + layer(panel.abline(a=c(0),b=c(.5)),packets=1)+layer(panel.abline(a=c(0),b=c(.6)),packets=2)+layer(panel.abline(a=c(0),b=c(.75)),packets=3) + layer(panel.abline(a=0,b=1/2,lty=3,lwd=3))

## ----treatList,results='asis'--------------------------------------------
load(file='dataZ.Rdata')
xx<-within(slist,{
    treat<-sprintf("%s%s%s",omega,ifelse(type=="2nd","II",""),ifelse(b=="0","+",""))
    range<-sprintf("[%s,%d]",omega,as.numeric(omega)+50)
})
xtable(rename(xx[,c(7,6,2:5)],c('treat'="Treatment",'range'='$[\\omegalow,\\omegaup]$','b'='$\\blow$','type'='\\begin{tabular}{c}auction\\\\type\\end{tabular}','indep'='indep.~observations')),align="ccccccc")

## ----convergence,cache=TRUE,fig.width=6.501------------------------------
load(file='dataZ.Rdata')
#
mTheme <- standard.theme(color = FALSE)
mTheme$superpose.line$lwd<-c(1,1,1,1,1,1,2)
#
dataL<-within(dataL[with(dataL,order(value,sid,Period)),],{
    dBid<-bid-c(NA,bid[-length(bid)])
    dBid[Period==1]<-NA
})
xx<-ddply(dataL,~Period+treat+sid,summarise,Acchange=ifelse(Period[1]==1,6,sum(dBid!=0)),Bpchange=mean(abs(dBid)))
xx<-melt(xx,id=c("Period","treat","sid"))
xx<-within(xx,{
           x<-Period
           y<-value
       })
xx2<-within(ddply(dataL,~treat+value+sid,summarise,vchange=mean(abs(dBid),na.rm=TRUE)),{
    variable<-"Cpchange"
    x<-value
    y<-vchange
})
xx3<-within(rbind.fill(xx,xx2),{v<-as.factor(variable);})
levels(xx3$v)<-c("\\# changes / period", "abs.~change / period", "abs.~change / value")
xx3b<-subset(xx3,!is.na(y))
xx4<-with(xx3b,split(xx3b,list(v,treat)))
xx4b<-lapply(xx4,function(xx) with(xx,{
    xrange<-seq(min(x),max(x),length.out=50);
    xx<-as.data.frame(list(y=predict(loess(y~x),newdata=as.data.frame(list(x=xrange))),x=xrange))
    xx$treat<-treat[1];xx$v<-v[1];xx}))
xx4c<-rbind.fill(xx4b)
key<-list(space="top",columns=4,lines=TRUE,points=FALSE)
xyplot(y~x |v,group=treat,type="l",data=xx4c,scale=list(x="free"),xlab="",ylab="",par.settings=mTheme,layout=c(3,1),auto.key=key,prepanel=function(x,y,...) {list(ylim=range(c(0,y)))})

## ----overbidTreat,fig.width=6.5------------------------------------------
qq<-c(25,50,75)
xx<-ddply(subset(dataL,Period>6),~value+treat,function(x) cbind(q=quantile(x$over,qq/100),qq=qq))
key=list(columns=3,lines=TRUE,points=FALSE,text=sprintf("%d\\%% quantile",unique(xx$qq)))
xyplot(q ~ value | treat,group=qq,data=xx,ylab="$b-\\RNBNE$",xlab="$x-\\omegalow$",
       type="l",layout=c(7,1),auto.key=key,
       par.settings=list(superpose.line=list(lty=c(3,1,2))))+layer(panel.refline(h=0),under=TRUE)

## ----overbidLH,cache=TRUE,fig.width=6.5----------------------------------
meanBin <- function (x) {
    bo<-with(x,aggregate(over ~ gid,FUN=mean)$over>0)
    bu<-!bo
    p<-1-pbeta(.5,1+sum(bo),1+sum(bu))
    odds<-p/(1-p)
    c(p=p,odds=odds,logI=10^(trunc(log10(odds)*2)/2))
}
pOver<-ddply(subset(dataL,Period>6),~treat+value,meanBin)
logIPool<-1/meanBin(subset(dataL,Period>6 & value==0 &!(treat %in% c("0+","50II+"))))[["logI"]]
r<-trunc(log10(range(pOver$odds)))
probs<-seq(r[1],r[2])
xyplot(odds ~ value| treat,data=pOver,t=c("a","o"),layout=c(7,1),
       panel=function(...) {
           panel.refline(h=probs)
           panel.abline(h=0)
           panel.xyplot(...)
       },
       yscale.components=yscale.components.log10.3,
       xlab="$x-\\omegalow$",scales=list(y=list(log=10)),
       ylab="\\begin{tabular}{c}\\footnotesize$o_\\qq$\\\\\\footnotesize overb.~more likely than underb.\\end{tabular}")
#
minOverOdds<-min(subset(pOver,value==50)[["logI"]])
minOverSmall<-signif(minUnderOddsPlus<-min(subset(pOver,value==0 & treat %in% c("0+","50II+"))[["logI"]]),1)
maxUnderSmall<-1/max(subset(pOver,value==0 & !(treat %in% c("0+","50II+")))[["logI"]])

## ----slopeBid-------------------------------------------------
dd<-subset(dataL,Period>6 & Type=="1st")
dd.00p<-within(subset(dd,x==0),compare<-min_eingabe!=0)
dd.00p.jags<-slopeDD(dd.00p)
#
dd.50p<- within(subset(dd,x==50),compare<-min_eingabe!=0)
dd.50p.jags<-slopeDD(dd.50p)
#
dd.25p<- within(subset(dd,x!=50),compare<-x== -25)
dd.25p.jags<-slopeDD(dd.25p)
save(file="slopeBid.Rdata",dd.00p.jags,dd.50p.jags,dd.25p.jags)

## ----slopeBid2,cache=TRUE,results='asis'---------------------------------
load("slopeBid.Rdata")
all.dd<-list(dd.00p.jags,dd.50p.jags,dd.25p.jags)
xx<-as.data.frame(t(sapply(all.dd,function(x) {
    z<-combine.mcmc(x)[,"beta[4]"]
    zOdd<-mean(z>0)/(1-mean(z>0))
    c(x[["summary"]][["quantiles"]]["beta[4]",c(3,1,5)],mean(z<=0),mean(z>=0))
})))
xx.stats<-sapply(all.dd,function(x) c(x$mcse$sseff["beta[4]"],psrf=x$psrf$psrf["beta[4]",1]))
xx<-cbind(c(0,50,"all except 0"),c("$\\blow\\neq0$","$\\blow\\neq0$","$\\omegalow = -25$"),xx)
colnames(xx)<-c("$\\omegalow$","$C$","$\\hat{\\beta}_3$","~~95\\%~-\\hspace{-2ex}","\\hspace{-2ex}C.I.~~","$\\Pr(\\beta_3\\le0)$","$\\Pr(\\beta_3\\ge0)$")
xtable(xx,digits=4,align=c("cccccccc"))
#
dd.stats<-cbind(xx[,1],as.data.frame(t(xx.stats)))
colnames(dd.stats)<-c("$\\omegalow$","eff. size","psrf")
#
h3slope<-sprintf("%.2f",xx[1,3])
h4slope<-sprintf("%.2f",xx[2,3])
h5slope<-sprintf("%.2f",xx[3,3])
h3perc<-sprintf("%.2f",xx[1,7]*100)
h4perc<-sprintf("%.2f",xx[2,7]*100)
h5perc<-sprintf("%.2f",xx[3,7]*100)
ssize<-dd.00p.jags$summary$nchain*dd.00p.jags$sample
dd.sseff<-signif(min(xx.stats[1,]),3)
dd.psrf<-round(max(xx.stats[2,]),5)
dd.ssize<-dd.00p.jags$summary$nchain*dd.00p.jags$sample

## ----indivBids,cache=TRUE------------------------------------------------
xx<-ddply(subset(dataL,Period>6),~sid+treat,function(x) {
    est<-lm(bid~value,data=x)
    c(coef(est),R2=summary(est)[["r.squared"]])})
names(xx)<-gsub("[()]","",names(xx))
xx<-within(xx,{ alpha<- - Intercept;streat<-factor(ifelse(!treat %in% c("50II+","0+"),"other 1st",as.character(treat)))})
xx2<-dlply(xx,~streat,function(x) kde(x[,c("value","alpha")]))
r2indivBid<-sprintf("%.4g",median(xx$R2),2)
#
panels<-c("0+","50II+","other 1st")
## ----mixBayes-------------------------------------------------
runMixture(); ## generates "mixG.Rdata"
load("mixG.Rdata") 
#(mixG.j.all,mixG.time,file="mixG.Rdata")
mix.All<-ldply(mixG.j.all,combine.mcmc)
mix.All.z<-within(mix.All,{
    levels(treat)<-c(levels(treat),"other 1st");
    treat[! (treat %in% c("0+","50II+"))]<-"other 1st"})
gg<-grep("^g",names(mix.All.z))
zFreq<-dlply(mix.All.z,~treat,function(x) apply(x[,gg],2,median))
#
zStats<-list()
for (tr in panels) 
    xx2[[tr]][["points"]]<-subset(cbind(cat=TYPES,freq=zFreq[[tr]],data.frame(getStatsFile(tr)$mu)),cat!="RAND")

## ----indivBidsFig,fig.width=6.5001---------------------------------------
library(plotrix)
layout(rbind(1:4),width=c(.1,.3,.3,.3))
par(mar=c(4,3,2,0),cex=1)
plot(NULL,xlim=1:2,ylim=1:2,xaxt='n',yaxt='n',ylab="$\\alpha$",xlab="",frame.plot=FALSE)
par(mar=c(4,0,2,0))
myRed<-hcl(0,60,60)
myTRed<-hcl(0,60,20)
qq<-sapply(c("other 1st","0+","50II+"),function(x) {
    plot(NULL,xlim=c(.1,1.2),ylim=c(-10,15),main=x,ylab="aaa",xlab="$\\beta$",yaxt=ifelse(x=="other 1st","s","n"));abline(h=0,v=c(.5,1),col="#A0A0A0")
    if(!is.null(xx2[[x]][["points"]])) with(xx2[[x]][["points"]],{points(X2,-X1,lwd=4,col=myRed,cex=8*sqrt(freq));thigmophobe.labels(X2,-X1,cat,col=myTRed,cex=.5)})
    plot(xx2[[x]],add=TRUE,labcex=.7)
})


## ----mixBayesEval--------------------------------------------------------
load("mixG.Rdata")
mixVar<-function(PAT) {
    x<-mix.All[,c("treat",grep(PAT,colnames(mix.All),value=TRUE))]
    names(x)[-1]<-TYPES
    x
}
mix.KAll<-mixVar("^kS\\[")
mix.gAll<-mixVar("^g\\[")
mix.mAll<-mix.KAll;PAT<-"^g\\[" ## <- here we take ci, not delta
mix1.All<-subset(mix.mAll,treat!="50II+")
mix1NN.All<-subset(mix.mAll,! (treat %in% c("50II+","0+")))
#
m.Stat <- function(x,ROUND=1) {
    x.sum<-apply(x[,-1],2,sum)
    round(100*x.sum/sum(x.sum),ROUND)
}
all.Stat<-m.Stat(mix1.All)
notM25.Stat<-m.Stat(subset(mix1.All,treat!="-25"))
t.Stat<-ddply(mix.mAll,~treat,m.Stat)
tNN.Stat<-m.Stat(mix1NN.All)
rownames(t.Stat)<-t.Stat[["treat"]]
#
x<-ddply(mix.mAll,~treat,m.Stat,10)
t.Rank<-numcolwise(function(x) rank(-x))(x)
rownames(t.Rank)<-x[["treat"]]
md25<-c("first","second","third","fourth")[t.Rank["-25","MD"]]
mdFirstName<-rownames(subset(t.Rank,MD==1))
mdFirstPerc<-t.Stat[mdFirstName,"MD"]
#
mix.Rel<-mix.mAll[,-1]/apply(mix.mAll[,-1],1,sum)
mix.Rel[["treat"]]<-mix.mAll[["treat"]]
mixM.Rel<-melt(mix.Rel,id="treat")
#

## ----mixBayesEval2,dependson="mixBayesEval",fig.width=6.5001-------------
TYPES2<-c("RAND","MD","BNE","FLEX")
my.stats<-function(x,...) list(stats=quantile(x,c(.05,.25,.5,.75,.95)),n=0,conf=NULL,out=NULL)
mixM2.Rel<-within(mixM.Rel,factor(variable<-factor(variable,TYPES2)))
bwplot(  value ~ treat | variable,data=mixM2.Rel,layout=c(4,1),ylab="relative frequency",stats=my.stats,scales=list(x=list(rot=90)),cex=.5)
##
sample<-mixG.j.all[[1]]$sample
thin<-mixG.j.all[[1]]$thin
burnin<-mixG.j.all[[1]]$burnin
nchain<-mixG.j.all[[1]]$summary$nchain

## ----mStats1-------------------------------------------------------------
load("mixG.Rdata")
mStats <- list()
for (tr in levels(dataL$treat)) 
    mStats[[tr]]<-getStatsFile(tr)

## ----mdStats,fig.width=6.5-----------------------------------------------
z<-within(ldply(mStats,function(x) cbind(cat=TYPES,as.data.frame(x$mu[,]))),
          {id<-relevel(relevel(factor(.id),"0+"),"-25")
           cat<-factor(cat,TYPES2)
           tCat<-relevel(factor(ifelse(cat==TYPES2[1],TYPES2[1],paste(TYPES2[-1],collapse=", "))),"RAND")})
save(z,file="z.Rdata")
xyplot(V1~V2 | tCat,group=id,data=z,xlab="$\\mathrm{mean}(\\beta_{1,i,c})$",
       ylab="$\\mathrm{mean}(\\beta_{0,i,c})$",scale=list(x="sliced",y="sliced"),
       auto.key=list(columns=7),layout=c(2,1),
       par.settings=list(superpose.symbol=list(pch=c(17,10,1,6,5,18,4))),
       panel=function(...) {
           panel.refline(v=0);
           panel.refline(v=1,h=0);
           panel.text(c(.72,1,0),c(0,-2.5,24),c("(\\BNE)","(\\MD)","(\\RA)"),srt=c(0,90,90),cex=.7)
           panel.xyplot(...)
       })

## ----dd.stats,results='asis'---------------------------------------------
xtable(dd.stats,digits=c(0,0,0,5))

## ----maxTab,results='asis',cache=FALSE-----------------------------------
sFun<-function(x) {
    a <- 100*round(quantile(x,c(.025,.5,.975)),3)
    sprintf("\\AA{%g}{%g}{%g}",a[1],a[2],a[3])
}
mixM.g<-melt(mix.gAll,id="treat")
xx<-ddply(mixM.g,.(treat,variable),summarize,r=sFun(value))
names(xx)[1]<-"treatment"
xx<-within(xx,variable<-factor(variable,TYPES2))
xtable(dcast(xx,treatment~variable),align=c("cccccc"))

## ----psrfStat,results='asis'---------------------------------------------
extStat <- function(n1,n2,pat,i=NULL) {
   x<-rbind.fill(lapply(mixG.j.all,function(x) {
       z<-x[[n1]][[n2]]
       if(!is.null(i)) z<-z[,i]
       as.data.frame(as.list(z))
   }))
   rownames(x)<-names(mixG.j.all)
   x<-x[,grep(pat,colnames(x))]
   colnames(x)<-TYPES[as.integer(gsub(pat,"",colnames(x)))]
   x
}
PAT<-"^g\\."
sseff<-extStat("mcse","sseff",PAT)
psrf<-extStat("psrf","psrf",PAT,i=1)
xtable(cbind(treatment=rownames(psrf),psrf[,TYPES2]))

## ----sseffStat,results='asis'--------------------------------------------
xtable(sseff[,TYPES2],digits=0)

## ------------------------------------------------------------------------

