#
# 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)