#' RSM predicted ROC-rating pdf for diseased cases
#' @param z The z-vector at which to evaluate the pdf.
#' @param mu The scalar RSM mu parameter.
#' @param lambda The scalar RSM lambda parameter.
#' @param nu The scalar RSM nu parameter.
#' @param lesDistr The lesion distribution 1D vector.
#'
#' @return pdf, density function for diseased cases
#'
#' @examples
#' RSM_pdfD(c(1,2),1,1,0.9, c(0.5, 0.5))
#' RSM_pdfD(c(1,2),1,1,0.5, c(0.2, 0.3, 0.5))
#'
#' @export
# this was the original form, simplified somewhat but otherwise
# identical to the 2017 version
RSM_pdfD <- function(z, mu, lambda, nu, lesDistr){
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
pdf <- 0
for (L in 1:length(lesDistr)){
a <- 1-nu/2+nu/2*erf_vect_cpp((z-mu)/sqrt(2))
b <- exp((-lambda/2)+lambda/2*erf_vect_cpp(z/sqrt(2)))
pdf <- pdf +
# This is identical to the pdfD2 form below obtained using Maple, only the factoring is different
lesDistr[L]*((a^(L-1)*b)*(L*nu/sqrt(2*pi)*exp(-(z-mu)^2/2)+a*lambda/sqrt(2*pi)*exp(-z^2/2)))
}
# the following two checks are included to test out another simplification of the Maple generated
# formulas, as well as one using calculus to get at the final expression (derivative of 1-TPF wrt
# z)
# the max(pdf ...) is needed as this function works with an input vector for z; but this causes occasional failures
# replaced max(pdf ...) with mean(pdf ...); median(pdf ...) might be even more resistant to outliers
# 4/2/21 changed stop criterion from 1e-16 to 1e-15
if (abs(mean(pdf - pdfD2(z, mu, lambda, nu, lesDistr))) > 1e-15) stop("Two forms disagree A")
if (abs(mean(pdf - pdfD3(z, mu, lambda, nu, lesDistr))) > 1e-15) stop("Two forms disagree B")
return (pdf)
}
#' RSM predicted ROC-rating pdf for non-diseased cases
#' @param z The z-vector at which to evaluate the pdf.
#' @param lambda The scalar RSM lambda parameter.
#'
#' @return pdf, density function for non-diseased cases
#'
#' @examples
#' RSM_pdfN(c(1,2),1)
#'
#' @export
#'
#'
RSM_pdfN <- function(z, lambda){
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
# verified using RSM_pdfN(1, 1) and (xROC_cpp(1, 1) - xROC_cpp(1 + 1e-8, 1)) / 1e-8
# following expression is identical to book equation 17.21
return(lambda * exp(-z^2/2) * exp(-lambda/2 * (1 - erf_vect_cpp(z/sqrt(2)))) / sqrt(2 * pi))
}
################################################################################
################################################################################
#' RSM predicted FROC abscissa
#' @param z The z-vector at which to evaluate the FROC abscissa.
#' @param lambda The scalar RSM lambda parameter.
#'
#' @return NLF, the abscissa of the FROC curve
#'
#' @examples
#' RSM_NLF(c(1,2),1)
#'
#' @export
#'
#'
RSM_NLF <- function(z, lambda){
# returns NLF, the abscissa of FROC curve
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
NLF <- lambda * (1 - pnorm(z))
return(NLF)
}
################################################################################
################################################################################
#' RSM predicted FROC ordinate
#' @param z The z-vector value at which to evaluate the FROC ordinate.
#' @param mu The scalar RSM mu parameter.
#' @param nu The scalar RSM nu prime parameter.
#'
#' @return LLF, the ordinate of the FROC curve
#'
#' @examples
#' RSM_LLF(c(1,2),1,0.5)
#'
#' @export
#'
#'
RSM_LLF <- function(z, mu, nu){
# bug fix 12/26/21
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
# returns LLF, the ordinate of FROC, AFROC curve
LLF <- nu * (1 - pnorm(z - mu))
return(LLF)
}
################################################################################
################################################################################
integrandFROC <- function(NLF, mu, lambda, nu){
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
zeta <- qnorm(1 - NLF / lambda)
LLF <- RSM_LLF(zeta, mu, nu)
return(LLF)
}
################################################################################
################################################################################
# y_AFROC_FPF is AFROC as a function of FPF + RSM parameters
# this function does not call any Cpp code
y_AFROC_FPF <- function(FPF, mu, lambda, nu){
# returns LLF, the ordinate of AFROC curve; takes FPF as the variable.
# AUC is calculated by integrating this function wrt FPF
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
tmp <- 1 / lambda * log(1 - FPF) + 1
tmp[tmp < 0] <- pnorm(-20)
zeta <- qnorm(tmp)
LLF <- RSM_LLF(zeta, mu, nu)
return(LLF)
}
################################################################################
################################################################################
# dpc 01/05/22 these comments were added while converting code to formula for chapter 21-optim-op-point
# needed formula for wAFROC ordinate, I know this is crazy, Einstein would never have done it this way :(
# finished see RJafrocFrocBook, search for rsm-pred-wafroc-curve 1/7/22
# checked from Console that following two give same results with following code
# UtilAnalyticalAucsRSM(mu = 2, lambda = 1, nu = 0.9, zeta1 = -3,
# lesDistr = c(0.1, 0.4, 0.4, 0.1), relWeights = c(0.2, 0.3, 0.1, 0.5))
# $aucROC
# [1] 0.9698827
# $aucAFROC
# [1] 0.8566404
# $aucwAFROC
# [1] 0.8448053
# following is R implementation
# see RJafrocFrocBook, search for rsm-pred-wafroc-curve 1/7/22
# see test_RSM-formulae.R
# contextStr <- "testing weights code with max 4 lesions per case: Cpp vs R"
# contextStr <- "testing weights code with max 4 lesions per case, random values: Cpp vs R"
# contextStr <- "testing weights code with max 10 lesions per case, random values: Cpp vs R"
# RSM_wLLF_R is wAFROC ordinate as a function of zeta + RSM parameters
# returns wLLF, the ordinate of wAFROC curve
# this has working C++ version named RSM_wLLF_cpp
# this is only here for me to check the C++ code
RSM_wLLF_R <- function(zeta, mu, nu, lesDistr, relWeights){
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
W <- UtilLesWghtsLD(lesDistr, relWeights)
f_L <- lesDistr
wLLF <- 0
for (L in 1:length(lesDistr)){
# outer loop sums over `lesionID`, abbreviated to L
wLLF_L <- 0
for (l_2 in 1:L){
# l_2 is the number of successes/hits
# inner loop sums over different numbers of hits,
# weighted by the probability of that many hits
# L is the trial size in dbinom
wLLF_L <- wLLF_L + W[L, l_2+1] * l_2 * dbinom(l_2, L, nu)
}
wLLF <- wLLF + f_L[L] * wLLF_L
}
wLLF <- wLLF * pnorm(mu - zeta)
return(wLLF)
}
#' RSM predicted wAFROC ordinate, cpp code
#' @param zeta The zeta-vector at which to evaluate the FROC ordinate.
#' @param mu The scalar RSM mu parameter.
#' @param nu The scalar RSM nu prime parameter.
#' @param lesDistr Lesion distribution vector.
#' @param relWeights The lesion weights matrix
#'
#' @return wLLF, the ordinate of the wAFROC curve
#'
#' @examples
#' RSM_wLLF(1, 1, 0.9, lesDistr = c(0.5, 0.4, 0.1), relWeights = c(0.7, 0.2, 0.1))
#' ## 0.34174
#'
#' @export
# This calls the C++ version named RSM_wLLF_cpp
RSM_wLLF <- function(zeta, mu, nu, lesDistr, relWeights){
return(RSM_wLLF_cpp(zeta, mu, nu, lesDistr, UtilLesWghtsLD(lesDistr, relWeights)))
}
################################################################################
################################################################################
#' RSM predicted ROC-abscissa as function of z
#' @param z The z-vector at which to evaluate the ROC-abscissa.
#' @param lambda The scalar RSM lambda parameter.
#'
#' @return FPF, the abscissa of the ROC
#'
#' @examples
#' RSM_FPF(c(-Inf,0.1,0.2,0.3),1)
#'
#' @export
RSM_FPF <- function(z, lambda) {
if (lambda < 0) stop("Incorrect value for lambda\n")
return(xROC_vect_cpp(z, lambda))
}
#' RSM predicted ROC-ordinate as function of z
#'
#' @param z The z-vector at which to evaluate the pdf.
#' @param mu The scalar RSM mu parameter.
#' @param lambda The scalar RSM lambda parameter.
#' @param nu The scalar nu parameter.
#' @param lesDistr The lesion distribution 1D vector.
#'
#' @return TPF, the ordinate of the ROC
#'
#' @examples
#' lesDistr <- c(0.1,0.3,0.6)
#' RSM_TPF(c(-Inf,0.1,0.2,0.3), 1, 1, 0.9, lesDistr)
#'
#' @export
RSM_TPF <- function(z, mu, lambda, nu, lesDistr) {
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
return(yROC_vect_cpp(z, mu, lambda, nu, lesDistr))
}
################################################################################
################################################################################
# y_wAFROC_FPF is wAFROC as a function of FPF + RSM parameters
y_wAFROC_FPF <- function(FPF, mu, lambda, nu, lesDistr, relWeights){
lesWghtDistr <- UtilLesWghtsLD(UtilLesDistr(lesDistr), relWeights)
# returns wLLF, the ordinate of AFROC curve; takes FPF as the variable.
# AUC is calculated by integrating this function wrt FPF
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
tmp <- 1 / lambda * log(1 - FPF) + 1
tmp[tmp < 0] <- pnorm(-20)
zeta <- qnorm(tmp)
# Cpp code
wLLF <- sapply(zeta, RSM_wLLF_cpp, mu = mu, nu = nu, lesDistr, lesWghtDistr)
return(wLLF)
}
# y_wAFROC_FPF_R is wAFROC ordinate as a function of FPF + RSM parameters
y_wAFROC_FPF_R <- function(FPF, mu, lambda, nu, lesDistr, W){
# returns wLLF, the ordinate of AFROC curve; takes FPF as the variable.
# AUC is calculated by integrating this function wrt FPF
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
tmp <- 1 / lambda * log(1 - FPF) + 1
tmp[tmp < 0] <- pnorm(-20)
zeta <- qnorm(tmp)
# R code
wLLF <- sapply(zeta, RSM_wLLF_R, mu = mu, nu = nu, lesDistr, W)
return(wLLF)
}
# error function
erf_R <- function (x) {
return(erf_vect_cpp(x)) # this implements the commented formula below
# return(2 * pnorm(sqrt(2) * x) - 1)
}
# added 12/22/20
# alternate form using Maple, Dec 2020
pdfD2 <- function(z, mu, lambda, nu, lesDistr){
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
pdf <- 0
for (L in 1:length(lesDistr)){
a <- 1-nu/2+nu/2*erf_vect_cpp((z-mu)/sqrt(2))
b <- exp((-lambda/2)+lambda/2*erf_vect_cpp(z/sqrt(2)))
pdf <- pdf +
lesDistr[L]*((a^(L-1))*L*nu*exp(-(z-mu)^2/2)+a^L*lambda*exp(-z^2/2))*b/sqrt(2*pi)
}
return (pdf)
}
# A is first term on rhs of book 17.22
A <- function(mu,nu,z,L)
{
# bug fix 12/26/21
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
return((1-nu/2+nu/2*erf_vect_cpp((z-mu)/sqrt(2)))^L)
}
# B is second term on rhs of book 17.22
B <- function(lambda,z)
{
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
return(exp(-lambda/2+lambda/2*erf_vect_cpp(z/sqrt(2))))
}
#dA is deriv. of A wrt z, Maple generated
dA <- function(mu,nu,z,L)
{
# bug fix 12/26/21
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
return((1-nu/2+nu/2*erf_vect_cpp((z-mu)/sqrt(2)))^(L-1)*L*nu*exp(-(z-mu)^2/2)/sqrt(2*pi))
}
#dB is deriv. of B wrt z, Maple generated
dB <- function(lambda,z)
{
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
return(lambda*exp(-z^2/2)*exp(-lambda/2+lambda/2*erf_vect_cpp(z/sqrt(2)))/sqrt(2*pi))
}
# added 12/22/20
# using Maple generated derivatives wrt z
pdfD3 <- function(z, mu, lambda, nu, lesDistr){
# bug fix 12/26/21
if (lambda < 0) stop("Incorrect value for lambda\n")
if (nu < 0) stop("Incorrect value for nu\n")
if (nu > 1) stop("Incorrect value for nu\n")
pdf <- 0
for (L in 1:length(lesDistr)){
# AB is the two terms in book 17.22
# dA is deriv. of A wrt z
# dB is deriv. of B wrt z
pdf <- pdf +
lesDistr[L]*(dA(mu,nu,z,L)*B(lambda,z)+A(mu,nu,z,L)*dB(lambda,z))
}
return (pdf)
}
################################################################################
# yROCVect_R <- function (zeta, mu, lambda, nu, lesDistr){
#
# TPF <- array(dim = length(zeta))
# for (il in 1:length(zeta)){
# for (i in 1:length(lesDistr)){
# TPF[il] = TPF[il] + lesDistr[i] *
# (1 - pow(1 - nu/2 + nu/2 * erf_cpp( (zeta[il] - mu) / sqrt(2.0) ) , (i+1)) * exp( (-lambda / 2) + 0.5 * lambda * erf_cpp(zeta[il] / sqrt(2.0))))
# }
# }
#
# return (TPF)
# }
# y_ROC_FPF_R <- function (FPF, mu, lambda, nu, lesDistr){
# for (il in 1:length(FPF)){
# temp = (1 / lambda) * log(1 - FPF[il]) + 1;
# if (temp <= 0){
# zeta[il] = -20;
# }else{
# zeta[il] = R::qnorm(temp, 0, 1, 1, 0);
# }
# }
#
# return (yROC_vect_cpp(zeta, mu, lambda, nu, lesDistr))
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.