Nothing
#' Full Bayes Power Prior with Correlated Weight Parameters
#'
#' @param x number of historical successes
#' @param n number historical patients
#' @param verbose Print messages
#' @param length Number of points to evaluate density at
#' @param mc.cores Number of cores for parallel
#' @param samples Number of Monte Carlo samples
#' @param p.prior.a shape1 parameter for initial beta prior on probability
#' @param p.prior.b shape2 parameter for initial beta prior on probability
#' @param focus List of triples specifying regions to focus on, eg peaks. Specified as list(c(lower, upper, length))
#' @param d.prior.cor Correlation parameter for transformed multivariate normal prior on weights
#'
#' @return A density function
#' @export
#'
#' @importFrom mvtnorm rmvnorm
#'
binom.PP.FB.MC.COR <- function(x, n, verbose=FALSE, length=200, d.prior.cor=0, p.prior.a=1, p.prior.b=1, mc.cores=1, samples=10000, focus ){
n.hist <- length(x)
#Where to calculate density at
p <- seq(0.0001, .9999, len=length)
if(!missing(focus)){
p2 <- unlist(sapply(focus, function(f) seq(from=f[1], to=f[2], length.out = f[3])))
p <- unique(sort(c(p,p2)))
}
dens <- mclapply(p, function(PROB){
# print(PROB)
eval.f <- function(d) dbeta(PROB, p.prior.a+sum(d*x), p.prior.b+sum(d*(n-x)))
#correlation matrix
sigma <- matrix(d.prior.cor, ncol=n.hist, nrow = n.hist)
diag(sigma) <- rep(1, n.hist)
mean(apply(
pnorm(rmvnorm(samples, mean=rep(0,n.hist), sigma=sigma)),
1,eval.f))
},mc.cores=mc.cores)
#
#
# a <- smooth.spline(p, dens)
# dens2 <- density(rep(p,unlist(dens)*10000),from=0,to=1, bw=.01,)
# opts <- optimr::optimr(c(2,3), function(Q) sum((dens2$y - dbeta(dens2$x, Q[1],Q[2]))^2),
# lower=c(1,1), upper=c(100,100), method="L-BFGS-B")
# dd <- density(rep(a$x,a$y*10000), from=0, to=1)
# lines(dd)
# points(a$x,a$y)
# plot(a$x,a$y, xlim=c(0,.5), ylim=c(0,1e-5))
#
# dens2$y[dens2$y<1e-15] <- dbeta(dens2$x[dens2$y<1e-15],opts$par[1],opts$par[2])
# f <- splinefun(smooth.spline(x=p, y=pmax(dens,0)))
# k <- integrate(f, 0,1)$value
# f <- splinefun(x=p, y=unlist(dens)/k)
# g <- function(p,X) ifelse(0<=p&p<=1, f(p),0)
#
# g <- approxfun(dens2, yleft=0, yright=0)
g <- approxfun(c(0,p,1), c(0,unlist(dens),0))
formals(g) <- alist(v=,...=)
return(g)
}
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.