# # svdfit.r # (c) Neil Gershenfeld 4/2/03 # SVD fits and errors # npts <- 100 nfits <- 100 std <- 0.5 # # SVD error # x = matrix(runif(npts, min=0, max=1),npts,1) zeta = rnorm(npts, mean=0, sd=std) y = 2 + 3*x + zeta M = cbind(matrix(1,npts,1),x) SVD = svd(M) U = SVD\$u W = SVD\$d V = SVD\$v ans = V %*% diag(1/W) %*% t(U) %*% y err1 = sqrt(sum(V[1,]^2/W^2)) * std err2 = sqrt(sum(V[2,]^2/W^2)) * std cat("SVD: a =",ans[1]," astd =",err1," b =",ans[2]," bstd =",err2,"\n") # # bootstrap error # x = matrix(runif(npts, min=0, max=1),npts,1) zeta = rnorm(npts, mean=0, sd=0.5) y = 2 + 3*x + zeta aboot = bboot = array(0,nfits) for (i in 1:nfits) { index = ceiling(runif(npts, min=0, max=npts)) xboot = x[index] yboot = y[index] M = cbind(matrix(1,npts,1),xboot) SVD = svd(M) U = SVD\$u W = SVD\$d V = SVD\$v coeff = V %*% diag(1/W) %*% t(U) %*% yboot aboot[i] = coeff[1] bboot[i] = coeff[2] } amean = mean(aboot) astd = sqrt(var(aboot)) bmean = mean(bboot) bstd = sqrt(var(bboot)) cat("bootstrap: a =",amean," astd =",astd," b =",bmean," bstd =",bstd,"\n") # # ensemble error # aens = bens = array(0,nfits) for (i in 1:nfits) { x = matrix(runif(npts, min=0, max=1),npts,1) zeta = rnorm(npts, mean=0, sd=0.5) y = 2 + 3*x + zeta M = cbind(matrix(1,npts,1),x) SVD = svd(M) U = SVD\$u W = SVD\$d V = SVD\$v coeff = V %*% diag(1/W) %*% t(U) %*% y aens[i] = coeff[1] bens[i] = coeff[2] } amean = mean(aens) astd = sqrt(var(aens)) bmean = mean(bens) bstd = sqrt(var(bens)) cat("ensemble: a =",amean," astd =",astd," b =",bmean," bstd =",bstd,"\n") boot = hist(c(aboot,bboot), seq(1,4,.05), plot=FALSE) ens = hist(c(aens,bens), seq(1,4,.05), add=TRUE, plot=FALSE) plot(boot\$mids, boot\$counts, type="b", pch=24, xlab="a,b", ylab="histogram") lines(ens\$mids, ens\$counts, type="b", lty="dashed", pch=21)