diff options
author | ranke <ranke@5fad18fb-23f0-0310-ab10-e59a3bee62b4> | 2006-05-10 15:44:14 +0000 |
---|---|---|
committer | ranke <ranke@5fad18fb-23f0-0310-ab10-e59a3bee62b4> | 2006-05-10 15:44:14 +0000 |
commit | 513dfbdcdda94a901b5901b486ff5500c7d158b1 (patch) | |
tree | fefbf7daadbd7da71add3ed63b3d3b07c4c8e4df /R/chemCal.R | |
parent | 8d30b2cd951c992e4f9aa3055054091e18b8b4f0 (diff) |
The inverse prediction works in a variety of cases and is
tested with Examples 7 and 8 from Massart!
I need to compare with the DIN and draper examples, and finish
the package vignette.
git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@6 5fad18fb-23f0-0310-ab10-e59a3bee62b4
Diffstat (limited to 'R/chemCal.R')
-rw-r--r-- | R/chemCal.R | 95 |
1 files changed, 0 insertions, 95 deletions
diff --git a/R/chemCal.R b/R/chemCal.R deleted file mode 100644 index fab7db4..0000000 --- a/R/chemCal.R +++ /dev/null @@ -1,95 +0,0 @@ -calm <- function(data) -{ - y <- data[[2]] - x <- data[[1]] - m <- lm(y ~ x) - s <- summary(m) - if (s$coefficients[1,4] > 0.05) - { - m <- lm(y ~ x - 1) - s <- summary(m) - m$intercept <- FALSE - } else { - m$intercept <- TRUE - } - class(m) <- "calm" - m$yname <- names(data)[[1]] - m$xname <- names(data)[[2]] - return(m) -} -predict.calm <- predict.lm -print.calm <- print.lm -summary.calm <- summary.lm -plot.calm <- function(x,..., - xunit="",yunit="",measurand="", - level=0.95) -{ - m <- x - x <- m$model$x - y <- m$model$y - newdata <- data.frame(x = seq(0,max(x),length=250)) - pred.lim <- predict(m, newdata, interval = "prediction",level=level) - conf.lim <- predict(m, newdata, interval = "confidence",level=level) - if (xunit!="") { - xlab <- paste("Concentration in ",xunit) - } else xlab <- m$xname - if (yunit=="") yunit <- m$yname - if (measurand!="") { - main <- paste("Calibration for",measurand) - } else main <- "Calibration" - plot(1, - xlab = xlab, - ylab = yunit, - type = "n", - main = main, - xlim = c(0,max(x)), - ylim = range(pred.lim) - ) - points(x,y, pch = 21, bg = "yellow") - matlines(newdata$x, pred.lim, lty = c(1, 4, 4), - col = c("black", "red", "red")) - matlines(newdata$x, conf.lim, lty = c(1, 3, 3), - col = c("black", "green4", "green4")) - - legend(min(x), - max(pred.lim, na.rm = TRUE), - legend = c("Fitted Line", "Confidence Bands", - "Prediction Bands"), - lty = c(1, 3, 4), - lwd = 2, - col = c("black", "green4", "red"), - horiz = FALSE, cex = 0.9, bg = "gray95") -} -predictx <- function(m,yobs,level=0.95) -{ - s <- summary(m) - xi <- m$model$x - yi <- m$model$y - n <- length(yi) - p <- length(yobs) - if (p > 1) { - varyobs <- var(yobs) - } else { - varyobs <- 0 - } - if (!m$intercept) { - b1 <- summary(m)$coef["x","Estimate"] - varb1 <- summary(m)$coef["x","Std. Error"]^2 - xpred <- mean(yobs)/b1 - varxpred <- (varyobs + xpred^2 * varb1) / b1^2 - sdxpred <- sqrt(varxpred) - } else - { - b0 <- summary(m)$coef["(Intercept)","Estimate"] - b1 <- summary(m)$coef["x","Estimate"] - S <- summary(m)$sigma - xpred <- (mean(yobs) - b0)/b1 - sumxxbar <- sum((xi - mean(xi))^2) - yybar <- (mean(yobs) - mean(yi))^2 - sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5 - } - t <- qt((1 + level)/2,n - 2) - confxpred <- t * sdxpred - - result <- c(estimate=xpred,sdxpred=sdxpred,confxpred=confxpred) -} |