## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4.5, out.width = "90%")
library(propagate)
M <- 1E5          # Monte Carlo trials; use 1E6 (or adaptive = TRUE) for final results

## ----typeA--------------------------------------------------------------------
I.obs <- c(4.358, 4.349, 4.353, 4.361, 4.347, 4.356, 4.351, 4.359, 4.352, 4.355)
c(mean = mean(I.obs), u = sd(I.obs)/sqrt(length(I.obs)))

## ----budget-------------------------------------------------------------------
BUD <- data.frame(name = c("V", "I", "dI", "phi"),
                  mean = c(230.1, mean(I.obs), 0, 0),
                  u    = c(NA, sd(I.obs)/sqrt(10), NA, NA),
                  U    = c(0.4, NA, NA, NA),
                  k    = c(2, NA, NA, NA),
                  a    = c(NA, NA, 0.01, 0.12),
                  df   = c(Inf, 9, Inf, Inf),
                  dist = c("normal", "t", "rect", "rect"))
BUD

## ----makeDat------------------------------------------------------------------
set.seed(123)
EXPR <- expression(P = V * (I + dI) * cos(phi))
DAT  <- makeDat(EXPR, data = BUD, layout = "budget", output = "sim", nsim = M)
head(DAT, 3)
attr(DAT, "budget")[, c("name", "mean", "u", "df", "dist")]

## ----propagate----------------------------------------------------------------
PROP <- propagate(EXPR, DAT, type = "sim", check = FALSE)

## ----results------------------------------------------------------------------
PROP$prop
PROP$sim
PROP$ci.short

## ----shift--------------------------------------------------------------------
c(GUF.2nd = unname(PROP$prop["mean.2"] - PROP$prop["mean.1"]),
  MCM     = unname(PROP$sim["mean.MC"] - PROP$prop["mean.1"]),
  approx  = unname(-PROP$prop["mean.1"] * 0.12^2/6))

## ----summary------------------------------------------------------------------
summary(PROP)

## ----sobol--------------------------------------------------------------------
S <- sobol(PROP)
round(rbind(first = S$S, total = S$ST), 3)

## ----sobolcheck---------------------------------------------------------------
P0 <- unname(PROP$prop["mean.1"]); a <- 0.12
Vphi <- P0^2 * a^4/45                         # exact for a rectangular phase angle
c(share.phi = Vphi/(unname(PROP$prop["u.1"])^2 + Vphi), MC.variance = unname(PROP$sim["u.MC"])^2,
  first.order.plus.phi = unname(PROP$prop["u.1"])^2 + Vphi)

## ----plot---------------------------------------------------------------------
plot(PROP)

## ----fitDistr-----------------------------------------------------------------
FIT <- fitDistr(PROP, nmax = 5000, verbose = FALSE)
head(FIT$stat[, c("Distribution", "k", "BIC", "dBIC", "weight")], 5)

## ----conformityMC-------------------------------------------------------------
CMC <- conformity(PROP, lower = 994, upper = 1006, rule = "all")
CMC

## ----conformityPlot, fig.height = 3.5-----------------------------------------
plot(CMC)

## ----conformityGUF------------------------------------------------------------
CGUF <- conformity(PROP, lower = 994, upper = 1006, method = "GUF", rule = c("stringent", "simple"))
CGUF

## ----report-------------------------------------------------------------------
f <- file.path(tempdir(), "heater-power.pdf")
REP <- report(PROP, file = f, title = "Power of a 1000 W heating element",
              measurand = "P", unit = "W")
REP$nonlinearity

## ----nls----------------------------------------------------------------------
DNase1 <- subset(DNase, Run == 1)
fm <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = DNase1,
          start = list(Asym = 3, xmid = 0, scal = 1))

## ----predictNLS1, results = "hide", message = FALSE---------------------------
set.seed(123)
PR1 <- predictNLS(fm, newdata = data.frame(conc = 2), nsim = M)
PR2 <- predictNLS(fm, newdata = data.frame(conc = 2), newerror = data.frame(conc = 0.2), nsim = M)

## ----predictNLS2--------------------------------------------------------------
round(rbind("u(conc) = 0" = PR1$summary[1, ], "u(conc) = 0.2" = PR2$summary[1, ]), 4)

## ----predictNLS3, results = "hide", message = FALSE---------------------------
CONC <- exp(seq(log(0.05), log(12.5), length.out = 25))
CB <- predictNLS(fm, newdata = data.frame(conc = CONC), nsim = M)$summary
PB <- predictNLS(fm, newdata = data.frame(conc = CONC), interval = "prediction", nsim = M)$summary

## ----predictNLS4--------------------------------------------------------------
plot(density ~ conc, data = DNase1, log = "x", pch = 16, xlab = "concentration", ylab = "optical density")
lines(CONC, CB[, "mean.1"], lwd = 2)
matlines(CONC, CB[, 11:12], lty = 2, col = 2)      # Monte Carlo confidence band
matlines(CONC, PB[, 11:12], lty = 3, col = 4)      # Monte Carlo prediction band
legend("topleft", c("fit", "95 % confidence band", "95 % prediction band"), lty = 1:3, col = c(1, 2, 4), bty = "n")

## ----MV-----------------------------------------------------------------------
H6 <- cbind(V   = c(5.007, 4.994, 5.005, 4.990, 4.999, 4.999),
            I   = c(19.663, 19.639, 19.640, 19.685, 19.678, 19.661)/1000,
            phi = c(1.0456, 1.0438, 1.0468, 1.0428, 1.0433, 1.0445))
EXMV <- expression(R = (V/I) * cos(phi), X = (V/I) * sin(phi), Z = V/I)
set.seed(123)
MV <- propagateMV(EXMV, H6, type = "raw", cov = "full", mc.margins = "mvt", nsim = M)
MV

## ----MVplot, fig.height = 6.5-------------------------------------------------
plot(MV)

## ----doGauge, eval = FALSE----------------------------------------------------
# G <- doGauge(c("JCGM100", "JCGM101", "JCGM102"), output = "html", outdir = tempdir())
# G$summary                     # PASS / NOTE / DIFF per demo
# browseURL(G$files[1])         # the full comparison of the first demo