knitr::opts_chunk$set(echo=T, comment=NA, error=T, warning=F, message = F, fig.align = 'center')
z.HY_derivative <- function(pres.pr, temp.pr, verbose = FALSE) { # Hall-Yarborough correlation modified to use the Newton-Raphson method tol <- 1E-13 f <- function(y) { - A * pres.pr + (y + y^2 + y^3 - y^4) / (1 - y)^3 - B * y^2 + C * y^D } fdot <- function(y) { (1 + 4 * y + 4 * y^2 - 4 * y^3 + y^4 ) / (1 - y)^4 - 2 * B * y + C * D * y^(D-1) } t <- 1 / temp.pr A <- 0.06125 * t * exp(-1.2 * (1 - t)^2) B <- t * (14.76 - 9.76 * t + 4.58 * t^2) C <- t * (90.7 - 242.2 * t + 42.4 * t^2) D <- 2.18 + 2.82 * t # first guess for y yk <- 0.0125 * pres.pr * t * exp(-1.2 * (1 - t)^2) delta <- 1 i <- 1 # itertations while (delta >= tol) { fyk <- f(yk) if (abs(fyk) < tol) break yk1 <- yk - f(yk) / fdot(yk) delta <- abs(yk - yk1) if (verbose) cat(sprintf("%3d %10f %10f %10f \n", i, delta, yk, fyk)) yk <- yk1 i <- i + 1 } if (verbose) cat("\n") y <- yk z <- A * pres.pr / y if (verbose) print(z) return(z) } z.HY_derivative(0.5, 1.3)
Now, we test with the the same test values in the Fatooreshi paper.
# test HY with 1st-derivative using the values from paper ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) tpr <- c(1.3, 1.5, 1.7, 2) tbl <- sapply(ppr, function(x) sapply(tpr, function(y) z.HY_derivative(pres.pr = x, temp.pr = y))) rownames(tbl) <- tpr colnames(tbl) <- ppr print(tbl) z_hy_deriv <- tbl library(ggplot2) plot(x = ppr, y = tbl[1,], type = "l", main = "z @ Tpr = 1.3", ylab = "z") plot(x = ppr, y = tbl[2,], type = "l", main = "z @ Tpr = 1.5", ylab = "z") plot(x = ppr, y = tbl[3,], type = "l", main = "z @ Tpr = 1.7", ylab = "z") plot(x = ppr, y = tbl[4,], type = "l", main = "z @ Tpr = 2.0", ylab = "z")
Values from Standing and Katz chart These values have been read from the SK-chart.
(0.5, 1.3) = 0.92 (1.5, 1.3) = 0.76 (2.5, 1.3) = 0.64 (3.5, 1.3) = 0.63 (4.5, 1.3) = 0.68 (5.5, 1.3) = 0.76 (6.5, 1.3) = 0.84 0.5, 1.5 = 0.94 1.5, 1.5 = 0.86 2.5, 1.5 = 0.79 3.5, 1.5 = 0.77 4.5, 1.5 = 0.79 5.5, 1.5 = 0.84 6.5, 1.5 = 0.89 (0.5, 1.7) = 0.97 (1.5, 1.7) = 0.92 (2.5, 1.7) = 0.87 (3.5, 1.7) = 0.86 (4.5, 1.7) = 0.865 (5.5, 1.7) = 0.895 (6.5, 1.7) = 0.94 (0.5, 2.0) = 0.985 (1.5, 2.0) = 0.957 (2.5, 2.0) = 0.941 (3.5, 2.0) = 0.938 (4.5, 2.0) = 0.945 (5.5, 2.0) = 0.97 (6.5, 2.0) = 1.01
ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) tpr <- c(1.3, 1.5, 1.7, 2) # data points from Standing-Katz chart z <- c( c(0.92, 0.76, 0.64, 0.63, 0.68, 0.76, 0.84), c(0.94, 0.86, 0.79, 0.77, 0.79, 0.84, 0.89), c(0.97, 0.92, 0.87, 0.86, 0.865, 0.895, 0.94), c(0.985, 0.957, 0.941, 0.938, 0.945, 0.97, 1.01) ) # create a matrix of SK chart points mx <- matrix(z, nrow = 4, ncol = 7, byrow = TRUE) rownames(mx) <- tpr colnames(mx) <- ppr print(mx) z_sk_chart <- mx save(z_hy_deriv, z_sk_chart, file = "../data/z_hy.rda") library(ggplot2) plot(x = ppr, y = mx[1,], type = "l", main = "z SK @ Tpr = 1.3", ylab = "z") lines(x = ppr, y = tbl[1,], col = "red") legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red")) plot(x = ppr, y = mx[2,], type = "l", main = "z SK @ Tpr = 1.5", ylab = "z") lines(x = ppr, y = tbl[2,], col = "red") legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red")) plot(x = ppr, y = mx[3,], type = "l", main = "z SK @ Tpr = 1.7", ylab = "z") lines(x = ppr, y = tbl[3,], col = "red") legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red")) plot(x = ppr, y = mx[4,], type = "l", main = "z SK @ Tpr = 2.0", ylab = "z") lines(x = ppr, y = tbl[4,], col = "red") legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red"))
We will compare the two tables:
# calculate the error in percentage (tbl - mx) / tbl * 100
HY is not recommended for $T_{pr} < 1$.
tbl mx
plot(x = tbl[1, ], y = mx[1, ])
x = tbl[1, ] y = mx[1, ] pvalues <- data.frame(x = x, y = y)
plot( # If this is a 2-column matrix, then first column pvalues[,1] will be by default the x-axis and second column pvalue[,2] the y-axis # Can be a matrix with 2 columns or a data.frame with 2 columns x = pvalues, # x = tbl[1, ], y = mx[1, ], # Analogous if you have two separated vectors instead of two columns, change to something like: # x = pvalues[,1], # x-axis observation values # y = pvalues[,2], # y-axis observation values # x- and y-axis ranges [0,1] # xlim=c(0, 1), # ylim=c(0, 1), xlim = c(0.6, 0.95), ylim = c(0.6, 0.95), # Select filled dots as the symbols pch=16, # Conditional color vector based on where the observation is located col=apply(pvalues, MARGIN=1, FUN = function(x) {ifelse(1*x[1] + 0.001 > x[2] & 1*x[1] - 0.001 < x[2], # Color for dots inside the area (semi-transparent black) "#00000075", # Color for dots outside the area (semi-transparent blue) "#0000FF75") } ), # Axis labels xlab="p values in differences condition 1", ylab="p values in differences condition 2" )
x <- ppr y1 <- tbl[1,] y2 <- mx[1,] plot( seq_along(y1), y1, type="l", col="red" ) par(new=TRUE) plot( seq_along(y2), y2, type="l", col="green" )
qplot(y=tbl[1,],x=mx[1,]) qplot(y=tbl[1,],x=mx[1,]) + geom_jitter(width = 0.01, height = 0.01)
qplot(y=tbl[1,],x=tbl[1,], geom_point(y =mx[1,], x=mx[1,], col = "green"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.