*(1) Create the mean function for Y as a Neural Network with one hidden layer **; 

Data A; keep x1 x2 h1 h2 h3 y p y1-y3 ymean; 
input b0 b1-b2 /  
      c0 c1-c2 /  
      d0 d1-d2 / 
	  a0 a1 a2 a3; 
Do X1 = 1 to 80; 
Do x2 = 1 to 80;
 
**********************
  First Hidden Layer
**********************;
  H1 = b0 + b1*X1+b2*X2; 
     h1=1/(1+exp(-h1));  
  H2 = c0 + c1*X1+c2*X2; 
     h2=1/(1+exp(-h2)); 
  H3 = d0 + d1*X1+d2*X2; 
     h3=1/(1+exp(-h3));
**********************************
  Output Layer 
	with & without activation fn.
*********************************;
Ymean = a0 + a1*H1 +a2*H2 + a3*H3;
p = 1/(1+exp(-Ymean)); ** For a binary target **; 
y1=a1*h1; y2=a2*h2; y3=a3*h3;  
output; 
end;   end;
cards; 
-9 -.5 .9
-15 0 .3 
22 -1 0 
-1.5 2 -1.5 2
 ; 
proc print; 
proc sort; by x1; 


* (2) Plot the surface **; 
proc g3d; plot x1*x2=Ymean; plot x1*x2=p; 
run;  

* (3) Add noise to Y - 3 replications **; 
data next; set A; 
  do replicate = 1 to 3; Y=Ymean + 1.5*normal(123); output; end;

/*(4) Plot the data (with jittering) **; 
  ** This can take a couple of minutes to render **; 
data plot; set next; Xjitter = X1 + .001*ranuni(123); 
proc g3d; scatter Xjitter *X2 = Y/noneedle shape="balloon" size = 0.3; run; 
************************************* */; 

* (5) Estimate the surface **; 
proc nlin data=next; ** estimate a nonlinaer model by least squares **; 
parms b0=-10 b1=-.4 b2=.8   c0=-13  c1=0 c2=.25   d0=20 d1=-.9 d2=.01  a0=-1 a1=3 a2=-1 a3=2.5;
 H1 = b0 + b1*X1+b2*X2; 
     h1=1/(1+exp(-h1));  
 H2 = c0 + c1*X1+c2*X2; 
     h2=1/(1+exp(-h2)); 
 H3 = d0 + d1*X1+d2*X2; 
     h3=1/(1+exp(-h3));
 model Y = a0 + a1*H1 +a2*H2 + a3*H3;
 output out=out1 predicted=pred; 

 run; 

* (6) Plot predictions in 3D **; 
proc g3d; plot x1*x2=Ymean; plot x1*x2=pred; 
run;

 * (7) Compare Ymean to predicted mean **; 
proc gplot; plot (pred ymean)*Ymean/overlay; 
symbol1 v=circle i=none c=black; 
symbol2 v=none i=join c=red; 
run;