##################################################
## simple lift function.
##################################################
#--------------------------------------------------
#plots lift function
liftf = function(yl,phatl,dopl=TRUE) {
if(is.factor(yl)) yl = as.numeric(yl)-1
oo = order(-phatl)
sy = cumsum(yl[oo])/sum(yl==1)
if(dopl) {
   par(mai=c(.8,1.2,.5,.5))
   ii = (1:length(sy))/length(sy)
   plot(ii,sy,type="l",lwd=2,col="blue",xlab="% tried",ylab="% of successes",cex.lab=2)
   abline(0,1,lty=2)
}
return(sy)
}

#deviance loss function
loss = function(y,phat,wht=0.0000001) {
#y should be 0/1
#wht shrinks probs in phat towards .5, don't log 0!
   if(is.factor(y)) y = as.numeric(y)-1
   phat = (1-wht)*phat + wht*.5
   py = ifelse(y==1,phat,1-phat)
   return(-2*sum(log(py)))
}

##################################################
### Section 1. The Neural Net Model for Numeric Y
###read in the data
zag = read.table("zagat.txt",header=T)
summary(zag)


###standardize the x's
minv = rep(0,3)
maxv = rep(0,3)
zagsc = zag
for(i in 1:3) {
minv[i] = min(zag[[i]])
maxv[i] = max(zag[[i]])
zagsc[[i]] = (zag[[i]]-minv[i])/(maxv[i]-minv[i])
}

### nn library
library(nnet)

###fit nn with just one x=food
set.seed(99)
znn = nnet(price~food,zagsc,size=3,decay=.1,linout=T)


###get fits, print summary,  and plot fit
fznn = predict(znn,zagsc)
plot(zagsc$food,zagsc$price)
oo = order(zagsc$food)
lines(zagsc$food[oo],fznn[oo],col="red",lwd=2)
abline(lm(price~food,zagsc)$coef)

summary(znn) # what does this mean?

### try 5 units
set.seed(99)
znn = nnet(price~food,zagsc,size=5,decay=.1,linout=T)
print(summary(znn))

### all three x's
znn = nnet(price~.,zagsc,size=5,decay=.1,linout=T)
fznn = predict(znn,zagsc)
zlm = lm(price~.,zagsc)
fzlm = predict(zlm,zagsc)
temp = data.frame(y=zagsc$price,fnn=fznn,flm=fzlm)
pairs(temp)
print(cor(temp))


##################################################
### Section 2. Size and Decay


### try four different fits

set.seed(14)
znn1 = nnet(price~food,zagsc,size=3,decay=.5,linout=T)
znn2 = nnet(price~food,zagsc,size=3,decay=.00001,linout=T)
znn3 = nnet(price~food,zagsc,size=50,decay=.5,linout=T)
znn4 = nnet(price~food,zagsc,size=50,decay=.00001,linout=T)
temp = data.frame(price = zagsc$price, food = zagsc$food)
znnf1 = predict(znn1,temp)
znnf2 = predict(znn2,temp)
znnf3 = predict(znn3,temp)
znnf4 = predict(znn4,temp)

### plot the fits

par(mfrow=c(2,2))
plot(zagsc$food,zagsc$price)
lines(zagsc$food[oo],znnf1[oo],lwd=2)
title("size=3, decay=.5")
plot(zagsc$food,zagsc$price)
lines(zagsc$food[oo],znnf2[oo],lwd=2)
title("size=3, decay=.00001")
plot(zagsc$food,zagsc$price)
lines(zagsc$food[oo],znnf3[oo],lwd=2)
title("size = 50, decay = .5")
plot(zagsc$food,zagsc$price)
lines(zagsc$food[oo],znnf4[oo],lwd=2)
title("size = 50, decay = .00001")

##################################################
### Section 3. Iterative Fitting and Random Starting Values

### you can control the number of iterations
set.seed(99)
znn3 = nnet(price~food,zagsc,size=50,decay=.5,linout=T)
znn3 = nnet(price~food,zagsc,size=50,decay=.5,linout=T,maxit=20)
znn3 = nnet(price~food,zagsc,size=50,decay=.5,linout=T,maxit=1000)

znnf3 = predict(znn3,temp)
par(mfrow=c(1,1))
plot(zagsc$food,zagsc$price)
lines(zagsc$food[oo],znnf3[oo],lwd=2)


### starting values are random !!!
## you need not converge to the same place

set.seed(23)
temp = nnet(price~food,zagsc,size=2,decay=.001)
summary(temp)

temp = nnet(price~food,zagsc,size=2,decay=.001)
summary(temp)

set.seed(23)
temp = nnet(price~food,zagsc,size=2,decay=.001)
summary(temp)


##################################################
### Section 4. How Does it Work?

### How does it work
## Zagat fit
print(summary(znn))
x = zagsc$food
y = zagsc$price

z1 = 4.35 -0.24 *x
z2 = -7.42 +21.41*x
z3 = -9.93 +13.28*x

f1 = 12.09*exp(z1)/(1+exp(z1))
f2 = 10.7*exp(z2)/(1+exp(z2))
f3 = 22.74*exp(z3)/(1+exp(z3))

oo = order(x)
plot(x,y-12.33)

lines(x[oo],f1[oo],col=2) 
lines(x[oo],f2[oo],col=3) 
lines(x[oo],f3[oo],col=4) 
lines(x[oo],(f1+f2+f3)[oo],col=5)

###how would you fit a bump

set.seed(23)
x = runif(1000)
x = sort(x)
y = exp(-80*(x-.5)*(x-.5)) + .05*rnorm(1000)
plot(x,y)
df = data.frame(y=y,x=x)

plot(x,y)
sz = 3
#Try various decay values.
for(i in 1:20) {
   nnsim = nnet(y~x,df,size=sz,decay = 1/2^i,linout=T,maxit=1000)
   simfit = predict(nnsim,df)
   lines(x,simfit,col=i,lwd=3)
   print(i)
   readline()
}

set.seed(99)
nnsim = nnet(y~x,df,size=3,decay=1/2^12,linout=T,maxit=1000)
thefit = predict(nnsim,df)
plot(x,y)
lines(x,thefit,col="blue",lwd=3,cex.axis=1.5,cex.lab=1.5)


z1 =  5.26 - 13.74*x
z2 = -6.58 + 13.98*x
z3 = -9.67 + 17.87

F = function(x) {return(exp(x)/(1+exp(x)))}

f1 = 2.21*F(z1)
f2 = 7.61*F(z2)
f3 = -5.40*F(z3)

 
rx=range(x)
ry = range(c(f1,f2,f3,y))
plot(rx,ry,type="n",xlab="x",ylab="fit",cex.axis=2,cex.lab=2)
points(x,y)
lines(x,f1,col=1,lwd=2)
lines(x,f2,col=2,lwd=2)
lines(x,f3,col=3,lwd=2)
lines(x,f1+f2+f3,col=4,lwd=4)

##################################################
### 5. Neural Nets for Binary Y

####################
### Beer data
##Fit nn to beer data , predict gender from nbeer
beer = read.csv("nbeer.csv")
beer$gender = as.factor(beer$gender)
print(summary(beer))

beer$nbeer = (beer$nbeer-min(beer$nbeer))/(max(beer$nbeer)-min(beer$nbeer))
set.seed(99)
nnbeer = nnet(gender~nbeer,beer,size=5,decay=.01,maxit=1000)

print(summary(nnbeer))

##get fits and plot
pbeer = predict(nnbeer,beer)
print(dim(pbeer)) #it is an array

plot(beer$nbeer,as.numeric(beer$gender)-1,xlab="nbeer",ylab="p female")
oo = order(beer$nbeer)
lines(beer$nbeer[oo],pbeer[oo],col=2,lwd=2)


####################
###Tabloid data

##read in data
tab = read.table("tabdat9n20.txt",header=T)
tab$purchase = as.factor(tab$purchase)

tab = tab[,1:3]
tab$nTab = tab$nTab/81
tab$moCbook = tab$moCbook/50
print(summary(tab))


#The usual three set stuff
nob = nrow(tab)
n1 = floor(.5*nob)
n2 = floor(.25*nob)
n3 = nob - n1 - n2
set.seed(19)
perm = sample(1:nob,nob)
set1 = tab[perm[1:n1],]
set2 = tab[perm[(n1+1):(n1+n2)],]
set3 = tab[perm[(n1+n2+1):nob],]

##lets' look at one and see if there is any fit
set.seed(99)
tempnn = nnet(purchase~.,set1,size=20,decay=.1,maxit=10000)
nnout = predict(tempnn,set2)[,1]
boxplot(nnout~set2$purchase)
dev.copy2pdf(file="check-tab-out.pdf",height=8,width=12)

#test loss
ltemp = loss(set2$purchase,nnout,wht=.001)
cat("ltemp: ",ltemp,"\n")

##let's try a bunch of decay valuse
#This takes a while !!
#I’m fitting length(decv)*nstart neural nets
#each with 10,000 observations.

if(file.exists("lossl.RData")) {
   cat("******reading in loosl from file lossl.RData\n")
   load("lossl.RData")
} else {
   cat("******running loop to compute lossl\n")
   decv = c(.5,.1,.01,.005,.0025,.001)
   nstart = 10
   lossl = list()
   
   set.seed(99)
   for(i in 1:length(decv)) {
      temploss<-rep(0,nstart)
      for(j in 1:nstart){
         cat("on dev: ",i," and start: ",j,"\n")
         tempnn = nnet(purchase~.,set1,size=20,decay=decv[i],maxit=10000)
         nnout = predict(tempnn,set2)[,1]
         temploss[j] = loss(set2$purchase,nnout,wht=.0001)
      }
      lossl[[i]] = temploss
   }
   save(lossl,decv,nstart,file="lossl.RData")
}

names(lossl) = as.character(decv)
boxplot(lossl,cex.lab=2,cex.axis=2)
dev.copy2pdf(file="tab-oos.pdf",height=10,width=14)

##combine sets 1 and 2, fit, then predict on 3

set12 = rbind(set1,set2)
nns12 = nnet(purchase~.,set12,size=20,decay =.01,maxit=10000)

nnfit12 = predict(nns12,set12)[,1]
nnfit3 = predict(nns12,set3)[,1]

##plot lift in and out of sample
par(mfrow=c(1,1))
sy12 = liftf(set12$purchase,nnfit12,dopl=FALSE)
ii12 = (1:length(sy12))/length(sy12)
sy3 = liftf(set3$purchase,nnfit3,dopl=FALSE)
ii3 = (1:length(sy3))/length(sy3)
plot(c(0,1),c(0,1),type="n",xlab="% used",ylab="% found",
               cex.axis=1.5,cex.lab=1.5)
lines(ii12,sy12,col="red",lty=2,lwd=2)
lines(ii3,sy3,col="blue",lty=3,lwd=2)
abline(0,1,lty=3)
legend("topleft",legend=c("in sample","out of sample"),
           col=c("red","blue"),lwd=c(2,2),lty=c(2,3))


dev.copy2pdf(file="tab-lift.pdf",height=8,width=14)