Homework Solutions to Problem 19 (Splus6) & 23 (SAS). ===================================================== (19) Of course we could simulate a multinomial(100,(a,(1+2a)/6, (1-a)/3, 1/2-a)) for a=.1, by: > table(sample(1:4,100,replace=T, prob=c(.1,.2,.3,.4))) 1 2 3 4 5 17 27 51 ## But that is not at all parallelizable, so to do it 10000 ## times would either involve doing a 10000x100 array or ## using a for-loop ! Instead, we simulate n1 ~ Binom(100,.1), n2 ~ Binom(100-n1,.2/.9), n3 ~ Binom(100-n1-n2, .3/.7), and n4=100-n1-n2-n3. nArr <- array(0, dim=c(1e4,4)) nArr[,1] <- rbinom(1e4,100,.1) nArr[,2] <- rbinom(1e4,100-nArr[,1],2/9) nArr[,3] <- rbinom(1e4,100-nArr[,1]-nArr[,2],3/7) nArr[,4] <- 100 - nArr[,1]-nArr[,2]-nArr[,3] T1vec <- (nArr[,1]-rep(10,1e4))^2/10 + (nArr[,2]-rep(20,1e4))^2/20 + (nArr[,3]-rep(30,1e4))^2/30 + (nArr[,4]-rep(40,1e4))^2/40 > summary(T1vec) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00000 1.20833 2.38333 2.99121 4.12500 19.45833 > summary(pchisq(T1vec,3)) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0000000 0.2489937 0.5032553 0.4991340 0.7517211 0.9997802 ### If chi-square with 3 df were correct distribution, these ### values would be: 0, .25, .5, .5, .75, 1 > rowlik <- function(nmat, avec) (nmat * log(cbind(avec,(1+2*avec)/6, (1-avec)/3, 0.5-avec))) %*% rep(1,4) rowlikpr <- function(nmat,avec) (nmat/cbind(avec,(1+2*avec)/6, (1-avec)/3, 0.5-avec)) %*% c(1,1/3,-1/3,-1) rowlikdpr <- function(nmat,avec) (nmat/cbind(avec,(1+2*avec)/6, (1-avec)/3, 0.5-avec)^2) %*% c(-1,-1/9,-1/9,-1) > aold <- rep(0,1e4) anew <- rep(.1,1e4) ctr <- 0 iset <- 1:10000 while( ctr < 20 & length(iset)) { aold[iset] <- anew[iset] anew[iset] <- anew[iset] - rowlikpr(nArr[iset,],aold[iset])/ rowlikdpr(nArr[iset,],aold[iset]) iset <- iset[abs(anew[iset]-aold[iset]) > 1.e-5] ctr <- ctr+1 } summary(rowlikpr(nArr,anew)) Min. 1st Qu. Median Mean 3rd Qu. Max. -0.01989776 0.0000 0.0000 -0.00000189 0.00000003 0.00000428 > sum(abs(rowlikpr(nArr,anew)) > 1.e-5) [1] 1 ### Actually the one remaining case, #1230 was FAR from converged !! > anew[abs(rowlikpr(nArr,anew)) > 1.e-5] [1] -5025.066 > nArr[1230,] [1] 4 6 41 49 > optimize(function(a) -rowlik(nArr[1230,],a), c(1e-6,.49))$min [1] 0.02962868 > anew[1230] <- 0.02962868 ### This is a fairly extreme value !! ### There were also 102 cases where a negative anew value was chosen! > for(i in (1:10000)[anew < 0]) anew[i] <- optimize(function(a) -rowlik(nArr[i,],a), c(1e-6,.49))$min > summary(anew) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0132118 0.0814582 0.1000000 0.1003715 0.1183398 0.2111745 > T2vec <- (nArr[,1]-100*anew)^2/(100*anew) + (nArr[,2]-100*(1/6+anew/3))^2/(100*(1/6+anew/3)) + (nArr[,3]-(100/3)*(1-anew))^2/((100/3)*(1-anew)) + (nArr[,4]-100*(0.5-anew))^2/(100*(0.5-anew)) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0000000 0.2612118 0.5057364 0.4993999 0.7444405 0.9998610 ### Again this is close to the desired values for Uniform[0,1] ### We finish off by testing the difference between empirical ### distribution function values for T1vec and T2vec and ### pchisq(.,3) and pchisq(.,3), at the points 1:10. > T1df <- numeric(10) -> T2df > for(i in 1:10) { T1df[i] <- sum(T1vec <= i)/1e4 T2df[i] <- sum(T2vec <= i)/1e4 } > cbind(T1val=T1df, T1CIwidth=round(1.96*sqrt(T1df*(1-T1df))/100,5), T2val=T2df, T2CIwidth=round(1.96*sqrt(T2df*(1-T2df))/100,5), Chi2 = round(pchisq(1:10,2),4), Chi3 = round(pchisq(1:10,3),4)) T1val T1CIwidth T2val T2CIwidth Chi2 Chi3 [1,] 0.1886 0.00767 0.3835 0.00953 0.3935 0.1987 [2,] 0.4262 0.00969 0.6345 0.00944 0.6321 0.4276 [3,] 0.6100 0.00956 0.7765 0.00817 0.7769 0.6084 [4,] 0.7358 0.00864 0.8650 0.00670 0.8647 0.7385 [5,] 0.8273 0.00741 0.9219 0.00526 0.9179 0.8282 [6,] 0.8897 0.00614 0.9542 0.00410 0.9502 0.8884 [7,] 0.9273 0.00509 0.9739 0.00312 0.9698 0.9281 [8,] 0.9546 0.00408 0.9835 0.00250 0.9817 0.9540 [9,] 0.9710 0.00329 0.9894 0.00201 0.9889 0.9707 [10,] 0.9814 0.00265 0.9932 0.00161 0.9933 0.9814 ## Very interesting results! Shows clearly that T1 is at least very close to Chi-square 3df distribution, while T2 is at least very close to Chi-square 2df distribution. There MAY be discrepancy in each case at distribution argument 1 (but nowehere else), and even there the most likely explanation may be some fault of the likelihood maximizations which resulted in the smallest values of the parameter a . =========================================================== (23) SAS version of the same problem: 337 data simdata (keep= x1-x4 chist1 chi3 a ctr chist2 chi2); 338 array x[4]; array xp[4]; 339 seed=5555551; 340 do i=1 to 10000; 341 chi2 = rangam(seed,1)*2; 342 chi3 = rangam(seed,1.5)*2; 343 x[1] = ranbin(seed, 100,0.1); 344 x[2] = ranbin(seed, 100-x[1], 2/9); 345 x[3] = ranbin(seed, 100-x[1]-x[2], 3/7); 346 x[4] = 100-x[1]-x[2]-x[3]; 347 chist1 = (x[1]-10)**2/10 + (x[2]-20)**2/20 + 348 (x[3]-30)**2/30 + (x[4]-40)**2/40 ; 349 a = .25; 350 gval = 1; 351 ctr = 0; 352 do while ( abs(gval) GT 1e-6 and ctr LT 50); 353 gval = x[1]/a + 2*x[2]/(1+2*a) - x[3]/(1-a) - 354 x[4]/(0.5-a); 355 gpr = -x[1]/(a*a) - 4*x[2]/(1+2*a)**2 - 356 x[3]/(1-a)**2 - x[4]/(0.5-a)**2; 357 a = a - gval/gpr; 358 ctr = ctr+1; 359 if(ctr gt 19 and a lt 0) then a=.1; 360 /* This allows a re-start for cases which would 361 not seem to converge to positive value a . */ 362 end; 363 xp[1] = 100*a; xp[2] = (100/6)*(1+2*a); 364 xp[3] = (100/3)*(1-a); xp[4] = 100*(.5-a); 365 chist1 = 0; chist2 = 0; 366 do j= 1 to 4; 367 chist1 = chist1 + (x[j]-j*10)**2/(10*j); 368 chist2 = chist2 + (x[j]-xp[j])**2/xp[j]; 369 end; 370 output; 371 end; 372 proc means min max mean std skewness kurtosis maxdec=3; 373 var chist1 chist2 chi3 chi2 a ctr ; 374 proc means q1 median q3 p95 maxdec=3; 375 var chist1 chist2 chi3 chi2 a ctr ; 376 run; Variable Minimum Maximum Mean Std Dev Skewness Kurtosis --------------------------------------------------------------------- chist1 0 23.608 3.037 2.504 1.773 4.890 chist2 0 20.061 2.009 2.025 2.130 7.081 chi3 0.002 22.275 3.050 2.497 1.600 3.612 chi2 0.001 18.271 2.018 1.994 1.949 5.592 a 0.010 0.216 0.100 0.028 0.207 0.062 ctr 2.000 10.000 4.693 0.594 -0.483 2.742 --------------------------------------------------------------------- Lower Upper Variable Quartile Median Quartile 95th Pctl --------------------------------------------------------- chi2 0.578 1.425 2.806 5.989 chi3 1.218 2.406 4.184 8.038 chist1 1.233 2.383 4.125 7.933 a 0.080 0.099 0.118 0.147 ctr 4.000 5.000 5.000 5.000 chist2 0.608 1.407 2.736 5.881 --------------------------------------------------------- ### This shows reasonably close agreement; also, there were no failures of convergence in a belonging to (0,.22) . NOTE: the close correspondence in all statistics between chist1 and chi3, and between chist2 and chi2, as expected. (But what `close' means in the case of skewness and kurtosis is considerably weaker than in other statistics, because these quantities involve higher moments and thus show more variability.