**Statisfaction » R**, and kindly contributed to R-bloggers)

Today I am going to introduce the moustache target distribution (*moustarget* distribution for brievety). Load some packages first.

library(wesanderson) # on CRAN library(RShapeTarget) # available on https://github.com/pierrejacob/RShapeTarget/ library(PAWL) # on CRAN

Let’s invoke the *moustarget* distribution.

shape <- create_target_from_shape( file_name=system.file(package = "RShapeTarget", "extdata/moustache.svg"), lambda=5) rinit <- function(size) matrix(rnorm(2*size), ncol = 2) moustarget <- target(name = "moustache", dimension = 2, rinit = rinit, logdensity = shape$logd, parameters = shape$algo_parameters)

This defines a target distribution represented by a SVG file using RShapeTarget. The target probability density function is defined on and is proportional to on the segments described in the SVG files, and decreases exponentially fast to away from the segments. The density function of the *moustarget* is plotted below, a picture being worth a thousand words.

ranges <- apply(shape$bounding_box, 2, range) gridx <- seq(from=ranges[1,1], to=ranges[2,1], length.out=300) gridy <- seq(from=ranges[1,2], to=ranges[2,2], length.out=300) grid.df <- expand.grid(gridx, gridy) grid.df$logdensity <- [email protected](cbind(grid.df$Var1, grid.df$Var2), [email protected]) names(grid.df) <- c("x", "y", "z") g2d <- ggplot(grid.df) + geom_raster(aes(x=x, y=y, fill=exp(z))) + xlab("X") + ylab("Y") g2d <- g2d + xlim(ranges[,1]) + ylim(ranges[,2]) pal <- wes.palette(name = "GrandBudapest", type = "continuous") g2d <- g2d + scale_fill_gradientn(name = "density", colours = pal(50)) g2d <- g2d + theme(legend.position = "bottom", legend.text = element_text(size = 10)) g2d

There are various interesting aspects to note about this distribution. First it is very multi-modal and strongly non-Gaussian, thus providing an interesting toy problem for testing MCMC algorithms. Furthermore, sampling from the *moustarget* can be made arbitrarily difficult by pulling the moustache down, thus separating the moustache mode from the remaining probability mass around the eyes, ears and hat. Finally, note that the colours chosen to represent the density above approximately match the principal colours used in Grand Budapest Hotel by Wes Anderson. This is thanks to the awesome wesanderson package on CRAN. Obviously it is now very tempting to launch the Wang-Landau algorithm on this target with a spatial binning strategy, in order to try out the various palettes provided in wesanderson.

mhparameters <- tuningparameters(nchains = 10, niterations = 10000, storeall = TRUE) getPos <- function(points, logdensity) points[,2] explore_range <- c(-700,0) ncuts <- 20 positionbinning <- binning(position = getPos, name = "position", binrange = explore_range, ncuts = ncuts, useLearningRate = TRUE, autobinning = FALSE) pawlresults <- pawl(target = moustarget, binning = positionbinning, AP = mhparameters, verbose = TRUE) pawlchains <- ConvertResults(pawlresults, verbose = FALSE) locations <- [email protected](pawlresults$finalbins, pawlchains$X2) pawlchains$locations <- factor(locations) g <- ggplot(subset(pawlchains), aes(x=X1, y = X2, alpha = exp(logdens), size = exp(logdens), colour = locations)) + geom_point() + theme(legend.position="none") + xlab("X") + ylab("Y") g <- g + geom_hline(yintercept = pawlresults$finalbins) pal <- wes.palette(name = "GrandBudapest", type = "continuous") print(g + scale_color_manual(values = pal(21)) + labs(title = "Moustarget in Grand Budapest colours")) pal <- wes.palette(name = "Darjeeling", type = "continuous") print(g + scale_color_manual(values = pal(21)) + labs(title = "Moustarget in Darjeeling colours")) pal <- wes.palette(name = "Zissou", type = "continuous") print(g + scale_color_manual(values = pal(21)) + labs(title = "Moustarget in Zissou colours"))

**leave a comment**for the author, please follow the link and comment on their blog:

**Statisfaction » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...