/*  Report on the Best Method of Proportioning the Excise upon Spirituous Liquors
    by Charles Blagden
    Philosophical Transactions of the Royal Society of London, Vol 80 (1790),pp321-345
*/

options nodate nonumber ps=54 ls=80;

data onea;
do temp=30 to 100 by 5;
do dilute=0 to 50 by 5;
input weight @@;
weight=weight+2000;
dilute0=0;
if dilute=0 then dilute0=1;
output;
end; end;
cards;
487.32 519.98 549.59 573.86 596.65 617.24 636.16 653.54 669.64 684.63 698.41
480.79 513.48 541.96 567.34 590.15 610.80 629.77 647.30 663.48 678.43 692.32
474.18 506.98 535.52 560.83 583.70 604.50 623.42 641.02 657.35 672.37 686.37
467.52 500.33 528.90 554.24 577.16 597.99 617.04 634.68 650.96 666.13 680.25
460.77 493.48 522.10 547.61 570.64 591.50 610.59 628.26 644.68 659.95 674.04
453.84 486.51 515.30 540.88 563.94 584.79 604.07 621.77 638.25 653.55 667.72
446.86 479.75 508.60 534.19 557.23 578.22 597.50 615.26 631.82 647.20 661.45
440.04 472.97 501.87 527.51 550.56 571.48 590.86 608.72 625.41 640.80 655.09
433.37 466.28 495.00 520.65 543.84 564.89 584.23 602.14 618.89 634.30 648.65
426.47 459.18 488.03 513.63 536.91 558.14 577.47 595.43 612.20 627.78 642.17
419.18 451.95 480.83 506.61 529.85 551.10 570.52 588.61 605.32 621.03 635.47
412.02 444.80 473.68 499.59 523.08 544.41 563.80 581.91 598.76 614.48 628.87
404.92 437.72 466.64 492.62 516.20 537.57 556.95 575.20 592.17 607.86 622.30
397.75 430.56 459.51 485.51 509.51 530.51 549.95 568.18 585.12 601.12 615.70
390.64 423.53 452.63 478.59 502.15 523.59 543.08 561.28 578.37 594.45 609.11
;
run;


data oneb;
do temp=30 to 100 by 5;
do dilute=55 to 100 by 5;
input weight @@;
weight=weight+2000;
dilute0=0;
output;
end; end;
cards;
711.19 723.00 733.84 744.19 753.67 762.61 771.26 779.21 786.47 793.26
705.08 716.96 727.78 738.24 747.90 756.97 765.47 773.53 780.75 787.59
699.09 711.02 721.90 732.45 742.18 751.35 759.85 767.78 775.15 782.06
692.97 704.82 715.98 726.38 736.21 745.47 754.13 762.03 769.55 776.40
686.81 698.63 709.92 720.37 730.27 739.52 748.22 756.25 763.67 770.62
680.57 692.44 703.67 714.27 724.20 733.47 742.25 750.32 757.82 764.72
674.31 686.16 697.44 708.18 718.26 727.56 736.26 744.32 751.87 758.82
667.97 680.00 691.22 701.99 712.06 721.47 730.27 738.35 745.93 752.82
661.67 673.68 685.02 695.66 705.87 715.40 724.22 732.42 740.00 746.88
655.19 667.32 678.60 689.34 699.57 709.08 717.95 726.25 733.92 740.83
648.47 660.69 672.02 682.77 693.03 702.57 711.50 719.78 727.49 734.49
641.85 654.19 665.54 676.40 686.77 696.33 705.37 713.69 721.47 728.60
635.38 647.61 659.07 670.09 680.60 690.22 699.10 707.44 715.22 722.32
628.83 641.10 652.56 663.64 674.20 683.79 692.81 701.18 708.91 716.04
622.22 634.38 645.95 657.14 667.61 677.25 686.36 694.76 702.50 709.75
;
run;

data one;
set onea oneb;
run;

data two;
set one;
file 'H:\public_html\data\spiritwt.dat';
put @1 temp 3.0 +4 dilute 3.0 +4 weight 7.2;
run;

quit;

proc glm;
model weight=temp dilute temp*dilute temp*temp dilute*dilute;
output out=sprtout p=yhat r=e;
run;

proc univariate normal;
var e;
run;

proc print;
run;

proc glm data=one;
model weight=temp dilute temp*dilute temp*temp dilute*dilute dilute0;
output out=sprtout0 p=yhat0 r=e0;
run;

proc univariate normal;
var e0;
run;

proc print;
quit;
/*
symbol1 c=black v=:;
proc gplot;
plot e*yhat / vref=0 frame;
run;

proc gchart;
vbar e / space=0;
run;

quit;
quit;