# NOT RUN {
## From summary data:
EXPR1 <- expression(x/y)
x <- c(5, 0.01)
y <- c(1, 0.01)
DF1 <- cbind(x, y)
RES1 <- propagate(expr = EXPR1, data = DF1, type = "stat",
do.sim = TRUE, verbose = TRUE)
RES1$prop
RES1$sim
## From raw data:
EXPR2 <- expression(x/y)
x <- c(2, 2.1, 2.2, 2, 2.3, 2.1)
y <- c(4, 4, 3.8, 4.1, 3.1, NA)
DF2 <- cbind(x, y)
RES2 <- propagate(expr = EXPR2, data = DF2, type = "raw",
do.sim = TRUE, verbose = TRUE)
RES2$prop
RES2$sim
## Compare to using a multivariate t-distribution
## because of low sample size => larger confidence
## intervals.
RES2b <- propagate(expr = EXPR2, data = DF2, type = "raw",
do.sim = TRUE, verbose = TRUE, dist.sim = "t")
RES2$sim
RES2b$sim
## Example using a recursive function:
## no Taylor expansion possible, only Monte-Carlo.
a <- c(5, 0.1)
b <- c(100, 2)
DAT <- cbind(a, b)
f <- function(a, b) {
N <- 0
for (i in 1:100) {
N <- N + i * log(a) + b^(1/i)
}
return(N)
}
propagate(f, DAT, nsim = 100000)
################# GUM 2008 (1) ########################
## Example in Annex H.1 from the GUM 2008 manual
## (see 'References'), an end gauge calibration
## study. At first, we will only use the first-order
## Taylor expansion.
EXPR3 <- expression(ls + d - ls * (da * the + as * dt))
ls <- c(50000623, 25)
d <- c(215, 9.7)
da <- c(0, 0.58E-6)
the <- c(-0.1, 0.41)
as <- c(11.5E-6, 1.2E-6)
dt <- c(0, 0.029)
DF3 <- cbind(ls, d, da, the, as, dt)
RES3 <- propagate(expr = EXPR3, data = DF3, type = "stat",
do.sim = TRUE, verbose = TRUE,
second.order = FALSE)
RES3$prop
RES3$sim
## propagate: sd.1 = 31.71
## GUM H.1.4/H.6c: u = 32
## Proof that covariance of Monte-Carlo
## simulated dataset is the same as from
## initial data.
RES3$covMat
cov(RES3$datSIM)
all.equal(RES3$covMat, cov(RES3$datSIM))
## Expanded uncertainty GUM H.1.6
## with 16 degrees of freedom.
qt(0.005, 16, lower.tail = FALSE) * 31.71
## propagate: 92.62
## GUM H.1.6: 93
## Second-order terms GUM H.1.7.
RES4 <- propagate(expr = EXPR3, data = DF3, type = "stat",
do.sim = TRUE, verbose = TRUE,
second.order = TRUE)
RES4$prop
RES4$sim
## propagate: sd.2 = 33.91115
## GUM H.1.7: u = 34.
## Also similar to the non-matrix-based approach
## in Wang et al. (2005, page 408): u1 = 33.91115.
## NOTE: After second-order correction, uncertainty is more
## similar to the value obtained from Monte Carlo simulation!
#################### GUM 2008 (2) #################
## Example in Annex H.2 from the GUM 2008 manual
## (see 'References'), simultaneous resistance
## and reactance measurement.
data(H.2)
## This gives exactly the means, uncertainties and
## correlations as given in Table H.2:
colMeans(H.2)
sqrt(colVarsC(H.2))/sqrt(5)
cor(H.2)
## H.2.3 Approach 1 using mean values and
## standard uncertainties:
EXPR6a <- expression((V/I) * cos(phi)) ## R
EXPR6b <- expression((V/I) * sin(phi)) ## X
EXPR6c <- expression(V/I) ## Z
MEAN6 <- colMeans(H.2)
SD6 <- sqrt(colVarsC(H.2))
DF6 <- rbind(MEAN6, SD6)
COV6ab <- cov(H.2) ## covariance matrix of V, I, phi
COV6c <- cov(H.2[, 1:2]) ## covariance matrix of V, I
RES6a <- propagate(expr = EXPR6a, data = DF6, type = "stat",
do.sim = TRUE, verbose = TRUE, use.cov = COV6ab,
second.order = TRUE)
RES6b <- propagate(expr = EXPR6b, data = DF6, type = "stat",
do.sim = TRUE, verbose = TRUE, use.cov = COV6ab,
second.order = TRUE)
RES6c <- propagate(expr = EXPR6c, data = DF6[, 1:2], type = "stat",
do.sim = TRUE, verbose = TRUE, use.cov = COV6c,
second.order = TRUE)
## This gives exactly the same values of mean and sd/sqrt(5)
## as given in Table H.4.
RES6a$prop
RES6b$prop
RES6c$prop
#################### GUM 2008 (3) #######################
## Example in Annex H.3 from the GUM 2008 manual
## (see 'References'), calibration of a thermometer.
## For this example, we can use the predict.lm function
## of R directly to return the error of predicted values.
data(H.3)
LM <- lm(bk ~ I(tk - 20), data = H.3)
## This will give exactly the same values as in H.3.3.
summary(LM)
## This will give exactly the same values as the
## fourth column ("Predicted correction") of Table H.6.
predict(LM)
## This will give exactly the same values as the
## fifth column ("Difference...") of Table H.6.
H.3$bk - predict(LM)
## Uncertainty in a predicted value. This will give
## exactly the values in H.3.4.
predict(LM, newdata = data.frame(tk = 30), se.fit = TRUE)
######### GUM 2008 Supplement 1 (1) #######################
## Example from 9.2.2 of the GUM 2008 Supplement 1
## (see 'References'), normally distributed input
## quantities. Assign values as in 9.2.2.1.
EXPR7 <- expression(X1 + X2 + X3 + X4)
X1 <- c(0, 1)
X2 <- c(0, 1)
X3 <- c(0, 1)
X4 <- c(0, 1)
DF7 <- cbind(X1, X2, X3, X4)
RES7 <- propagate(expr = EXPR7, data = DF7, type = "stat",
do.sim = TRUE, verbose = TRUE, nsim = 1E5)
## This will give exactly the same results as in
## 9.2.2.6, Table 2.
RES7$prop
RES7$sim
######### GUM 2008 Supplement 1 (2) #######################
## Example from 9.3 of the GUM 2008 Supplement 1
## (see 'References'), mass calibration.
## Formula 24 in 9.3.1.3 and values as in 9.3.1.4, Table 5.
EXPR8 <- expression((Mrc + dMrc) * (1 + (Pa - Pa0) * ((1/Pw) - (1/Pr))) - Mnom)
Mrc <- rnorm(1E5, 100000, 0.050)
dMrc <- rnorm(1E5, 1.234, 0.020)
Pa <- runif(1E5, 1.10, 1.30) ## E(Pa) = 1.2, (b-a)/2 = 0.1
Pw <- runif(1E5, 7000, 9000) ## E(Pw) = 8000, (b-a)/2 = 1000
Pr <- runif(1E5, 7950, 8050) ## E(Pr) = 8000, (b-a)/2 = 50
Pa0 <- 1.2
Mnom <- 100000
DF8 <- cbind(Mrc, dMrc, Pa, Pw, Pr, Pa0, Mnom)
RES8 <- propagate(expr = EXPR8, data = DF8, type = "sim",
do.sim = TRUE, verbose = TRUE, nsim = 1E5)
## This will give exactly the same results as in
## 9.3.2.3, Table 6
RES8$prop
RES8$sim
######### GUM 2008 Supplement 1 (3) #######################
## Example from 9.4 of the GUM 2008 Supplement 1
## (see 'References'), comparioson loss in microwave
## power meter calibration, zero covariance.
## Formula 28 in 9.4.1.5 and values as in 9.4.1.7.
EXPR9 <- expression(X1^2 - X2^2)
X1 <- c(0.050, 0.005)
X2 <- c(0, 0.005)
DF9 <- cbind(X1, X2)
RES9a <- propagate(expr = EXPR9, data = DF9, type = "stat",
do.sim = TRUE, verbose = TRUE, nsim = 1E5)
## This will give exactly the same results as in
## 9.4.2.2.7, Table 8, x1 = 0.050.
RES9a$prop
RES9a$sim
## Using covariance matrix with r(x1, x2) = 0.9
## We convert to covariances using cor2cov.
COR9 <- matrix(c(1, 0.9, 0.9, 1), nrow = 2)
COV9 <- cor2cov(COR9, c(0.005^2, 0.005^2))
colnames(COV9) <- c("X1", "X2")
rownames(COV9) <- c("X1", "X2")
RES9b <- propagate(expr = EXPR9, data = DF9, type = "stat",
use.cov = COV9, do.sim = TRUE,
verbose = TRUE, nsim = 1E5)
## This will give exactly the same results as in
## 9.4.3.2.1, Table 9, x1 = 0.050.
RES9b$prop
RES9b$sim
######### GUM 2008 Supplement 1 (4) #######################
## Example from 9.5 of the GUM 2008 Supplement 1
## (see 'References'), gauge block calibration.
## Assignment of PDF's as in Table 10 of 9.5.2.1.
EXPR10 <- expression(Ls + D + d1 + d2 - Ls *(da *(t0 + Delta) + as * dt) - Lnom)
Lnom <- 50000000
Ls <- propagate:::rst(100000, mean = 50000623, sd = 25, df = 18)
D <- propagate:::rst(100000, mean = 215, sd = 6, df = 25)
d1 <- propagate:::rst(100000, mean = 0, sd = 4, df = 5)
d2 <- propagate:::rst(100000, mean = 0, sd = 7, df = 8)
as <- runif(100000, 9.5E-6, 13.5E-6)
t0 <- rnorm(100000, -0.1, 0.2)
Delta <- propagate:::rarcsin(100000, -0.5, 0.5)
da <- propagate:::rctrap(100000, -1E-6, 1E-6, 0.1E-6)
dt <- propagate:::rctrap(100000, -0.050, 0.050, 0.025)
DF10 <- cbind(Ls, D, d1, d2, as, t0, Delta, da, dt, Lnom)
RES10 <- propagate(expr = EXPR10, data = DF10, type = "sim",
use.cov = FALSE, verbose = TRUE,
alpha = 0.01)
RES10
## This gives the same results as in 9.5.4.2, Table 11.
## However: results are exacter than in the GUM 2008
## manual, especially when comparing sd(Monte Carlo)
## sd(second-order Taylor expansion)!
## GUM 2008 gives 32 and 36, respectively.
RES10$prop
RES10$sim
########## Comparison to Pythons 'soerp' ###################
## Exactly the same results as under
## https://pypi.python.org/pypi/soerp !
EXPR11 <- expression(C * sqrt((520 * H * P)/(M *(t + 460))))
H <- c(64, 0.5)
M <- c(16, 0.1)
P <- c(361, 2)
t <- c(165, 0.5)
C <- c(38.4, 0)
DAT11 <- makeDat(EXPR11)
RES11 <- propagate(expr = EXPR11, data = DAT11, type = "stat",
do.sim = TRUE, verbose = TRUE)
RES11
# }
Run the code above in your browser using DataLab