# Chapter 6

#Sample Data
dat.signal=c(3,7,14,24,32,12)
dat.noise=c(12,16,21,15,10,6)
NS=sum(dat.signal)
NN=sum(dat.noise)

#General 10-parameter multinumial

#Estimate Parameters
parest.signal=dat.signal/NS
parest.noise=dat.noise/NN
#why are there 12 estimates for a 10-parameter model?

#evaluate nll for model
#negative log likelihood for one stimulus
nll.mult.1=function(p,y) -sum(y*log(p))

nll.general=nll.mult.1(parest.signal,dat.signal)+
nll.mult.1(parest.noise,dat.noise)

#Signal Detection Model
#d' + 5 criteria.

sdprob.1=function(mean,sd,bounds)
{
cumulative=c(0,pnorm(bounds,mean,sd),1)
p.ij=diff(cumulative)
return(p.ij)
}

nll.sigdet=function(par,y)
#negative log likelihood of free-variance signal detection
#for confidence interval paradigm
#par=d,sigma,bounds
#y=y_(1,s),..,y_(I,s),y_(1,n),..,y_(I,n)
{
I=length(y)/2
d=par[1]
sigma=par[2]
bounds=par[3:length(par)]
p.noise=sdprob.1(0,1,bounds)
p.signal=sdprob.1(d,sigma,bounds)
nll.signal=nll.mult.1(p.signal,y[1:I])
nll.noise=nll.mult.1(p.noise,y[(I+1):(2*I)])
return(nll.signal+nll.noise)
}

y=c(dat.signal,dat.noise)
par=c(1,1,-1,0,1,2,3)
nll.sigdet(par,y)

a=optim(par,nll.sigdet,y=y)
a1=optim(a\$par,nll.sigdet,y=y)

2*(a1\$value-nll.general)
#Signal Detection Model Fits Quite Well :)

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

biz=c(103,2,46,1,7,22)
com=c(80,1,65,3,9,23)

#General Binom:
Nbiz=sum(biz)
Ncom=sum(com)

#General 10-parameter multinumial

#Estimate Parameters
parest.biz=biz/Nbiz
parest.com=com/Ncom
#why are there 12 estimates for a 10-parameter model?

nll.gen.biz=nll.mult.1(parest.biz,biz)
nll.gen.com=nll.mult.1(parest.com,com)

#Storage Retrieval Model (negative log likelihood)
#y=#E1, #E2,..#E6
#par=c(a,r,s,u) logit transformed
nll.sr.1=function(par,y)
{
a=plogis(par[1]) #transform variables to 0,1
r=plogis(par[2])
s=plogis(par[3])
u=plogis(par[4])
p=1:6
p[1]=a*(r+(1-r)*s^2)
p[2]=2*a*(1-r)*s*(1-s)
p[3]=a*(1-r)*(1-s)^2
p[4]=(1-a)*u^2
p[5]=2*(1-a)*u*(1-u)
p[6]=(1-a)*(1-u)^2
-sum(y*log(p))
}

#fit to bizare & common
par=qlogis(c(.8,.5,.1,.1))
biz.sr=optim(par,nll.sr.1,y=biz)
biz.sr.est=plogis(biz.sr\$par)
com.sr=optim(par,nll.sr.1,y=com)
com.sr.est=plogis(com.sr\$par)

#common R model (non-nested, 7 parameters)
#a1,r,s1,u1,a2,s2,u2

nll.common.R=function(par,biz,com)
{
nll.sr.1(par[1:4],biz)+nll.sr.1(par[c(5,2,6,7)],com)
}

par=qlogis(c(.9,.6,.1,.1,.9,.1,.1))
common.R=optim(par,nll.common.R,biz=biz,com=com)
common.R=optim(common.R\$par,nll.common.R,biz=biz,com=com)
common.R=nlm(nll.common.R,ommon.R\$par,biz=biz,com=com)

plogis(common.R\$estimate)