%%% GLOBAL VARIABLES FOR YOU TO FIDDLE WITH %%%%%

%%% ASSUMES numcases AND inputstandev HAVE BEEN SET  %%%%
%%% ASSUMES makeallvecs HAS BEEN CALLED %%%%

tiny=exp(-100);
inputstandev=4;
prior=ones(9^4,1)/(9^4); %%a uniform prior

%%% now generate some random inputs with gaussian distributions %%%

data=inputstandev*randn(numcases,2);  %%there are 2 inputs to the net.

targets = .001 + .998*(rand(1,numcases)<applyweights(wteacher, data)); 
        %% not quite binary to stop the log blowing up.
%% the random teacher defines the right answer! %%
%% targets is a row vector with one entry per training case;

reptargets = repmat(targets,9^4,1); %%now we have one row per gridpoint

predictions = applyweights(allwvecs, data);
%% every weightvector makes its own predictions %%

predictions=predictions*(1-2*tiny)+tiny;  %%eliminates 1's and 0's

costs = - sum(reptargets.*log(predictions)+(1-reptargets).*log(1-predictions),2);
%% computes the sum of the cross-entropy costs on all training cases.
%% sums the rows to get total prediction cost for each possible weightvector.


posterior = prior.*exp(-costs);
posterior = posterior/sum(posterior);

[bestposterior bestnet] = max(posterior);
bestpredictions = predictions(bestnet,:);

bayespredictions = sum(predictions.*repmat(posterior,1,numcases));

bayestrainerr = - sum(targets.*log(bayespredictions) + ...
                      (1-targets).*log(1-bayespredictions))  ...
                + sum(targets.*log(targets) + ...
                      (1-targets).*log(1-targets));
%% we subtract off the entropy of the targets so that we can see how
%% much worse we are doing than the best possible performance.

fprintf(1,'Creating new training and test sets but with same teacher net \n');
fprintf(1,'Bayes error in bits per case on training set %4.4f \n', ...
        bayestrainerr/(numcases*log(2)));

figure(1);
%%subplot(1,2,1);
hinton([targets;bayespredictions;bestpredictions]', 1);
title(sprintf('OUTPUTS ON TRAINING DATA \n teacher, bayes, best'));
drawnow;

figure(3);
hist(posterior,100);
title(sprintf('histogram of posterior probabilities \n even the biggest values are small'));

%%%%%%%%%%%%%%%%% END OF TRAINING %%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%% NOW REUSE THE SAME VARIABLES FOR TESTING %%%%%%%%%%%%%
testnumcases=100;

data = inputstandev*randn(testnumcases,2);
targets =  .001 + .998*(rand(1,testnumcases)<applyweights(wteacher, data)); 

predictions = applyweights(allwvecs, data);
predictions=predictions*(1-2*tiny)+tiny;  %%eliminates 1's and 0's
posterior = posterior/sum(posterior);

%% we already decided which net is best during training.
bestpredictions = predictions(bestnet,:);

bayespredictions = sum(predictions.*repmat(posterior,1,testnumcases));

bayestesterr = - sum(targets.*log(bayespredictions)+...
                    (1-targets).*log(1-bayespredictions)) ...
                + sum(targets.*log(targets)+...
                    (1-targets).*log(1-targets));

fprintf(1,'Bayes error in bits per case on test set %4.4f \n', ...
            bayestesterr/(testnumcases*log(2)));

besttesterr = - sum(targets.*log(bestpredictions)+...
                    (1-targets).*log(1-bestpredictions)) ...
                + sum(targets.*log(targets)+...
                    (1-targets).*log(1-targets));

fprintf(1,'Best net error in bits per case on test set %4.4f \n', ...
            besttesterr/(testnumcases*log(2)));

figure(2);
hinton([targets;bayespredictions;bestpredictions]', 1);
title(sprintf('OUTPUTS ON TEST DATA \n teacher, bayes, best'));