e_dist | R Documentation |
e_dist
computes the energy distance between points D
and a target distribution (or big dataset) F
. The cross-term E[||X-X'||]
, X,X'~F
is NOT computed in e_dist
for computational efficiency, since this is not needed for optimizing D
. The target distribution or big dataset can be set using dist.str
or dist.samp
, respectively.
e_dist(D, dist.str=NA, dist.param=vector("list",ncol(D)),
nsamp=1e6, dist.samp=NA)
D |
An |
dist.str |
A |
dist.param |
A
|
nsamp |
Number of samples to draw from |
dist.samp |
An |
Szekely, G. J. and Rizzo, M. L. (2013). Energy statistics: A class of statistics based on distances. Journal of Statistical Planning and Inference, 143(8):1249-1272.
## Not run:
#############################################################
# Generate 25 SPs for the 2-d i.i.d. N(0,1) distribution
#############################################################
n <- 25 #number of points
p <- 2 #dimension
D <- sp(n,p,dist.str=rep("normal",p))
Drnd <- matrix(rnorm(n*p),ncol=p)
e_dist(D$sp,dist.str=rep("normal",p)) #smaller
e_dist(Drnd,dist.str=rep("normal",p))
#############################################################
# Support points for big data reduction: Franke's function
#############################################################
# library(MHadaptive) # Package archived, but you can use your favorite MCMC sampler
# #Generate MCMC samples
# li_func <- franke2d #Desired log-posterior
# ini <- c(0.5,0.5) #Initial point for MCMc
# NN <- 1e5 #Number of MCMC samples desired
# burnin <- NN/2 #Number of burn-in runs
# mcmc_franke <- Metro_Hastings(li_func, pars=ini, prop_sigma=0.05*diag(2),
# iterations=NN, burn_in=burnin)
data(mcmc_franke) # Loading MCMC sample from data
#Use modified Franke's function as posterior
franke2d <- function(xx){
if ((xx[1]>1)||(xx[1]<0)||(xx[2]>1)||(xx[2]<0)){
return(-Inf)
}
else{
x1 <- xx[1]
x2 <- xx[2]
term1 <- 0.75 * exp(-(9*x1-2)^2/4 - (9*x2-2)^2/4)
term2 <- 0.75 * exp(-(9*x1+1)^2/49 - (9*x2+1)/10)
term3 <- 0.5 * exp(-(9*x1-7)^2/4 - (9*x2-3)^2/4)
term4 <- -0.2 * exp(-(9*x1-4)^2 - (9*x2-7)^2)
y <- term1 + term2 + term3 + term4
return(2*log(y))
}
}
#Generate ncur SPs
ncur <- 50
D <- sp(ncur,2,dist.samp=mcmc_franke$trace)$sp
Drnd <- mcmc_franke$trace[sample(1:nrow(mcmc_franke$trace),n,F),]
e_dist(D,dist.samp=mcmc_franke$trace) #smaller
e_dist(Drnd,dist.samp=mcmc_franke$trace)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.