# # 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 a = a2 <- b <- b2 <- 0 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 a = a + coeff[1] a2 = a2 + coeff[1]^2 b = b + coeff[2] b2 = b2 + coeff[2]^2 } amean = a/nfits astd = sqrt(a2/nfits - amean^2) bmean = b/nfits bstd = sqrt(b2/nfits - bmean^2) cat("bootstrap: a =",amean," astd =",astd," b =",bmean," bstd =",bstd,"\n") # # ensemble error # a = a2 = b = b2 = 0 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 a = a + coeff[1] a2 = a2 + coeff[1]^2 b = b + coeff[2] b2 = b2 + coeff[2]^2 } amean = a/nfits astd = sqrt(a2/nfits - amean^2) bmean = b/nfits bstd = sqrt(b2/nfits - bmean^2) cat("ensemble: a =",amean," astd =",astd," b =",bmean," bstd =",bstd,"\n") plot(x,y) yfit = coeff[1] + coeff[2]*x lines(x,yfit)