## Internal functions
iunifpdf <- function(x, y, sigma, assume = "Bayes"){
a <- (1 - sigma)
b <- (1 + sigma)
## Bayes Law inversion of uniform, a la Sethi
if(assume == "Bayes")
out <- (x * log(b / a)) ^ -1
## Inverse Uniform Distribution, a la Springborn
else
out <- y * (x ^ -2) / (b - a)
out[x < y / b] <- 0
out[x > y / a] <- 0
out
}
iunifcdf <- function(x, y, sigma, assume = "Bayes"){
a <- (1 - sigma)
b <- (1 + sigma)
## Bayes Law inversion of uniform, a la Sethi:
if(assume == "Bayes")
out <- log(x * b / y) / log(b / a)
## Inverse uniform distribution CDF, a la Springborn
else
out <- (b - y / x ) / (b - a)
out[x < y / b] = 0
out[x > y / a] = 1
out
}
ilognpdf <- function(x, y, sigma, assume = "Bayes"){
if(assume == "Bayes"){
# out <- exp(- (log(y) - log(x)) ^ 2 / (2 * sigma^2) ) / (x * sigma * sqrt(2 * pi))
out <- dlnorm(x, y, sigma)
} else {
out <- (y / (x * sigma * sqrt(2 * pi))) * exp(-(log(x) - y)^2 / (2 * sigma^2))
}
out[is.nan(out)] <- 0
out
}
ilogncdf <- function(x, y, sigma, assume = "Bayes"){
if(assume == "Bayes"){
N <- plnorm(x, y, sigma)
} else {
suppressWarnings(N <- qlnorm(x, y, sigma))
N[is.na(N)] <- 1
}
N
}
snap_to_grid <- function(x, grid) sapply(x, function(x) grid[which.min(abs(grid - x))])
pdfn <- function(p, mu, s, grid, pdf){
if(mu <= 0){
out <- as.integer(p == 0)
} else if(s == 0){ ## delta spike if s = 0
out <- as.numeric(as.integer(snap_to_grid(p, grid) == snap_to_grid(mu, grid)))
} else if(s > 0){ ## Evaluate pdf only for mu, s > 0
out <- pdf(p, mu, s)
} else { # all other cases. perhaps should be warning/error instead
out <- 0
}
out
}
pdf_matrix <- function(i_grid, j_grid, sigma, pdf, cdf){
n_i <- length(i_grid)
n_j <- length(j_grid)
out <- array(0, c(n_i, n_j))
A <- array(0, c(n_i, n_j))
for(i in 1:n_i){
## compute pdf on extended grid as loop
A[i, ] <- pdfn(j_grid, i_grid[i], sigma, j_grid, pdf)
## handle any all-zero rows
if(sum(A[i,]) == 0 || i_grid[i] == 0){
A[i, ] <- c(1, numeric(dim(A)[2]-1))
} else if(sigma > 0){
## normalize row
N <- cdf(j_grid[n_j], i_grid[i], sigma)
A[i, ] <- A[i, ] %*% t(N) / sum(A[i,])
## pile on boundary any density from extended grid
A[i,n_j] <- 1 - N + A[i, n_j]
}
}
out = A
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.