*** SAS file generated by sastangle ***; options nocenter nodate nonumber ls=80 formdlim=''; *** Code chunk number 1 ***; data test; p1 = .1; p2 = .2; p3 = .3; do i = 1 to 10000; x = uniform(0); mycat1 = (x ge 0) + (x gt p1) + (x gt p1 + p2) + (x gt p1 + p2 + p3); mycat2 = rantbl(0, .5, .4, .05); mycat3 = rand("TABLE", .3, .3, .4); output; end; run; *** Code chunk number 2 ***; proc freq data=test; tables mycat1 mycat2 mycat3; run; *** Code chunk number 3 ***; data test; intercept = -1; beta = .5; do i = 1 to 5000; xtest = normal(12345); linpred = intercept + (xtest * beta); prob = exp(linpred)/ (1 + exp(linpred)); ytest = uniform(0) lt prob; output; end; run; *** Code chunk number 4 ***; ods select parameterestimates; proc logistic data=test; model ytest(event='1') = xtest; run; ods select all; *** Code chunk number 5 ***; data sim; sigbsq=4; beta0=-2; beta1=1.5; beta2=0.5; beta3=-1; n=1500; do i = 1 to n; x1 = (i lt (n+1)/2); randint = normal(0) * sqrt(sigbsq); do x2 = 1 to 3 by 1; x3 = uniform(0); linpred = beta0 + beta1*x1 + beta2*x2 + beta3*x3 + randint; expit = exp(linpred)/(1 + exp(linpred)); y = (uniform(0) lt expit); output; end; end; run; *** Code chunk number 6 ***; options ls=64; ods select parameterestimates; proc nlmixed data=sim qpoints=50; parms b0=1 b1=1 b2=1 b3=1; eta = b0 + b1*x1 + b2*x2 + b3*x3 + bi1; mu = exp(eta)/(1 + exp(eta)); model y ~ binary(mu); random bi1 ~ normal(0, g11) subject=i; predict mu out=predmean; run; ods select all; *** Code chunk number 7 ***; ods select parameterestimates covparms; proc glimmix data=sim order=data method=laplace; nloptions maxiter=100 technique=dbldog; model y = x1 x2 x3 / solution dist=bin; random int / subject=i; run; ods select all; *** Code chunk number 8 ***; data test; p1=.15; p2=.25; corr=0.4; p1p2=corr*sqrt(p1*(1-p1)*p2*(1-p2)) + p1*p2; do i = 1 to 10000; cat=rand('TABLE', 1-p1-p2+p1p2, p1-p1p2, p2-p1p2); y1=0; y2=0; if cat=2 then y1=1; else if cat=3 then y2=1; else if cat=4 then do; y1=1; y2=1; end; output; end; run; *** Code chunk number 9 ***; options ls = 68; proc corr data=test; var y1 y2; run; *** Code chunk number 10 ***; data simcox; beta1 = 2; beta2 = -1; lambdat = 0.002; *baseline hazard; lambdac = 0.004; *censoring hazard; do i = 1 to 10000; x1 = normal(0); x2 = normal(0); linpred = exp(-beta1*x1 - beta2*x2); t = rand("WEIBULL", 1, lambdaT * linpred); * time of event; c = rand("WEIBULL", 1, lambdaC); * time of censoring; time = min(t, c); * time of first?; censored = (c lt t); * 1 if censored; output; end; run; *** Code chunk number 11 ***; options ls = 68; ods select censoredsummary parameterestimates; proc phreg data=simcox; model time*censored(1) = x1 x2; run; *** Code chunk number 12 ***; data mh; burnin=5000; numvals=5000; thin=10; x = normal(0); do i = 1 to (burnin + (numvals * thin)); y = normal(0) + x; switchprob = min(1, exp(-y**4 + x**4) * (1 + abs(y))**3 * (1 + abs(x))**(-3)); if uniform(0) lt switchprob then x = y; * if we don't change x, the previous value is kept in the data step by default-- no code needed; if (i gt burnin) and mod(i-burnin, thin) = 0 then output; * only start saving if we're past the burn-in period, then keep every 10th to thin; end; run; *** Code chunk number 13 ***; ods select none; proc kde data=mh; univar x / out=mhepdf; run; ods exclude none; *** Code chunk number 14 ***; data mh2; set mhepdf; true = 6.809610784**(-1) * exp(-value**4) * (1 + abs(value))**3; run; *** Code chunk number 15 ***; goptions hsize = 4 in vsize = 4 in; legend1 position=(bottom center inside) across=1 noframe label=none value=(h=1.5); axis1 label=(angle=90 "Density"); symbol1 i=j l=1 w=3 c=black; symbol2 i=j l=21 w=3 c=black; proc gplot data=mh2; plot (density true)*value / overlay vaxis = axis1 legend=legend1; label value="x" density="Simulated" true="True"; run; *** Code chunk number 16 ***; data test; do sim = 1 to 1000; do i = 1 to 500; x = ranexp(42); output; end; end; run; *** Code chunk number 17 ***; ods select none; ods output conflimits=ttci; proc ttest data=test h0=1 sides=2; by sim; var x; run; ods exclude none; *** Code chunk number 18 ***; data summtt; set ttci; lower = (lowerclmean > 1); upper = (upperclmean < 1); run; proc means data=summtt mean; var lower upper; run; *** Code chunk number 19 ***; data dips; do sim = 1 to 10000; do diploma = 1 to 650; randorder = uniform(0); output; end; end; run; proc sort data=dips; by sim randorder; run; *** Code chunk number 20 ***; data match; set dips; student = mod(_n_, 650); match = (student eq diploma); run; proc print data = match (obs = 6); run; *** Code chunk number 21 ***; proc summary data=match; by sim; var match; output out=summatch sum=totalmatches; run; *** Code chunk number 22 ***; proc freq data=summatch; tables totalmatches / nocum; run; *** Code chunk number 23 ***; proc means data=summatch mean std var; var totalmatches; run; *** Code chunk number 24 ***; data mh; do i = 1 to 10000; prize = rand("TABLE", .333, .333); * Monty puts the prize behind a random door; initial_guess = rand("TABLE", .333, .333); * We make a random initial guess; * if the initial guess is right, Monty can open either of the others; * which means that player can switch to either of the other doors; if initial_guess eq prize then do; new_guess = initial_guess; do until (new_guess ne initial_guess); new_guess = rand("TABLE", .333, .333); end; end; * If the initial guess was wrong, we must switch to the right door; if initial_guess ne prize then new_guess = prize; output; end; run; *** Code chunk number 25 ***; data mh2; set mh; win_by_keep = (initial_guess eq prize); win_by_switch = (new_guess eq prize); run; proc means data=mh2 mean; var win_by_keep win_by_switch; run;