# # Program to Illustrative Singular Value Decomposition # # Draws n points from bivariate normal distribution # and draws graph of ratio of 1st over 2nd # Singular Values # nrow <- 1000 ncol <- 2 nitr <- 39 w <- rep(0,nitr*4) dim(w) <- c(nitr,4) xinc <- 2/(nitr+1) rcor <- .95 # # Create Variance-Covariance Matrix # # 1.0 0.9 # 0.9 1.0 # iii <- 1 while (iii <= nitr) { Sigma <- matrix(c(1,rcor,rcor,1),2,2) # # Call Bivariate Normal with zero mean and # Sigma Var-Cov Matrix, Place in X # X <- mvrnorm(n=nrow,rep(0,2),Sigma) # # Perform Singular Value Decomposition # xsvd <- svd(X) # # The Two Lines Below Put the Singular Values in a # Diagonal Matrix -- The first one creates an # identity matrix and the second command puts # the singular values on the diagonal # Lambda <- diag(ncol) diag(Lambda) <- xsvd$d # # Compute U*LAMBDA*V' for check below # XX <- xsvd$u %*% Lambda %*% t(xsvd$v) # # Compute Fit of SVD -- This is just the sum of squared # error -- Note that ssesvd should be zero! # # i <- 0 j <- 0 ssesvd <- 0 while (i < nrow) { i <- i + 1 j <- 0 while (j < ncol) { j <- j + 1 ssesvd <- ssesvd + (X[i,j] - XX[i,j])**2 } } w[iii,1] <- rcor w[iii,2] <- ssesvd w[iii,3] <- xsvd$d[1] w[iii,4] <- xsvd$d[2] rcor <- rcor - xinc iii <- iii+1 # # End of Big Loop } plot(w[,1],w[,3]/w[,4],xlab="Correlation Between Dimensions",ylab="Ratio of 1st to 2nd Singular Value",pch=16,col="blue") #lines(lowess(w[,1],w[,3]/w[,4],f=.2),lwd=3) mtext(side=3,line=1.5,"Bivariate Normal Random Draws \nRatio of Singular Values",font=2)