LOG FOR IN-CLASS RANDOM-NUMBER DEMONSTRATION, 2/6/08 ===================================================== > Xmat = matrix(sample(1:5,size=1e6,replace=T, prob=c(.1,.2,.4,.25,.05)), ncol=1000) > summary(Xmat) Min. 1st Qu. Median Mean 3rd Qu. Max. 1 2 3 2.949 4 5 > table(Xmat) 1 2 3 4 5 100021 200121 400591 249574 49693 > chisq = function (xnum, pvec) sum((xnum-sum(xnum)*pvec)^2/(sum(xnum)*pvec)) > unix.time({tmp = summary(apply(Xmat,2, function(xcol) chisq(table(xcol),rep(.2,5))))}) [1] 14.00000000 0.04000008 14.00000000 0.00000000 0.00000000 ### USER SYSTEM ELAPSED : now these times are 10 times ### faster than they were a couple of years ago !! > tmp Min. 1st Qu. Median Mean 3rd Qu. Max. 264.3 353.3 378.5 380.1 405.6 517.5 > tmpv = numeric(1000) > unix.time({for(i in 1:1000) tmpv[i] = chisq(table(Xmat[,i]), rep(.2,5))}) [1] 13.04999924 0.05999994 13.00000000 0.00000000 0.00000000 > summary(tmpv) ### same as tmp before ### Now re-do calculation with the true prob vector, ### compare distribution with chisq_4 > p0 = c(.1,.2,.4,.25,.05) > tmpv = apply(Xmat,2, function(xcol) chisq(table(xcol),p0)) > summary(tmpv) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.089 1.836 3.246 3.991 5.258 19.01 > tmphst = hist(tmpv, breaks=c(0, qchisq((1:9)*0.1,4),1.e3),plot=F) > tmphst$counts [1] 97 111 111 104 99 94 91 79 111 103 > chisq(tmphst$counts,rep(.1,10)) [1] 9.56 > 1-pchisq(.Last.value,9) [1] 0.3872644 TIMING RELATED TO VECTORIZATION OF LINEAR-ALGEBRA OPERATIONS ============================================================ > dim(Xmat) [1] 1000 1000 > unix.time(Xmat %*% rep(1,1000)) [1] 0.96999359 0.04999995 1.00000000 0.00000000 0.00000000 > unix.time(apply(Xmat,1,sum)) [1] 3.45000458 0.08000016 3.00000000 0.00000000 0.00000000 > unix.time(for(i in 1:1000) { sm = 0 for (j in 1:1000) sm = sm+Xmat[i,j] }) [1] 3.18 0.00 3.20 0.00 0.00 #### These are times in seconds as of 2008 > unix.time(for(i in 1:5) { sm = 0 for (j in 1:1000) sm = sm+Xmat[i,j] }) [1] 0.02 0.00 0.01 0.00 0.00 ### Again, timing as of 2008 > sum(abs(c(Xmat %*% rep(1,1000)) - apply(Xmat,1,sum))) ## = 0