# Chapter 5

#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))

#ARE ALL THE GENERAL 10 parameter model the same?

################################
#10 parameter Binomial
h=y[,1]/(y[,1]+y[,2])
m=1-h
f=y[,3]/(y[,3]+y[,4])
c=1-f
p=as.matrix(cbind(h,m,f,c))
nll.bin= -sum(y*log(p))

#single high
d.=(h-f)/(1-f)
g=f

ht.h=d+(1-d)*g
ht.f=g

#double high
d=(h-f)
g=f/(1-d)

ht2.h=d+(1-d)*g
ht2.f=(1-d)*g

#sigdet

dprime=qnorm(h)-qnorm(f)
crit=-qnorm(f)

sd.h=pnorm(dprime-crit)
sd.f=pnorm(-crit)

######################################
######################################
#General High Threshold, 7 parameters

#negative log likelihood for general threshold model, 1 condition
#det=c(ds,dn)
#y=c(hit,miss,fa,cr)
#g=guessing
nll.ght.1=function(g,y,det)
{
ds=det[1]
dn=det[2]
p=1:4 # reserve space
p[1]=ds+(1-ds)*g #probability of a hit
p[2]=1-p[1] # probability of a miss
p[3]=(1-dn)*g # probability of a false alarm
p[4] = 1-p[3] #probability of a correct rejection
return(-sum(y*log(p)))
}

#zdet is qlogis(c(ds,dn))
nll.ght=function(zdet,dat)
{
det=plogis(zdet)
return(
optimize(nll.ght.1,interval=c(0,1),y=dat[1,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[2,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[3,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[4,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[5,],det=det)\$objective ) }

zdet=rep(0,2) #starting values of zdet
gen.ht=optim(zdet,nll.ght,dat=y)
plogis(gen.ht\$par)

##############################################
#Double High Threshold Model, 6 parameters
nll.2ht=function(d,dat)
{
det=rep(d,2)
return(
optimize(nll.ght.1,interval=c(0,1),y=dat[1,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[2,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[3,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[4,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[5,],det=det)\$objective ) }

dht=optimize(nll.2ht,interval=c(0,1),dat=y)

############################################
#High Threshold Model, 6 parameters

nll.ht=function(d,dat)
{
det=c(d,0)
return(
optimize(nll.ght.1,interval=c(0,1),y=dat[1,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[2,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[3,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[4,],det=det)\$objective+
optimize(nll.ght.1,interval=c(0,1),y=dat[5,],det=det)\$objective ) }

ht=optimize(nll.ht,interval=c(0,1),dat=y)

##############################################
#Low Threshold Model, 7 parameters

amount.of.bias=function(b,p)
ifelse(b>0,(1-p)*b,p*b)

#det=c(ds,dn)
#y=c(hits,misses,fa’s,cr’s)
nll.lt.1=function(b,y,det)
{
ds=det[1]
dn=det[2]
p=1:4 # reserve space
p[1]=ds+amount.of.bias(b,ds) #probability of a hit
p[2]=1-p[1] # probability of a miss
p[3]=dn+amount.of.bias(b,dn) # probability of a false alarm
p[4] = 1-p[3] #probability of a correct rejection
return(-sum(y*log(p)))
}

#zdet is qlogis(c(ds,dn))
nll.lt=function(zdet,dat)
{
det=plogis(zdet)
return(
optimize(nll.lt.1,interval=c(-1,1),y=dat[1,],det=det)\$objective+
optimize(nll.lt.1,interval=c(-1,1),y=dat[2,],det=det)\$objective+
optimize(nll.lt.1,interval=c(-1,1),y=dat[3,],det=det)\$objective+
optimize(nll.lt.1,interval=c(-1,1),y=dat[4,],det=det)\$objective+
optimize(nll.lt.1,interval=c(-1,1),y=dat[5,],det=det)\$objective ) }

zdet=c(-.5,-4) #starting values of zdet
lt=optim(zdet,nll.lt,dat=y)
plogis(lt\$par)

########################################################3
#Signal Detection, 7 parameters
#negative log likelihood of free-variance signal detection model
#par=c(d,sigma)
#y=c(hits,misses,fa’s,cr’s)
nll.fvsd.1=function(c,y,par)
{
d=par[1]
sigma=par[2]
p=1:4 # reserve space
p[1]=1-pnorm(c,d,sigma) #probability of a hit
p[2]=1-p[1] # probability of a miss
p[3]=1-pnorm(c) # probability of a false alarm
p[4] = 1-p[3] #probability of a correct rejection
return(-sum(y*log(p)))
}

nll.fvsd=function(par,y)
{
return(
optimize(nll.fvsd.1,interval=c(-5,5),y=y[1,],par=par)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[2,],par=par)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[3,],par=par)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[4,],par=par)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[5,],par=par)\$objective
) }

par=c(1,1)
fvsd=optim(par,nll.fvsd,y=y)

#sd, equal variance, 6 parameters
nll.evsd=function(par,y)
{
pars=c(par,1)
return(
optimize(nll.fvsd.1,interval=c(-5,5),y=y[1,],par=pars)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[2,],par=pars)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[3,],par=pars)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[4,],par=pars)\$objective+
optimize(nll.fvsd.1,interval=c(-5,5),y=y[5,],par=pars)\$objective
) }

evsd=optimize(nll.evsd,interval=c(-5,5),y=y)

#model order: binomial,ght,ht,2ht,lt,fvsd,evsd

dev=2*c(nll.bin,gen.ht\$value,dht\$objective,ht\$objective,lt\$value,fvsd\$value,evsd\$objective)

#bin to GHT
dev[2]-dev[1]
#bin to LT
dev[5]-dev[1]
#bin to FVSD
dev[6]-dev[1]

#ght to dht
dev[3]-dev[2]
#ght to ht
dev[4]-dev[2]
#fvsd to evsd
dev[7]-dev[6]

#AIC
num.par=c(10,7,6,6,7,7,6)
penalty=2*num.par
aic=dev+penalty

aic=aic-min(aic)