R/uncertainty-internal.R

Defines functions .internal_mc .internal_kragten .internal_gso .internal_gfo .internal_dof_eff .internal_build_envir

.internal_build_envir <-
function(name, value) {
  n <- length(name)
  res <- paste(name[1], "=", value[1], sep = "")
  if (n > 1) {
    for (i in 2:n) {
      res <- paste(res, ", ", paste(name[i], "=", value[i], sep = ""), sep = "")
    }
  }
  res <- paste("list(", res, ")", sep = "")
  return(res)
}
.internal_dof_eff <-
function(c_u, dof, cor = diag(1, length(c_u))) {
	ct <- sqrt( t(c_u) %*% cor %*% c_u )
	res <- 0
	rr_1 <- sum(c_u ^ 4 / dof)
	rr_2 <- 0
	rr_3 <- 0
	n <- length(c_u)
	# this implements the Castrup.2010 algorithm
	for (i in 1:(n - 1)) {
		for (j in (i + 1):n) {
			rr_2 <- rr_2 + cor[i, j] ^ 2 * (c_u[i] ^ 2 / dof[i]) *
			  (c_u[j] ^ 2 / dof[j]) * (dof[i] + dof[j] + 1 / 2)
			rr_3 <- rr_3 + cor[i, j] * (c_u[i] ^ 3 / dof[i] * c_u[j] +
			                              c_u[i] * c_u[j] ^ 3 / dof[j])
		}
	}
	res <- ct ^ 4 / (rr_1 + rr_2 + 2 * rr_3)
	return (res)
}


.internal_gfo <-
function(name, label = name, model, ub, alpha = 0.05) {
	envir_str <- .internal_build_envir(ub$name, ub$mean)
	n <- length(ub$name)
	
	if (any(is.na(ub$cor))) {
	  ub$cor <- diag(1, n)
	}

	data_cov <- matrix(ub$u, n, n) * ub$cor * t(matrix(ub$u, n, n))

	# get the partial derivatives of the expression wrt each of the variables
	# compute the contribution to the uncertainty by evaluating the list of
	# variables as parameter of the environment.
	# this table contains the variables, the gradient of the function wrt the
	# variable and the uncertainty contribution for that variable.
	res <- data.frame(var = rep(NA, n), partial = rep(NA, n),
	                  coefficient = rep(NA, n), contribution = rep(NA, n))

	res$var <- ub$name
	for (i in 1:n) {
		res$partial[i] <- as.expression(D(parse(text = model), ub$name[i]))
		res$coefficient[i] <- eval(res$partial[i],
		                           envir = eval(parse(text = envir_str)))
	}

	res$contribution <- res$coefficient * ub$u

	# correlation contributions
	JJ <- res$coefficient %*% t(res$coefficient)

	GUM_mu <- eval(parse(text = model), envir = eval(parse(text = envir_str)))
	GUM_u <- sqrt(sum(JJ * data_cov))

	cor_contribution <- GUM_u ^ 2 - sum(res$contribution ^ 2)

	GUM_dof <- .internal_dof_eff(c_u = res$contribution, dof = ub$dof,
	                             cor = ub$cor)

#	print(GUM_dof)
	
	GUM_U <- qt(1 - alpha / 2, GUM_dof) * GUM_u

	su <- NA
	if (GUM_dof > 2) {
		su <- GUM_u * sqrt(GUM_dof / (GUM_dof - 2))
	}

	result_GUM <- list(method = "GFO",
		call = match.call(),
		measurand_name = name,
		measurand_model = model,
		measurand_label = label,
		mean = GUM_mu,
		sd = GUM_u,
		u = GUM_u,
		alpha = alpha,
		dof = GUM_dof,
		U = GUM_U,
		lcl = GUM_mu - GUM_U,
		ucl = GUM_mu + GUM_U,

		variables = list(name = ub$name, label = ub$label, mean = ub$mean, u = ub$u,
		                 dof = ub$dof, coeff = res$coefficient, cor = ub$cor),
		contribution = res$contribution,
		cor_contribution = cor_contribution,

		# additional ouput
		partial = res$partial,
		coeff = res$coefficient,
		su = su)

	return(result_GUM)
}
.internal_gso <-
function(name, label = name, model, ub, alpha = 0.05) {
	envir_str <- .internal_build_envir(ub$name, ub$mean)
	n <- length(ub$name)

	if (any(is.na(ub$cor))) ub$cor <- diag(1, n)

	data_cov <- matrix(ub$u, n, n) * ub$cor * t(matrix(ub$u, n, n))

	# get the partial derivatives of the expression wrt each of the variables.
	# compute the contribution to the uncertainty by evaluating the list of
	# variables as parameter of the environment.
	# this table contains the variables, the gradient of the function wrt the
	# variable and the uncertainty contribution for that variable.
	res <- data.frame(var = rep(NA, n), partial = rep(NA, n),
	                  coefficient = rep(NA, n), contribution = rep(NA, n))

	res$var <- ub$name
	for (i in 1:n) {
		res$partial[i] <- as.expression(D(parse(text = model), ub$name[i]))
		res$coefficient[i] <- eval(res$partial[i],
		                           envir = eval(parse(text = envir_str)))
	}
	res$contribution <- res$coefficient * ub$u

	JJ <- res$coefficient %*%t (res$coefficient)
	GFO_u <- sqrt(sum(JJ * data_cov))

	J <- rep(NA, n)
	H <- matrix(NA, n, n)
	for (ii in 1:n) {
		J[ii] <- eval(res$partial[ii], envir = eval(parse(text = envir_str)))
		for (jj in 1:n) {
			H[ii, jj] <- eval(D(res$partial[ii], ub$name[jj]),
			                 envir = eval(parse(text = envir_str)))
		}
	}

	# second order expected value
	mu_so <- eval(parse(text = model), envir = eval(parse(text = envir_str))) +
	  sum(diag(H %*% data_cov)) / 2
	# second order variance
	u_so <- sqrt(GFO_u ^ 2 + sum(diag(H %*% data_cov %*% H %*% data_cov)) / 2)

	dof_so <- .internal_dof_eff(c_u = res$contribution, dof = ub$dof,
	                            cor = ub$cor)

	cor_contribution <- u_so ^ 2 - sum(res$contribution ^ 2)

	U_so <- qt(1 - alpha / 2, dof_so) * u_so

	res <- list(method = "GSO",
		call = match.call(),
		measurand_name = name,
		measurand_model = model,
		measurand_label = label,
		mean = mu_so,
		sd = u_so,
		u = u_so,
		alpha = alpha,
		dof = dof_so,
		U = U_so,
		lcl = mu_so - U_so,
		ucl = mu_so + U_so,

		variables = list(name = ub$name, label = ub$label, mean = ub$mean, u = ub$u,
		                 dof = ub$dof, coeff = res$coefficient, cor = ub$cor),
		contribution = res$contribution,
		cor_contribution = cor_contribution,

		# additional ouput
		partial = res$partial,
		coeff = res$coefficient,
		su = u_so * sqrt(dof_so / (dof_so - 2)))

	return(res)
}
.internal_kragten <-
function(name, label = name, model, ub, alpha = 0.05) {
	envir_str <- .internal_build_envir(ub$name, ub$mean)
	n <- length(ub$name)

	f <- rep(NA, n + 1)

	f[n + 1] <- eval(parse(text = model), envir = eval(parse(text = envir_str)))

	for (i in 1:n) {
		# build a new point at (x[1],...,x[i]+u[i],...,x[n]) to evaluate the model
		new_mean <- ub$mean
		new_mean[i] <- new_mean[i] + ub$u[i]
		new_envir_str <- .internal_build_envir(ub$name, new_mean)
		f[i] <- eval(parse(text = model), envir = eval(parse(text = new_envir_str)))
	}

	Kragten_mean <- f[n + 1]
	Kragten_contribution <- (f[ - (n + 1)] - Kragten_mean)
	Kragten_u <- sqrt(sum(Kragten_contribution %*% ub$cor %*%
	                        Kragten_contribution))
	Kragten_cor_contribution <- Kragten_u ^ 2 - sum(Kragten_contribution ^ 2)
	Kragten_dof <- .internal_dof_eff(c_u = Kragten_contribution, dof = ub$dof,
	                                 cor = ub$cor)
	Kragten_U <- qt(1 - alpha / 2, Kragten_dof) * Kragten_u
	Kragten_coeff <- Kragten_contribution
	Kragten_coeff[ub$u > 0] <- Kragten_contribution[ub$u > 0] / ub$u[ub$u > 0]
	Kragten_coeff[ub$u == 0] <- NA

	result_Kragten <- list(method = "Kragten",
		call = match.call(),
		measurand_name = name,
		measurand_model = model,
		measurand_label = label,
		mean = Kragten_mean,
		sd = Kragten_u,
		u = Kragten_u,
		alpha = alpha,
		dof = Kragten_dof,
		U = Kragten_U,
		lcl = Kragten_mean - Kragten_U,
		ucl = Kragten_mean + Kragten_U,

		variables = list(name = ub$name, label = ub$label, mean = ub$mean, u = ub$u,
		                 dof = ub$dof, coeff = Kragten_coeff, cor = ub$cor),
		contribution = Kragten_contribution,
		cor_contribution = Kragten_cor_contribution,

		su = Kragten_u * sqrt(Kragten_dof / (Kragten_dof - 2))

		# no additional output
		)

	return(result_Kragten)
}
.internal_mc <-
function(name, label = name, model, ub, alpha = 0.05, B = 2000) {

	n <- length(ub$name)

	rvalues <- matrix(NA, B, n)

	var_env <- new.env()

	if (sum(abs(ub$cor - diag(1,n))) < sqrt(.Machine$double.eps)) {
	  # variables are uncorrelated
		# get the simulated variables
		for (i in 1:n) {
			if (ub$distribution[i] == "bernoulli") {
				rvalues[, i]<- rbinom(B, size = 1, ub$mean[i])
			} else if (ub$distribution[i] == "beta") {
				a <- ub$mean[i] * (ub$mean[i] * (1 - ub$mean[i]) - ub$u[i] ^ 2) /
				  ub$u[i] ^ 2
				b <- (1 - ub$mean[i]) * (ub$mean[i] * (1 - ub$mean[i]) - ub$u[i] ^ 2) /
				  ub$u[i] ^ 2
				rvalues[,i] <- rbeta(B, a, b)
			} else if (ub$distribution[i] == "binomial") {
				rvalues[, i] <- rbinom(B, size = ub$mean[i], ub$u[i])
			} else if (ub$distribution[i] == "cauchy") {
				rvalues[, i] <- rcauchy(B, location = ub$mean[i], scale = ub$u[i])
			} else if (ub$distribution[i] == "chisq") {
				rvalues[, i] <- rchisq(B, ub$dof[i])
			} else if (ub$distribution[i] == "exp") {
				rvalues[, i] <- rexp(B, rate = 1 / ub$mean[i])
			} else if (ub$distribution[i] == "f") {
				# this uses two dof as parameters
				rvalues[, i] <- rf(B, ub$mean[i], ub$u[i])
			} else if (ub$distribution[i] == "gamma") {
				rvalues[, i] <- rgamma(B, shape = (ub$mean[i] / ub$u[i]) ^ 2,
				                       scale = ub$u[i] ^ 2 / ub$mean[i])
			} else if (ub$distribution[i] == "lognormal") {
				rvalues[, i] <- exp(rnorm(B, log(ub$mean[i]) - 1 / 2 *
				                            log(1 + ub$u[i] ^ 2 / ub$mean[i] ^ 2),
				                          sqrt(log(1 + ub$u[i] ^ 2 / ub$mean[i] ^ 2))))
			} else if (ub$distribution[i] == "poisson") {
				rvalues[, i] <- rpois(B, ub$mean[i])
			} else if (ub$distribution[i] == "normal") {
				rvalues[, i] <- rnorm(B, ub$mean[i], ub$u[i])
			} else if (ub$distribution[i] == "unif") {
				rvalues[, i] <- runif(B, ub$mean[i] - sqrt(3) * ub$u[i],
				                      ub$mean[i] + sqrt(3) * ub$u[i])
			} else if (ub$distribution[i] == "t") {
				rvalues[, i] <- ub$mean[i] + rt(B, ub$dof[i]) * ub$u[i] /
				  sqrt(ub$dof[i] / (ub$dof[i] - 2))
			} else if (ub$distribution[i] == "triangular") {
				rvalues[, i] <- rtriangle(B, ub$mean[i] - sqrt(6) * ub$u[i],
				                          ub$mean[i] + sqrt(6) * ub$u[i], ub$mean[i])
			} else if (ub$distribution[i] == "weibull") {
				rvalues[, i] <- rweibull(B, shape = ub$mean[i], scale = ub$u[i])
			} else if (ub$distribution[i] == "arcsine") {
				rvalues[, i] <- ub$mean[i] + ub$u[i] * sin(2 * pi * runif(B))
			} else {
				rvalues[, i] <- rnorm(B, ub$mean[i], ub$u[i])
			}
			eval(parse(text = paste(ub$name[i]," <- rvalues[,", i, "]", sep = "")),
			     envir = var_env)
		}
	}  else {
	  # variables are correlated, generate approximate random variables
		z <- rmvnorm(B, mean = rep(0, n), sigma = ub$cor)
		z <- pnorm(z)

		for (i in 1:n) {
			if (ub$distribution[i] == "bernoulli") {
				rvalues[, i] <- qbinom(z[, i], size = 1, ub$mean[i])
			} else if (ub$distribution[i] == "beta") {
				a <- ub$mean[i] * (ub$mean[i] * (1 - ub$mean[i]) - ub$u[i] ^ 2) /
				  ub$u[i] ^ 2
				b <- (1 - ub$mean[i]) * (ub$mean[i] * (1 - ub$mean[i]) - ub$u[i] ^ 2) /
				  ub$u[i] ^ 2
				rvalues[, i] <- qbeta(z[, i], a, b)
			} else if (ub$distribution[i] == "binomial") {
				rvalues[, i] <- qbinom(z[, i], size = ub$mean[i], ub$u[i])
			} else if (ub$distribution[i] == "cauchy") {
				rvalues[, i] <- qcauchy(z[, i], location = ub$mean[i], scale = ub$u[i])
			} else if (ub$distribution[i] == "chisq") {
				rvalues[, i] <- qchisq(z[, i], ub$dof[i])
			} else if (ub$distribution[i] == "exp") {
				rvalues[, i] <- qexp(z[, i], rate = 1 / ub$mean[i])
			} else if (ub$distribution[i] == "f") {
				rvalues[, i] <- qf(z[, i], ub$mean[i], ub$u[i])
			} else if (ub$distribution[i] == "gamma") {
				rvalues[, i] <- qgamma(z[, i], shape = (ub$mean[i] / ub$u[i]) ^ 2,
				                       scale = ub$u[i] ^ 2 / ub$mean[i])
			} else if (ub$distribution[i] == "lognormal") {
				rvalues[, i] <- exp(qnorm(z[, i], log(ub$mean[i]) - 1 / 2 *
				                            log(1 + ub$u[i] ^ 2 / ub$mean[i] ^ 2),
				                          sqrt(log(1 + ub$u[i] ^ 2 / ub$mean[i] ^ 2))))
			} else if (ub$distribution[i] == "poisson") {
				rvalues[, i] <- qpois(z[, i], ub$mean[i])
			} else if (ub$distribution[i] == "normal") {
				rvalues[, i] <- qnorm(z[, i], ub$mean[i], ub$u[i])
			} else if (ub$distribution[i] == "unif") {
				rvalues[, i] <- qunif(z[, i], ub$mean[i] - sqrt(3) * ub$u[i],
				                      ub$mean[i] + sqrt(3) * ub$u[i])
			} else if (ub$distribution[i] == "t") {
				rvalues[, i] <- ub$mean[i]+qt(z[, i], ub$dof[i]) * ub$u[i] /
				  sqrt(ub$dof[i] / (ub$dof[i] - 2))
			} else if (ub$distribution[i] == "triangular") {
				rvalues[, i] <- qtriangle(z[, i], ub$mean[i] - sqrt(6) * ub$u[i],
				                          ub$mean[i] + sqrt(6) * ub$u[i], ub$mean[i])
			} else if (ub$distribution[i] == "weibull") {
				rvalues[, i] <- qweibull(z[, i], shape = ub$mean[i], scale = ub$u[i])
			} else if (ub$distribution[i] == "arcsine") {
				rvalues[, i] <- ub$mean[i] + ub$u[i] * sin(2 * pi * z[,i])
			} else {
				rvalues[, i] <- qnorm(z[, i], ub$mean[i], ub$u[i])
			}
			eval(parse(text = paste(ub$name[i], " <- rvalues[,", i, "]", sep = "")),
			     envir = var_env)
		}
	}

	f_mc <- eval(parse(text = model), envir = var_env)

	# infinite values are handled as NA
	f_mc[f_mc == Inf] <- NA
	f_mc[f_mc == -Inf] <- NA

	res_mean <- mean(f_mc, na.rm = TRUE)
	res_u <- sd(f_mc, na.rm = TRUE)

	######################################
	## uncertainty contribution estimation
	######################################
	res_contribution <- rep(1, n)
	# backup the simulated inputs with extension ".rv"
	for (i in 1:n) {
		eval(parse(text = paste(ub$name[i], ".rv", "<- ", ub$name[i], sep = "")),
		     envir = var_env)
	}

	# for each variable compute its contribution
	for (i in 1:n) {
		# restore the scalar values to the environment variables
		for (j in 1:n) {
			eval(parse(text = paste(ub$name[j], "<- ", ub$mean[j])), envir = var_env)
		}
		# now replace the ith variable with the simulated values into the environment
		eval(parse(text = paste(ub$name[i], "<- ", ub$name[i], ".rv", sep = "")),
		     envir = var_env)
		# compute the function while changing the ith variable
		f_i <- eval(parse(text = model), envir = var_env)
		# compute the sd of the values as the contribution for that variable
		res_contribution[i] <- sd(f_i - res_mean, na.rm = TRUE)
		if (is.na(res_contribution[i])) res_contribution[i] <- 0
	}

	res_cor_contribution <- res_u^2 - sum(res_contribution^2)
	res_coeff <- res_contribution
	res_coeff[ub$u > 0] <- res_contribution[ub$u > 0] / ub$u[ub$u > 0]
	res_coeff[ub$u == 0] <- NA

	## estimate the dof for the estimated uncertainty
	## update this code with the new criteria
	bb <- floor(sqrt(B))
	vi <- rep(NA, bb)
	for (j in 1:bb) {
		vi[j] <- var(f_mc[ ( (j - 1) * bb + 1):(j * bb)], na.rm = TRUE)
	}
	res_dof_ws <- 2 * res_u^4 / var(vi)
	res_U <- (f_mc[order(f_mc)[round(B * (1 - alpha / 2))]] -
	            f_mc[order(f_mc)[round(B * (alpha / 2))]]) / 2

	result_MC <- list(method = "MC",
		call = match.call(),
		measurand_name = name,
		measurand_model = model,
		measurand_label = label,
		mean = res_mean,
		sd = res_u,
		u = res_u,
		alpha = alpha,
		dof = res_dof_ws,
		U = res_U,
		lcl = f_mc[order(f_mc)[round(B * (alpha / 2))]],
		ucl = f_mc[order(f_mc)[round(B * (1 - alpha / 2))]],

		variables = list(name = ub$name, label = ub$label, mean = ub$mean, u = ub$u,
		                 dof = ub$dof, coeff = res_coeff, cor = ub$cor),
		contribution = res_contribution,
		cor_contribution = res_cor_contribution,

		# additional ouput
		ncycles = B,
		na_count = sum(is.na(f_mc)),
		simulated_values = f_mc,
		simulated_data = rvalues)

	return(result_MC)
}

Try the uncertainty package in your browser

Any scripts or data that you put into this service are public.

uncertainty documentation built on May 2, 2019, 7:03 a.m.