# Double & Single HT Models

Paper to read is at http://pcl.missouri.edu/sites/default/files/rouder.etal_.pnas_.2008.pdf

Problems 3.4.1 and 3.5.1

Here are the HT and DHT Analysis for the data:

y=matrix(nrow=5,ncol=4,byrow=T,
c( 404, 96, 301, 199,
348, 152, 235, 265,
287, 213, 183, 317,
251, 249, 102, 398,
148, 352, 20, 480))

###HIGH THRESHOLD ANALYSIS

nll.ht=function(par,y)
{
p=1:4
d=par[1]
g=par[2]
p[1]=d+(1-d)*g
p[2]=1-p[1]
p[3]=g
p[4]=1-p[3]
return(-sum(y*log(p)))
}

nll.ht.gen=function(par,y)
{
nll.ht(par[1:2],y[1,])+
nll.ht(par[3:4],y[2,])+
nll.ht(par[5:6],y[3,])+
nll.ht(par[7:8],y[4,])+
nll.ht(par[9:10],y[5,])
}

nll.ht.one.d=function(par,y)
{
nll.ht(par[c(1,2)],y[1,])+
nll.ht(par[c(1,3)],y[2,])+
nll.ht(par[c(1,4)],y[3,])+
nll.ht(par[c(1,5)],y[4,])+
nll.ht(par[c(1,6)],y[5,])
}

#pick good starting values for convergence!
h=y[,1]/(y[,1]+y[,2])
f=y[,3]/(y[,3]+y[,4])
d=(h-f)/(1-f)
g=f

#general
par=c(d[1],g[1],d[2],g[2],d[3],g[3],d[4],g[4],d[5],g[5])
gen.ht=optim(par,nll.ht.gen,y=y)

#restricted
par=c(mean(d),g)
one.d.ht=optim(par,nll.ht.one.d,y=y)

g2.ht=2*(one.d.ht\$value-gen.ht\$value)

#########################################################

###DOUBLE-HIGH THRESHOLD ANALYSIS

nll.dht=function(par,y)
{
p=1:4
d=par[1]
g=par[2]
p[1]=d+(1-d)*g
p[2]=1-p[1]
p[3]=(1-d)*g
p[4]=1-p[3]
return(-sum(y*log(p)))
}

nll.dht.gen=function(par,y)
{
nll.dht(par[1:2],y[1,])+
nll.dht(par[3:4],y[2,])+
nll.dht(par[5:6],y[3,])+
nll.dht(par[7:8],y[4,])+
nll.dht(par[9:10],y[5,])
}

nll.dht.one.d=function(par,y)
{
nll.dht(par[c(1,2)],y[1,])+
nll.dht(par[c(1,3)],y[2,])+
nll.dht(par[c(1,4)],y[3,])+
nll.dht(par[c(1,5)],y[4,])+
nll.dht(par[c(1,6)],y[5,])
}

#pick good starting values for convergence!
d=(h-f)
g=f/(1-d)

#general
par=c(d[1],g[1],d[2],g[2],d[3],g[3],d[4],g[4],d[5],g[5])
gen.dht=optim(par,nll.dht.gen,y=y)

#restricted
par=c(mean(d),g)
one.d.dht=optim(par,nll.dht.one.d,y=y)

g2.dht=2*(one.d.dht\$value-gen.dht\$value)

#Quick graph

par(pty='s')
range=c(0,1)
plot(f,h,pch=21,cex=1.2,bg='yellow',xlim=range,ylim=range,
ylab="Hit Rate",xlab="False Alarm Rate")
abline(0,1,lty=3)
abline(one.d.ht\$par[1],1-one.d.ht\$par[1],col='darkred',lwd=2)
abline(one.d.dht\$par[1],1,col='darkblue',lwd=2)