/*  Where did you drill on the square grid? (X,Y)
    Did you find water?  W=(Y,N)          */
%let libname=work; options symbolgen; 

Data a; 
  ** Compute logit of water based on (X,Y); 
 do x = 0 to 100; do Y = 0 to 100; 
  L = -1 + .001*((x-10)**2 - .05*(x-30)**3  
        + 2*(x-40)*(Y-30) - 1.9*(y-30)**2);
  ** Compute probability of water & plot **; 
  p = exp(L)/(1+exp(L));output; end; end; 
proc g3d; plot x*Y=L;plot x*Y=p;
title "Theoretical Probability";  
label L="Logit";
run; 

*(2) Compute 0-1 data for Pr{water}=f(X,Y)
     Create a data set for scoring as well (id=2); 
Data water; set a; W = (ranuni(123)<p); 
  XX=X*X; YY=Y*Y; XY=X*Y; cval = "blue"; 
  if W=0 then cval="red";
  keep W X Y XX YY XY id cval; id=1; output; 
  w=.; id=2; output; 

  
 *(3)  ****  Graphs *****; 
proc g3d data=water; 
  scatter X*Y=W/shape="balloon" size = 0.4 noneedle color=cval;
  goptions targetdevice=pscolor; 
  axis1 length=4 in minor=none; 
  title "Water Data"; 
 
proc gplot data=water; plot x*Y=W/haxis=axis1 vaxis=axis1; 
  where id=1; 
  symbol1 v=dot i=none c="red" h=.3; 
  symbol2 v=dot i=none c="blue" h=.3;
  title2 c=red "0 = dry well " c=blue " 1 = water "; 
run;  

     ****  Take 20% for training data, 
           output score data too ***; 
data &libname..watertrain; set water; if w ne .; 
 if ranuni(12345) < .2; * about 20%; 
data &libname..waterscore; set water; if w=.; 
proc print data=&libname..watertrain(obs=10); 
proc print data=&libname..waterscore(obs=10); 
proc gplot data=&libname..watertrain; 
  plot x*Y=W/haxis=axis1 vaxis=axis1; 
  title "Training Data";
run;