Nothing
#####################################################################################
#
#
# FUNCTIONS THAT IMPLEMENT VARIOUS BAESIAN MODELS
#
#
#####################################################################################
# get.posterior(par, data, skel, knots.x, Gr, p.bkg)
#
# returns: wrapper to the target function.
# called from DEoptim
# arguments
# par: DEoptim parameter
# data: an object of type data
# skel: DEoptim parameter
# knots.x: the spline knot x-values.
# Gr: low-r G(r) information
# p.bkg: the probability that a single pixel contains "only" background.
get.posterior <- function(par, data, skel=NA, knots.x, Gr=NA, p.bkg=.5){
par <- relist(par, skel)
knots.y <- if(is.null(par$knots.y)) NA else par$knots.y
alpha <- if(is.null(par$alpha)) 1 else par$alpha
ADP <- if(is.null(par$ADP)) NA else par$ADP
knots.n <- length(knots.x)
pars <- if(is.null(par$pars)) NA else par$pars
if(!is.na(ADP)){
# recalculate coherent baseline
n.atoms <- data$fitADP$n.atoms
ADP <- rep(ADP, length(n.atoms))
scatter.length <- data$fitADP$scatter.length
N_total <- sum(n.atoms)
f.av2 <- (sum(n.atoms*scatter.length)/N_total)^2
f2.av <- sum(n.atoms*scatter.length^2)/N_total
L <- (f.av2-f2.av)/f.av2
expADP <- 0
for(i in 1:length(n.atoms))
expADP <- expADP + n.atoms[i]*scatter.length[i]^2*exp(-ADP[i]*data$x^2)/N_total/f2.av
data$SB <- 1-expADP*(1-L)
}
if(!is.na(pars[1]))
posterior <- logPosteriorAnalyt(data=data, pars=pars, p.bkg=p.bkg)
else{
posterior <- logPosterior(data=data, alpha=alpha, knots.x=knots.x,
knots.y=knots.y, Gr=Gr, p.bkg=p.bkg)
}
return(posterior)
}
####################################################################################
# logPosterior(data, alpha, knots.x, knots.y, Gr, p.bkg)
#
# returns: the optimization function from Fischer et al.: the negative log of
# the posterior.
# arguments
# data: an object of type data. See reference manual.
# alpha: scale parameter.
# knots.x: the spline knot x-values.
# knots.y: the spline knot y-values ('c' in Fischer et al.).
# Gr: low-r G(r) information
# p.bkg: the probability that a single pixel contains "only" background.
logPosterior <- function(data, alpha, knots.x, knots.y, Gr=NA, p.bkg=0.5) {
Phi <- basisMatrix(x=data$x, knots.x=knots.x)
# 1. Prior
psi.prior <- logPriorF(knots.x, knots.y)$prior
# 2. Likelihood:
psi.likelihood.bkg <- logLikelihoodBkg( knots.y=knots.y, y=data$y-data$SB, Phi=Phi, p.bkg=p.bkg, sigma=data$sigma)
psi.likelihood.signal <- logLikelihoodSignal(knots.y=knots.y, y=data$y-data$SB, Phi=Phi, p.bkg=p.bkg, sigma=data$sigma, lambda=data$lambda)
max.log <- pmax(psi.likelihood.bkg, psi.likelihood.signal) # parallel max --> vector of max values
# avoiding inf values:
psi.likelihood <- max.log + log(exp(psi.likelihood.bkg - max.log) + exp(psi.likelihood.signal - max.log))
psi.likelihood <- -sum(psi.likelihood)
# 3. Gr=-4*Pi*rho*r restriction
psi.gr <- 0
if(!is.na(Gr[1]))
psi.gr <- logProbabilityBkgR(y=data$y-data$SB, Gr=Gr, alpha=alpha, Phi=Phi, knots.y=knots.y)
# 4. Together
psi <- psi.prior + psi.likelihood + psi.gr
return(psi)
}
logPosteriorAnalyt <- function(data, pars, p.bkg=0.5) {
# 2. Likelihood:
psi.likelihood.bkg <- logLikelihoodBkgAnalyt(pars=pars, x=data$x, y=data$y-data$SB, p.bkg=p.bkg, sigma=data$sigma)
psi.likelihood.signal <- logLikelihoodSignalAnalyt(pars=pars, x=data$x, y=data$y-data$SB, p.bkg=p.bkg, sigma=data$sigma, lambda=data$lambda)
max.log <- pmax(psi.likelihood.bkg, psi.likelihood.signal) # parallel max --> vector of max values
# avoiding inf values:
psi.likelihood <- max.log + log(exp(psi.likelihood.bkg - max.log) + exp(psi.likelihood.signal - max.log))
psi.likelihood <- -sum(psi.likelihood)
return(psi.likelihood)
}
####################################################################################
# logProbabilityBkgR(y, knots.y, alpha, Phi, Gr)
#
# returns: the negative log of G(r)-part of the likelihood.
# arguments
# y: signal values.
# knots.y: the spline knot y-values ('c' in Fischer et al.).
# alpha: scale parameter.
# Phi: basis matrix.
# Gr: low-r G(r) information.
logProbabilityBkgR <- function(y, knots.y, alpha, Phi, Gr){
if(is.na(Gr$type1))
psi.gr.r1 <- 0
else if(Gr$type1=="gaussianNoise")
psi.gr.r1 <- logLikelihoodGrGauss(y=y, knots.y=knots.y, alpha=alpha, Phi=Phi, bkg.r=Gr$bkg.r,
sigma.r=Gr$sigma.r, matrix.FT=Gr$matrix.FT1)$likelihood
else if(Gr$type1=="correlatedNoise")
psi.gr.r1 <- logLikelihoodGrCorr(knots.y=knots.y, Phi=Phi, bkg.r=Gr$bkg.r,
KG.inv=Gr$KG.inv, matrix.FT=Gr$matrix.FT1)$likelihood
if(is.na(Gr$type2))
psi.gr.r2 <- 0
else if(Gr$type2=="secondDeriv")
psi.gr.r2 <- logPriorBkgRSmooth(bkg.r=Gr$matrix.FT2 %*% (Phi %*% t(t(knots.y))), D=Gr$D)$likelihood
else if(Gr$type2=="gaussianProcess")
psi.gr.r2 <- logPriorBkgRGP(bkg.r=Gr$matrix.FT2 %*% (Phi %*% t(t(knots.y))), covMatrix=Gr$covMatrix)$likelihood
return(psi.gr.r1 + psi.gr.r2)
}
########################################
# logPriorF(knots.x, knots.y)
#
# returns: the negative log of the prior (due to Fischer et al.)
# arguments:
# knots.x: the spline knot x-values ('xi' in Fischer et al.)
# knots.y: the spline knot y-values ('c' in Fischer et al.)
logPriorF <- function(knots.x, knots.y, Hessian=FALSE){
E <- length(knots.y)
if (length(knots.x) != E) {
stop("Inconsistent lengths of parameters knots.y and knots.x!")
}
D <- DMatrix(knots.x=knots.x)
detD <- D$det
D <- D$matrix
cDc <- as.vector(t(knots.y) %*% D %*% t(t(knots.y)))
# Prior: separate terms according to which variables contribute
prior.E <- (0.5 * E * log(pi) - log(gamma(0.5 * E))) - sum(log(1:E))
# "Depends" on knots.x, but constant as long as knots.x OK.
prior.E.xi <- -0.5 * log(detD)
prior.E.xi.cc <- 0.5 * E * log(cDc)
prior <- prior.E + prior.E.xi + prior.E.xi.cc
hess <- NA
if(Hessian==TRUE)
hess <- E/2 * (2*(D / cDc) - (4 * D %*% t(t(knots.y)) %*% t(knots.y) %*% D) / (cDc ^ 2) )
return (list(prior=prior, grad=(t(knots.y) %*% D) * E / cDc, hess=hess))
}
#####################################################################################
# logPriorBkgRSmooth(bkg.r, D, Hessian, Phi, matrix.FT, knots.y)
#
# returns: prior in the r-space
# arguments
logPriorBkgRSmooth <- function(bkg.r, D, Hessian=FALSE, Phi=NA, matrix.FT=NA, knots.y=NA){
E <- length(bkg.r)
cDc <- as.vector(t(bkg.r) %*% D %*% t(t(bkg.r)))
# Prior: separate terms according to which variables contribute
prior.E <- 0.5 * E * log(pi) - log(gamma(0.5 * E)) - sum(log(1:E))
# "Depends" on knots.x, but constant as long as knots.x OK.
# prior.E.xi <- -0.5 * log(detD)
prior.E.xi <- 0
prior.E.xi.cc <- 0.5 * E * log(cDc)
prior <- prior.E + prior.E.xi + prior.E.xi.cc
hess <- NA
if(Hessian==TRUE){
D <- t(Phi) %*% t(matrix.FT) %*% D %*% matrix.FT %*% Phi
hess <- E/2 * (2*(D / cDc) - (4 * D %*% t(t(knots.y)) %*% t(knots.y) %*% D) / (cDc ^ 2) )
}
return (list(likelihood=prior, hess=hess))
}
#####################################################################################
# logPriorBkgRGP(bkg.r, covMatrix, Hessian)
#
# returns: prior in the r-space due to covariance matrix
# arguments
logPriorBkgRGP <- function(bkg.r, covMatrix, Hessian=FALSE){
N <- length(bkg.r)
f <- covMatrix$factor
term.volume <- 0.5*N*log(2*pi)
term.det <- 0.5*log(covMatrix$det) - 0.5*N*log(f)
term.bkg <- f*0.5*t(bkg.r) %*% covMatrix$inv %*% t(t(bkg.r))
prior <- term.volume + term.det + term.bkg
hess <- NA
if(Hessian==TRUE)
hess <- f*covMatrix$inv
return (list(likelihood=prior, hess=hess))
}
#####################################################################################
# logLikelihoodBkg(knots.y, y, Phi, p.bkg, sigma)
#
# returns: the background-only contribution to the likelihood
# arguments
# knots.y: 'c' in Fischer et al.: vector of spline knot y-values
# y: vector of datapoints
# Phi: spline matrix taking 'c' into 'y'
# p.bkg: the probability that a single pixel contains "only" background.
# sigma: vector of experimental uncertainties
logLikelihoodBkg <- function(knots.y, y, Phi, p.bkg, sigma) {
deviation <- (y - Phi %*% t(t(knots.y))) # 'y-bkg'_i should be eps_i
deviation.norm <- deviation / sigma
likelihood <- log(p.bkg) - 0.5 * log(2 * pi) - log(sigma) - 0.5 * deviation.norm ^ 2
likelihood
}
logLikelihoodBkgAnalyt <- function(pars, x, y, p.bkg, sigma) {
deviation <- (y - bkg.analyt(pars,x)) # 'y-bkg'_i should be eps_i
deviation.norm <- deviation / sigma
likelihood <- log(p.bkg) - 0.5 * log(2 * pi) - log(sigma) - 0.5 * deviation.norm ^ 2
likelihood
}
#####################################################################################
# logLikelihoodSignal(knots.y, y, Phi, p.bkg, sigma, lambda)
#
# returns: the signal-containing contribution to the likelihood
# arguments
# knots.y: 'c' in Fischer et al.: vector of spline knot y-values
# y: vector of datapoints; 'd' at Fisher et al.
# Phi: spline matrix taking 'c' into 'y'
# p.bkg: The probability that a single pixel contains "only" background.
# sigma: vector of experimental uncertainties
# lambda: mean signal magnitude for signal-containing pixels. Either
# vector of length y or a scalar value
logLikelihoodSignal <- function(knots.y, y, Phi, p.bkg, sigma, lambda) {
rho <- sigma / lambda
z <- (y - (Phi %*% t(t(knots.y)))) / lambda
qq <- z / rho - rho
likelihood <- log(1 - p.bkg) - log(lambda) + pnorm(log.p=TRUE, q=qq) - z + 0.5 * (rho ^ 2)
likelihood
}
logLikelihoodSignalAnalyt <- function(pars, x, y, p.bkg, sigma, lambda) {
rho <- sigma / lambda
z <- (y - bkg.analyt(pars,x)) / lambda
qq <- z / rho - rho
likelihood <- log(1 - p.bkg) - log(lambda) + pnorm(log.p=TRUE, q=qq) - z + 0.5 * (rho ^ 2)
likelihood
}
#####################################################################################
# logLikelihoodGrCorr(knots.y, Phi, bkg.r, KG.inv, matrix.FT, Hessian)
#
# returns: the 'G(r)=-4*Pi*rho*r restriction' contribution to the likelihood.
# arguments
logLikelihoodGrCorr <- function(knots.y, Phi, bkg.r, KG.inv, matrix.FT, Hessian=FALSE){
M.Phi <- matrix.FT %*% Phi
KG.inv.M.Phi <- KG.inv %*% M.Phi
mu <- t(KG.inv.M.Phi) %*% bkg.r # KG.inv is a symmetric matrix
J <- t(KG.inv.M.Phi) %*% M.Phi
psi.gr <- 0.5 * (t(knots.y) %*% J) %*% t(t(knots.y)) - t(mu) %*% t(t(knots.y))
psi.gr <- sum(psi.gr)
hess <- NA
if(Hessian==TRUE)
hess <- J
return (list(likelihood=psi.gr, hess=hess))
}
#####################################################################################
#logLikelihoodGrGauss(y, knots.y, alpha, Phi, bkg.r, sigma.r, matrix.FT, Hessian=FALSE)
#
# returns: the 'G(r)=-4*Pi*rho*r restriction' contribution to the likelihood.
# arguments
logLikelihoodGrGauss <- function(y, knots.y, alpha, Phi, bkg.r, sigma.r, matrix.FT, Hessian=FALSE){
############
# y = SIGNAL - SB
deviation <- bkg.r/alpha + matrix.FT %*% ( (1-1/alpha)*y - Phi %*% t(t(knots.y)) )
deviation.norm <- deviation / sigma.r
psi.gr <- 0.5 * log(2 * pi) + log(sigma.r) + 0.5 * deviation.norm^2
psi.gr <- sum(psi.gr)
hess <- NA
if(Hessian==TRUE){
Mprime <- matrix.FT%*%Phi / (sqrt(2)*sigma.r)
hess <- 2*t(Mprime) %*% Mprime
}
return (list(likelihood=psi.gr, hess=hess))
}
#########################################################
# basisSpline(x, knots.x, knots.i, deriv)
#
# returns: the i'th spline basis function (i.e. that is nonzero (=1)
# only in the ith knot) or its derivative at points x.
# arguments
# x: x-values where we evaluate the basis functions.
# knots.x: knot positions
# knots.i: which basis function to return
# deriv: which derivative of the spline function (0 to 3)
basisSpline <- function(x, knots.x, knots.i, deriv=0){
y <- 0*knots.x
y[knots.i] <- 1
basisSpline <- splinefun(x=knots.x, y=y, method="natural")
return(basisSpline(x=x, deriv=deriv))
}
#########################################################
# basisMatrix(x, knots.x,deriv)
#
# returns: matrix whose columns are spline basis functions
# values at points x; if an arbitrary function is
# known only at points knots.x, multiplying it by
# basisMatrix results in function spline approximations
# at points x
# arguments
# x: x-values where we evaluate the basis functions.
# knots.x: knot positions
# deriv: which derivative of the spline function?
basisMatrix <- function(x, knots.x, deriv=0) {
E <- length(knots.x)
N <- length(x)
bM <- matrix(nrow=N, ncol=E)
for (i in 1:E) {
bM[, i] <- basisSpline(x=x, knots.x=knots.x, knots.i=i, deriv=deriv)
}
return(bM)
}
########################################################################################
# DMatrix(knots.x, only.trace, robust.factor)
#
# returns: the second derivative overlap matrix ('D' in Fischer et al.
# (1999)) of the spline basis functions with knots at 'knots.x'.
# arguments
# knots.x: spline knot positions
# onlyTrace: if true, returns trace(D) (a single number) instead of D.
# robustFactor: We add a constant to all eigenvalues to promote stability;
# this is the ratio of that constant to the smallest
# eigenvalue which is SUPPOSED to be nonzero.
DMatrix <- function(knots.x, onlyTrace=FALSE, robustFactor=1e-12) {
# cat("Calculating basis matrix...")
N <- length(knots.x)
ddM <- basisMatrix(x=knots.x, knots.x=knots.x, deriv=2)
# cat(" done!\n")
x.mat <- matrix(rep(knots.x, N), nrow=N)
dy <- apply(ddM, 2, diff)
dx <- apply(x.mat, 2, diff)
# slope[i, j] is the slope of the i'th segment of the j'th basis function;
# intercept[i, j] is its y-intercept
slope <- dy / dx
intercept <- ddM[-1, ] - slope * x.mat[-1, ]
d.x3 <- diff(knots.x ^ 3)
d.x2 <- diff(knots.x ^ 2)
# The for-loop code below reduces to this if we only care about diagonal
# elements (which is the case for computing the trace).
if (onlyTrace) {
return (sum((slope ^ 2) * d.x3 / 3.0 + (slope * intercept) * d.x2 +
intercept ^ 2 * diff(knots.x)))
}
# If we're this far, we need to compute the whole matrix.
DD <- matrix(0, nrow=N, ncol=N)
# computing basis function support:
N2 <- floor(N/2)
phi.N2 <- basisSpline(x=knots.x, knots.x=knots.x, knots.i=N2, deriv=2)
max.phi <- max(phi.N2)
x3.factor <- (max(knots.x)^3-min(knots.x)^3)/3
supp <- min(length(which(abs(phi.N2)>1e-9*max.phi/x3.factor)), N-1) # only i-supp to i+supp matters
# cat("Calculating overlap intergals...\n")
for (i in 1:N) {
# if (i %% 100 == 0)
# cat("...x=", knots.x[i], "\n")
relevant <- min(i+supp, N)
phi.i.supp.min <- max(1, i-ceiling(supp/2))
phi.i.supp.max <- min(N-1, i+floor(supp/2))
for (j in i:relevant) {
phi.j.supp.min <- max(1, j-ceiling(supp/2))
phi.j.supp.max <- min(N-1, j+floor(supp/2))
phi.ij.supp <- max(phi.i.supp.min, phi.j.supp.min):min(phi.i.supp.max, phi.j.supp.max)
prod.ss <- slope[phi.ij.supp, i] * slope[phi.ij.supp, j] / 3.0
prod.si <- 0.5 * (slope[phi.ij.supp, i] * intercept[phi.ij.supp, j] + slope[phi.ij.supp, j] * intercept[phi.ij.supp, i])
prod.ii <- intercept[phi.ij.supp, i] * intercept[phi.ij.supp, j]
integrals <- prod.ss * d.x3[phi.ij.supp] + prod.si * d.x2[phi.ij.supp] + prod.ii * diff(knots.x)[phi.ij.supp]
DD[i, j] <- DD[j, i] <- sum(integrals)
}
}
# cat("...done!\n")
# cat("Calculating D-matrix eigenvalues...")
# "modified" determinant (i.e., product of highest (N-2) eigenvalues)
eigenvals <- eigen(x=DD, symmetric=TRUE, only.values=TRUE)$values
# cat(" done!\n")
detD <- prod(rev(eigenvals)[-(1:2)])
# Make sure we avoid negative eigenvalues!
lowest.true.eigenval <- sort(eigenvals)[3]
return (list(matrix=DD + lowest.true.eigenval * robustFactor * diag(nrow(DD)), det=detD))
}
#########################################################
# get.bkg(x, knots.x, knots.y)
#
# returns: natural spline approximation of function knots.y
# arguments
# x: points at which function should be approximated
# knots.x: knot position
# knots.y: knot values
get.bkg <- function(x, knots.x, knots.y){
bM <- basisMatrix(x=x, knots.x=knots.x)
get.bkg <- c(bM %*% t(t(knots.y)))
return(get.bkg)
}
bkg.analyt <- function(pars, x){
bkg.analyt <- pars[1]*exp(-pars[2]*x)*x^pars[3] + pars[4]/((x-pars[5])^2+pars[6]^2)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.