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