/*------------------------------------------------
|  Task: given trimodal data, find the NARROWEST  |
| interval that contains 50% of the data.         |
 ------------------------------------------------*/
goptions reset=all; 
footnote " "; 
Data Life A; 
Do Age = 0 to 99 by .5; 
Y1 = .03*max(1-age/8,0);  
Y2 = .001*max(age-70,0); 
Y3 = .07*exp(-.004*(age-78)**2); 
fy = (y1+y2+y3+.005)/5.52; 
output Life; age1=age; age = age1 + round (ranuni(123)-.5, .1); drop age1; 
do i=1 to floor(1000*fy); x=ranuni(1827655);  
age = max(age,0.5); 
output A; end; end; 
proc means mean; var fy; run; 

proc gplot data=Life; 
plot (Y1 Y2 Y3)*AGE/overlay vaxis = axis1; 
plot fy*AGE/overlay vaxis = axis1; 
symbol1 v=none i=join; 
run; 

 proc sort data=A; by X;
 data A; set A; id=_n_; 
 keep age id; 
proc means data=A n mean max min q1 q3; Var age;
run;
 proc gplot data=A; plot id*age /href = 65.5  88.5;  
symbol1 v=plus i=none h=0.5; 
footnote "Quartiles 65.5 and 88.5"; 
 ** plot the 50% interval
 (on 2nd run)**;
 proc gchart data=A; 
vbar age/midpoints=0.5 to 99.5 by 4;   ** check trimodality **;

proc sort data=A; by age; 
 data B; set A;
 high=age; low=lag781(age); D= high-low;
 drop age;
 proc sort data=B; by D;
 proc print data=B (obs=5); where D ne .;
 proc sort data=B; by low; 
 data b; set b;
 if d>78 or d<19.41 then do; 
output;  save=high; high=low; 
output; high=save; end; output;  
 proc gplot data=b; plot d*(high low)/overlay href=65.5 88.5 chref=green;
 symbol1 v=dot i=join h=0.5 ; 
 label high="50% age containment interval"; 
 run;