data loopex4;
seed=11023;
do stickprov=1 to 1000;
    do obs=1 to 5;
    x=100*ranuni(seed);
    output;
    end;
end;
drop obs;
drop seed;
run;

proc summary data=loopex4;
    by stickprov;
    id x;
    output out=WTPmax (drop= _TYPE_ _FREQ_ ) max(x)=;
run;

data WTPmax;
set WTPmax;
    e=(5/((5+1))*100);
    std=SQRT((5*(100*100))/(((5+1)*(5+1))*(5+2)));
    xstd=(x-e)/std;
run;

data WTPmax;
set WTPmax;
    if xstd < -2.5 then xgroup=1;
    if xstd < -1.25 & xstd > -2.5 then xgroup=2;
    if xstd < -0.5 & xstd > -1.25 then xgroup=3;
    if xstd < 0 & xstd > -0.5 then xgroup=4;
    if xstd < 0.5 & xstd > 0 then xgroup=5;
    if xstd < 1 & xstd > 0.5 then xgroup=6;
    if xstd > 1 then xgroup=7;
run;

proc print data=WTPmax;
run;

proc format;
value groupformat
    1 = '<-2.5'
    2 = '-2.5 - -1.25'
    3 = '-1.25 - -0.5'
    4 = '-0.5 - 0'
    5 = '0 - 0.5'
    6 = '0.5 - 1'
    7 = '>1'
    ;
run;

proc sgplot data=WTPmax;
vbar xgroup / stat=percent barwidth=1 limits=both;
format xgroup groupformat.;
run;

proc freq data=WTPmax;
/weight;
by xgroup;
1*2*3*4*5*6*7;
output out=percent;
run;

proc sort data=WTPmax;
      by xgroup;
   run;

   ods graphics on;
   proc freq data=WTPmax order=data;
      tables xgroup / chisq testp=(0.0258 0.09685 0.13578 0.14344 0.20106 0.27451 0.12268)
                    plots(only)=deviationplot(type=dot);
      weight percent;
      by xgroup;
      title 'Chi square of theoretical distribution';
   run;
   ods graphics off;