is the powerful tool to do Statistical Computing. This is the first R code for me.
calibrate <- function(x,y,newy, confidence=0.95)
{
n = length(x)
if(n != length(y))
stop("x and y must be the same length")
lm.out <- lm(y~x)
lm.sum <- summary(lm.out)
beta <- lm.out$coef
sigma <- lm.sum$sigma
print(sigma)
xbar <- mean(x)
ybar <- mean(y)
m <- length(newy) # no. measurements on the unknown
newybar <- mean(newy)
confP <- (1+confidence)/2
# estimate the unknown x
xhat <- (newybar - beta[1])/beta[2]
ndf <- n-2
tmult <- qt(confP,ndf)
sxx <- sum((x-xbar)^2)
# print.out <- c(t.dist, sxx)
# print(print.out)
g <- tmult / (abs(beta[2]) / (sigma / sqrt(sxx)))
if(g>=1)
{
lo <- -Inf
hi <- Inf
}
else
{
xhatxbar <- xhat-xbar
gfactor <- 1-g^2
updown <- (xhatxbar*g^2 + (tmult * sigma / abs(beta[2])) *
sqrt(xhatxbar^2 / sxx + gfactor*(1/m+1/n)))/gfactor
lo <- xhat - updown
hi <- xhat + updown
a <- c(xhatxbar,gfactor,updown, lo, hi)
}
result <- c(xhat,lo,hi)
names(result) <- c("est","lo","hi")
print(result)
result
}
dipoleFunc <- function(ID, x, y){
plotcolors <- c("red","blue", "grey")
plot(x,y,
,xlab="I [A]", ylab="B [kGauss]"
,xlim=c(0,190), ylim=c(0,10), col=plotcolors[1], lwd=0.1
)
index = ID -1
title(main = substitute(paste("Dipole 2",index," Field Measurement")), font.main=1)
fm <- lm(y~x)
abline(coef(fm), lwd=0.1, col=plotcolors[1])
cat("Dipole 2",ID,"Field measurement\n")
prelow =predict(fm, data.frame(x = x), interval = "prediction", level=.95)[ , "lwr"]
prehigh=predict(fm, data.frame(x = x), interval = "prediction", level=.95)[ , "upr"]
a <- c(prelow,prehigh)
print(a)
curve(predict(fm, data.frame(x = x), interval = "confidence", level=.95)[ , "lwr"], add = TRUE, col=plotcolors[2], lty=2, lwd=0.1)
curve(predict(fm, data.frame(x = x), interval = "confidence", level=.95)[ , "upr"], add = TRUE, col=plotcolors[2], lty=2, lwd=0.1)
curve(predict(fm, data.frame(x = x), interval = "prediction", level=.95)[ , "lwr"], add = TRUE, col=plotcolors[3], lty=1, lwd=0.1)
curve(predict(fm, data.frame(x = x), interval = "prediction", level=.95)[ , "upr"], add = TRUE, col=plotcolors[3], lty=1, lwd=0.1)
axis(1,c(0,25,75,125,175))
# 1: solid, 2: dashed, 3: dotted, 4: dotdash, 5: longdash, 6: twodash
legend(5, 9, c("fitted data","confidence interval", "prediction interval"), cex=1,
col=c(plotcolors), lty=1:2, bty="n");
n <- length(x)
ndf <- n-2
xbar <- mean(x)
ybar <- mean(y)
SSxx <- sum(x^2) - n*xbar*xbar
SSyy <- sum(x^2) - n*ybar*ybar
SSxy <- sum(x*y) - n*xbar*ybar
bhat <- SSxy/SSxx
ahat <- ybar - bhat*xbar
yhat <- ahat + bhat*x
SSR <- sum((yhat-ybar)^2)
SSE <- sum((y - yhat)^2)
SDyhat <- sqrt(SSE/ndf);
SEbhat <- SDyhat/sqrt(SSxx)
SEahat <- SDyhat*sqrt(1/n + xbar*xbar/SSxx);
# test <- c(ahat, bhat, SEahat, SEbhat)
# print(test)
}
dipole <- function()
{
# postscript("dipole.field.measurement.eps")#, pointsize=10)
postscript("/home/jhlee/documents/tex/thesis/eps/chicane/dipole.field.measurement.eps", height=15, width=8)#, pointsize=10)
par(mfrow=c(2,2), mar=c(5.1,4.1,4.1,2.1), oma=c(0,0,0,0))
B1 <- c(1.0624000, 2.1099000, 3.1600000, 4.2130000, 5.2830000,
6.3600000, 7.4270000, 8.4890000, 7.4590000, 6.3930000,
5.3140000, 4.2430000, 3.1880000, 2.1365000, 1.0812000,
1.0643000, 2.1126000, 2.7848000, 3.1590000, 4.2130000,
8.4980000, 6.3930000, 4.2420000, 2.8097000, 2.1367000,
1.0812000)
I1 <- c(17.6822440, 37.6822440, 57.6822440, 77.6822440, 97.6822440,
117.6822440, 137.6822440, 157.6822440, 137.6822440, 117.6822440,
97.6822440, 77.6822440, 57.6822440, 37.6822440, 17.6822440,
17.6822440, 37.6822440, 50.4822440, 57.6822440, 77.6822440,
157.6822440, 117.6822440, 77.6822440, 50.4822440, 37.6822440,
17.6822440
)
B2 <- c(1.0644000, 2.1082000, 2.7785000, 3.1660000, 4.2550000,
5.3340000, 6.4010000, 7.4600000, 8.5050000, 7.4990000,
6.4430000, 5.3700000, 4.2880000, 3.2020000, 2.8140000,
2.1415000, 1.0877000)
I2 <- c(18.2545973, 38.2545973, 51.0545973, 58.2545973, 78.2545973,
98.2545973, 118.2545973, 138.2545973, 158.2545973, 138.2545973,
118.2545973, 98.2545973, 78.2545973, 58.2545973, 51.0545973,
38.2545973, 18.2545973
)
B3 <- c(1.0638000, 2.1079000, 2.7780000, 3.1530000, 4.2150000,
5.2940000, 6.3620000, 7.4230000, 8.4700000, 7.4600000,
6.4020000, 5.3290000, 4.2480000, 3.1860000, 2.8119000,
2.1399000, 1.0862000)
I3 <- c(17.4394640, 37.4394640, 50.2394640, 57.4394640, 77.4394640,
97.4394640, 117.4394640, 137.4394640, 157.4394640, 137.4394640,
117.4394640, 97.4394640, 77.4394640, 57.4394640, 50.2394640,
37.4394640, 17.4394640
)
B4 <- c(1.0624000, 2.1080000, 2.7790000, 3.1520000, 4.2040000,
5.2520000, 6.3250000, 7.3900000, 8.4430000, 7.4200000,
6.3590000, 5.2810000, 4.2330000, 3.1840000, 2.8083000,
2.1351000, 1.0805000)
I4 <- c(16.7124446, 36.7124446, 49.5124446, 56.7124446, 76.7124446,
96.7124446, 116.7124446, 136.7124446, 156.7124446, 136.7124446,
116.7124446, 96.7124446, 76.7124446, 56.7124446, 49.5124446,
36.7124446, 16.7124446
)
dipoleFunc(1,I1, B1)
dipoleFunc(2,I2, B2)
dipoleFunc(3,I3, B3)
dipoleFunc(4,I4, B4)
newy <- c(4.5,4.6,4.55,4.54)
dev.off()
}