/********************************************************************
GAUSS Code for One Hidden Layer Feedforward Neural Networks

Written by:	Langche Zeng		 
	  	lzeng@fas.harvard.edu
	  	November 30, 1998


This file contains two procedures for training and applying single
hidden layer feedforward neural networks for continuous or binary
output data. The network is trained with or without weight decay
(regularization in the form of Gaussian priors for the parameters.)
The hidden layer activation function is logit, and output layer
activation function is linear if continuous data, logit if binary. A
couple of examples using these procedures are also included. The
training proc requires the Gauss Maximum Likelihood application
module.

The example programs can be executed simply by running this file in
Gauss.  For general use of the nnet procs, strip the examples from the
file, put into the (GAUSSHOME)/src subdirectory, and add the entry

     NNET.src
         NNET_TRAIN  : proc
         NNET_APPLY  : proc

to the maxlik.lcg library file.

The code is provided without guarantee for public, non-commercial
use. The author thanks Gary King for helpful discussions, and NSF for
supporting her visit at Harvard where the code was written. 

*********************************************************************/




/********************** PROC NNET_TRAIN ****************************

-- Purpose: train a neural network (using maxlik)

-- Format: {b,f,g,cov,retcode} = nnet_train(dataset, out_var, hidden,
                              prior_sd, start, terminate, "logfile") 

-- Input: 

dataset - name of data matrix. First column is dependent var, and the
next k columns are independent vars. IMPORTANT: data should be
normalized for training to be effective/efficient. See the example
program for more details.

out_var - scalar, output var type: 1=continuous, 2=binary

hidden	 - scalar, number of hidden neurons
          
prior_sd - scalar, standard deviation of the prior distribution for
the parameters (i.i.d. Gaussain with 0 means). Setting this to 0
causes the default sqrt(2) to be used, and setting this to -1 disables
the use of priors.

start - 2x1 vector, the first element is the starting value for error
standard deviation, and the second the size of the random starting
values for all other parameters (drawn from N(0, start[2]). Setting
this to scalar 0 causes the default {1,.1} to be used. (When
out_type==2, the first element is ignored.)

terminate - scalar, termination condition: the value for the maxlik
global _max_GradTol. Setting it to 0 causes the default 1e-4 to be
used.

logfile - name of the file that logs results from print commands (such
as results from maxprt). Set to "null" if do not want any logfile.

-- Output:  

b - estimated parameters (or most probable parm values). A vector of
length (k+2)*m+2 for continuous data, (k+2)*m+1 for binary data. The
first element is error standard deviation (not present if out_type=2),
the next (k+1)*m are input to hidden parameters, and the remaining m+1
are hidden to output parameters, where k is the number of input vars
and m the number of hidden neurons (set with "hidden".)

f - scalar, function value at convergence (maximum log-likelihood).

g - gradient vector, same dimention as b.

cov - covariance matrix of b.

retcode - return code. Same as that returned by the maxlik procedure. 


-- Globals

Maxlik globals are relevant.  See MAXLIK.SRC for their description. 

*********************************************************************/


#include maxlik.ext;

PROC (5) = NNET_TRAIN(dataset, out_var, hidden, prior_sd, start, 
                      terminate, logfile);

  local b0, start_sigma, start_size, b,f,g,cov,ret;

  declare matrix k,m,prior,out_type;
  
  out_type=out_var;
  m=hidden; 			@ hidden neurons @

  /* prior */

  if prior_sd==-1;
     prior=0;
  elseif prior_sd==0;
     prior=.5;
  else;
     prior=1/prior_sd^2;
  endif;

  /* starting values for estimation */

  if start==0;
     start_sigma=1;		@ starting sigma value @
     start_size=.1;		@ starting parm values: N(0, start_size) @
  else;
     start_sigma=start[1];		 
     start_size=start[2];		
  endif;

  /* termination condition */

  if terminate==0;
     _max_gradtol=.0001;				 
  else;
     _max_gradtol=terminate;
  endif;

  /* parameter estimation */

  __title=" *** One Hidden Layer Feedforward Neural Network Model ***"; 

  k=cols(dataset)-1;	                         @ No. of input vars @

  if out_type==1;
  b0=start_sigma|rndn((k+2)*m+1,1)*start_size;   @ starting values of parms.@
  else;
  b0=rndn((k+2)*m+1,1)*start_size;   @ starting values of parms.@  
  endif;

  if logfile$/="null";
      output file=^logfile reset;
      {b,f,g,cov,ret}=maxprt(maxlik(dataset,0,&llik_nn,b0)); 
      output off;
  else;
      {b,f,g,cov,ret}=maxprt(maxlik(dataset,0,&llik_nn,b0)); 
  endif;

/* return of nnet_train */
   retp(b,f,g,cov,ret);
endp;


 /* llik proc for nnet */

  proc llik_nn(b,dta);                       @ log-likelihood function @

     external matrix k,m,prior, out_type;

     local bmat,gmat,s,s2,xb,zg,x,y,z,w,llik;
     
     w=(prior/2)*sumc(b^2)/rows(dta);     @ weight decay term @ 
     y=dta[.,1];                          @ output var @
     x=ones(rows(dta),1)~dta[.,2:cols(dta)]; @ inputs plus const. n by k+1 @ 
     if out_type==1;
       bmat=reshape(trimr(b,1,m+1),k+1,m);@ inputs to hidden coeffs, k+1 by m @
       gmat=trimr(b,(k+1)*m+1,0);         @ hidden to output coeffs, m+1 by 1 @
       s=b[1];                              @ error sd; @
       s2=s^2;                              @ error variance @
       xb=x*bmat;                           @ input for hidden layer, n by m @
       z=1./(1+exp(-xb));                   @ hidden neurons, n by m @
       zg=(ones(rows(z),1)~z)*gmat;         @ input for output layer, n by 1 @
       llik=-0.5*(ln(s2)+(y-zg).*(y-zg)/s2)-w;
     else;
       bmat=reshape(trimr(b,0,m+1),k+1,m);@ inputs to hidden coeffs, k+1 by m @
       gmat=trimr(b,(k+1)*m,0);           @ hidden to output coeffs, m+1 by 1 @
       xb=x*bmat;                           @ input for hidden layer, n by m @
       z=1./(1+exp(-xb));                   @ hidden neurons, n by m @
       zg=(ones(rows(z),1)~z)*gmat;         @ input for output layer, n by 1 @
       llik=-y.*ln(1+exp(-zg))-(1-y).*ln(1+exp(zg))-w;
     endif;
     retp(llik);
  endp;

/********************** PROC NNET_APPLY ****************************

-- Purpose: generate in-sample or out-of-sample estimated output from
            a trained neural network model.

-- Format:  yhat = nnet_apply(dataset, out_type, hidden, b)  

-- Input: 

dataset - name of data matrix. Should be in the same format as the
training dataset (i.e., normalized to the same scale; first column is
dependent var, the next k columns are independent vars.)

out_type - scalar, output type: 1=continuous, 2=binary. (Must set to
the same value as for NNET_TRAIN.)

hidden - scalar, number of hidden neurons (Must set to the same value
as for NNET_TRAIN.)

b - the estimated parameter vector from NNET_TRAIN.

-- Output: 

yhat - estimated y (if out_type==1) or prob(y=1) (if out_type==2) from
a trained neural network.

*********************************************************************/

PROC (1) = NNET_APPLY(dataset, out_type, hidden, b);

  local k, m, x, bmat, gmat, xb, z, zg, yhat;
  
  m=hidden;
  k=cols(dataset)-1;
  x=ones(rows(dataset),1)~dataset[.,2:cols(dataset)]; @inputs plus const.@
  if out_type==1;
   bmat=reshape(trimr(b,1,m+1),k+1,m);   @ inputs to hidden coeffs, k+1 by m @
   gmat=trimr(b,(k+1)*m+1,0);		@ hidden to output coeffs, m+1 by 1 @
   xb=x*bmat;    			@ input for the hidden layer, n by m @
   z=1./(1+exp(-xb));			@ hidden neurons, n by m @
   zg=(ones(rows(z),1)~z)*gmat;		@ input for the output layer, n by 1 @
   yhat=zg;				@ nnet prediction of y @
  else;
   bmat=reshape(trimr(b,0,m+1),k+1,m);   @ inputs to hidden coeffs, k+1 by m @
   gmat=trimr(b,(k+1)*m,0);		@ hidden to output coeffs, m+1 by 1 @
   xb=x*bmat;    			@ input for the hidden layer, n by m @
   z=1./(1+exp(-xb));			@ hidden neurons, n by m @
   zg=(ones(rows(z),1)~z)*gmat;		@ input for the output layer, n by 1 @
   yhat=1./(1+exp(-zg));		@ nnet prediction of prob(y=1)@
  endif;
  retp(yhat);
endp;


/**********************************************************************
**
**  Sample Programs Using NNET_TRAIN and NNET_APPLY
**
**  -- Purpose: To illustrate the use of the NNET procedures with
**              simulated data. 
**
**********************************************************************/



/************** example 1: continuous data *******************/

library maxlik;
#include maxlik.ext;
maxset;

rndseed 2078;

n=1000;				@ sample size @
n1=500; 			@ training set size @
sigma=0.2; 			@ error term standard deviation @

x=rndn(n,2); 			@ input vars @
e=rndn(n,1).*sigma;		@ error @
y0=ones(n,1)+x[.,1]+x[.,1]^2+x[.,1].*x[.,2]; @systematic part@
y=y0+e;				@ continuous dep. var data @

/* normalize data for effective/efficient training. This step is VERY
IMPORTANT, as the priors, starting values for parms etc make sense only
when the data are properly scaled. The estimated parameters themselves
are relatively more meaningful, too when all input vars are of similar
scales. */

mx=meanc(x)';
sx=stdc(x)';
my=meanc(y);
sy=stdc(y);

xn=(x-mx)./sx;		@ normalized x, having mean 0 and s.d. 1 @
yn=(y-my)./sy;	        @ normalized y. Do this ONLY FOR continuous data @

/* dataset for NNET, using the first n1 obs for training and the rest
for test */ 

data_train=trimr(yn~xn,0,n1);
data_test=trimr(yn~xn,n1,0);

/* training a nnet with 5 hidden neurons, using default values for
"prior", "start" and "terminate", and sending printout to "test1.out" */

{b,f,g,cov,retcode} = nnet_train(data_train, 1, 5, 0, 0, 0, "test1.out"); 


/* generating estimated y values from the trained nnet, for both
training and test data */

yhat_train= nnet_apply(data_train, 1, 5, b);
yhat_test= nnet_apply(data_test, 1, 5, b);

/* convert the estimates to the original y scale */

yhat_train=yhat_train.*sy + my; 
yhat_test=yhat_test.*sy + my; 

/* Now compare nnet performance between training and test data, and
with that from the "true" model. You may want to do other things, like
calculating R^2, comparing with a simple linear model, etc. */

ytrue=y0;			@ best possible est. given error @ 

output file=test1.out on;

"";

"=============";
"mean and s.d. of absolute error, true model (training data part):";
meanc(abs(y[1:n1]-ytrue[1:n1]))~stdc(abs(y[1:n1]-ytrue[1:n1]));
"mean and s.d. of absolute error, true model (test data part):";
meanc(abs(y[(n1+1):n]-ytrue[(n1+1):n]))~stdc(abs(y[(n1+1):n]-ytrue[(n1+1):n]));

"mean and s.d. of absolute error, nnet (training data):";
meanc(abs(y[1:n1]-yhat_train))~stdc(abs(y[1:n1]-yhat_train));
"mean and s.d. of absolute error, nnet (test data):";
meanc(abs(y[(n1+1):n]-yhat_test))~stdc(abs(y[(n1+1):n]-yhat_test));

"true sigma:";
sigma;
"estimated sigma:";
b[1]*sy;
output off;


/******************* example 2: binary data ********************/

library maxlik;
#include maxlik.ext;
maxset;

rndseed 2078;

n=1000;				@ sample size @
n1=500; 			@ traing set size @;
sigma=0.2; 			@ error term standard deviation @

x=rndn(n,2); 			@ input vars @
e=rndn(n,1).*sigma;		@ error @
y0=x[.,1]+2.*x[.,1].*x[.,2];    @ systematic part of "utility function" @
   y=((y0+e) .> 0);             @ binary data @


/* normalize data for effective/efficient training. DO NOT normalize y
for binary data. */

mx=meanc(x)';
sx=stdc(x)';

xn=(x-mx)./sx;		@ normalized x, having mean 0 and s.d. 1 @

/* dataset for NNET, using the first 500 obs for training and the rest
for test */ 

data_train=trimr(y~xn,0,n1);
data_test=trimr(y~xn,n1,0);

/* training a nnet with 5 hidden neurons, using default values for
"prior", "start" and "terminate", and sending printout to "test2.out" */

{b,f,g,cov,retcode} = nnet_train(data_train, 2, 5, 0, 0, 0, "test2.out"); 


/* generating estimated prob(y=1) from the trained nnet, for both
training and test data. No conversion needed here as the original y is
not scaled. */

phat_train= nnet_apply(data_train, 2, 5, b);
phat_test= nnet_apply(data_test, 2, 5, b);

/* Now compare nnet performance between training and test data, and
with that from the "true" model. You may want to do other things, too,
like looking at the mean probabilities for the 1's and 0's, comparing
with a simple logit model, etc. */

ptrue=1-cdfn(-y0/sigma);       @ best possible est. given error @

correct_true_train=((y[1:n1].==1 .and ptrue[1:n1].>0.5) .or
(y[1:n1].==0 .and ptrue[1:n1].<=0.5));

correct_true_test=((y[(n1+1):n].==1 .and ptrue[(n1+1):n].>0.5) .or
(y[(n1+1):n].==0 .and ptrue[(n1+1):n].<=0.5));

correct_train=((y[1:n1].==1 .and phat_train.>0.5) .or (y[1:n1].==0 .and phat_train.<=0.5));
correct_test=((y[(n1+1):n].==1 .and phat_test.>0.5) .or (y[(n1+1):n].==0 .and phat_test.<=0.5));

output file=test2.out on;

"";

"=================";
"percentage of correct prediction, true model: (training data part)";
   sumc(correct_true_train)/rows(correct_true_train);
"percentage of correct prediction, true model: (test data part)";
   sumc(correct_true_test)/rows(correct_true_test);
"percentage of correct prediction, nnet: (training)";
   sumc(correct_train)/rows(correct_train);
"percentage of correct prediction, nnet: (test)";
   sumc(correct_test)/rows(correct_test);

output off;
end;