MLP Network in Regression Problem with 2 Inputs, 'demo_2input'

The demonstration program solves a simple 2-dimensional regression problem using an MLP network with Bayesian lerning. The demonstration is discussed completely in the documentation. The data in the demonstration has been used in the work by Christopher J. Paciorek and Mark J. Schervish (2004). The training data has an unknown Gaussian noise and can be seen in the figure 1.

An MLP network with 10 hidden units and one output is created. The weight and bias structure is hierarchical with three levels of hyper-parameters and automatic relevance determination (ARD). The sampling for hyperpameters is done with Gibbs sampling and for weights with Hybrid Monte Carlo sampling. We sampled for 2500 samples from which 250 burn in samples were thrown away. After this We still thinned the sample chain by taking only every 138'th sample. This was done to ensure the independence of all samples. The convergence, independence and explicitly of samples was examined by Markov Chain Monte Carlo diagnostics tools. Figure 2 shows the prediction from the network.

The code needed to get the results is shown below and is treated more carefully in the documentation. The noisy data and the prediction of network is shown in the bottom.

	      function demo_2input

	      % Load the data.
	      data=load('dat.1');
	      x = [data(:,1) data(:,2)];
	      y = data(:,3);
     
	      % Draw the data. The data is not sorted.
	      figure
	      title({'The noisy training data'});
	      [xi,yi,zi]=griddata(data(:,1),data(:,2),data(:,3),-1.8:0.01:1.8,[-1.8:0.01:1.8]');
	      mesh(xi,yi,zi)
	      
	      % Create an MLP network for regression model. function mlp 
	      % initializes weights to zero.
	      nin=size(x,2);
	      nhid=10;
	      nout=1;
	      
	      net = mlp2('mlp2r', nin, nhid, nout);
	      
	      % Create a Gaussian multivariate hierarchical prior with ARD 
	      % for network...
	      net=mlp2normp(net, {{repmat(0.1,1,net.nin) 0.05 0.5 -0.05 1} ...
	                         {0.1 0.05 0.5} ...
                                 {0.1 -0.05 0.5} ...
                                 {1}});

	      % ...and the same for residual.
	      net.p.r = norm_p({0.05 0.05 0.5});

	      % Set the options for Monte Carlo sampling
	      opt=mlp2r_mcopt;
	      opt.repeat=50;
	      opt.plot=0;
	      opt.hmc_opt.steps=40;
	      opt.hmc_opt.stepadj=0.1;
	      hmc2('state', sum(100*clock));

	      % Initialize weights to zero and sample for the first time.
	      % The opt.gibbs = 0 for this round of MC-sampling, so the 
	      % sampling is done only for weigths. This way we can reach the 
	      % equilibrium distribution more quickly.
	      net = mlp2unpak(net, mlp2pak(net)*0);
	      [r1,net1]=mlp2r_mc(opt, net, x, y);

	      % Now sample also for hyper-priors with Gibbs sampling.
	      opt.hmc_opt.stepadj=0.2;
	      opt.hmc_opt.steps=60;
	      opt.repeat=70;
	      opt.gibbs=1;
	      [r2,net2]=mlp2r_mc(opt, net1, x, y, [], [], r1);

	      % Now we should be close to equilibrium already and we can start 
	      % the sampling with windowing and persistence.
	      opt.hmc_opt.stepadj=0.3;
	      opt.hmc_opt.steps=100;
	      opt.hmc_opt.window=5;
	      opt.hmc_opt.persistence=1;
	      [r3,net3]=mlp2r_mc(opt, net2, x, y, [], [], r2);
	      
	      % Sample for the posterior 2500 samples from which every second 
	      % is taken after 41 burn-in samples (the function thin).
	      opt.repeat=60;
	      opt.hmc_opt.steps=100;
	      opt.hmc_opt.stepadj=0.5;
	      opt.nsamples=2500;
	      [r,net]=mlp2r_mc(opt, net3, x, y, [], [], r3);
	      mlp=thin(r,150,25);

	      % Create new data with the right scale checked above
	      figure
	      [p1,p2]=meshgrid(-1.8:0.05:1.8,-1.8:0.05:1.8);
	      p=[p1(:) p2(:)];
	      tms=squeeze(mlp2fwds(mlp,p))';
	      g=mean(tms);
	      
	      %Plot the new data
	      gp=zeros(size(p1));
	      gp(:)=g;
	      mesh(p1,p2,gp);
	      axis on;