Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Code to generate wallpaper (bivariate normal distribution)
library(mvtnorm)
### correlation for the bivariate normal distribution
rho <- 0.4
### number of points for x and y at which to evaluate the density of the bivariate normal distribution
n <- 501
sigma <- matrix(c(1,rho,rho,1), nrow=2)
x <- seq(-3, 3, length=n)
y <- seq(-3, 3, length=n)
z <- matrix(NA, nrow=n, ncol=n)
for (i in seq_along(x)) {
for (j in seq_along(y)) {
z[i,j] <- dmvnorm(c(x[i],y[j]), sigma=sigma)
}
}
### used width and height equal to twice the resolution, which looks nice (but adapt as needed)
jpeg("bivariate_normal.jpg", width=2732, height=1536, quality=95, bg="gray10", type="cairo")
### contour plot
par(mar=c(0,0,15,0))
lvls <- c(seq(.02, .20, by=.01))
cols <- colorRampPalette(c("gray25", "gray95"))(length(lvls))
contour(z=z, axes=F, drawlabels=FALSE, levels=lvls, col=cols)
par(new=TRUE)
### bivariate normal surface
par(mar=c(13,0,9,0))
nrz <- nrow(z)
ncz <- ncol(z)
jet.colors <- colorRampPalette(c(rgb(36,36,36,maxColorValue=255), "gray80"))
nbcol <- 1000
color <- jet.colors(nbcol)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
facetcol <- cut(zfacet, nbcol)
persp(x, y, z, theta=45, phi=30, r=10, box=FALSE, shade=0.3, col=color[facetcol], border=NA, expand=0.6, ltheta=-180, lphi=-10)
### add pdf text
text(0, 0.07, col="gray70", pos=1, cex=2.6,
expression(frac(1, 2 ~ pi ~ sigma[x] ~ sigma[y] ~ sqrt(1 - rho^2)) ~
exp ~ bgroup("[", -frac(1,2*(1-rho^2)) ~ bgroup("(", frac((x-mu[x])^2, sigma[x]^2) ~ + ~
frac((y-mu[y])^2, sigma[y]^2) ~ - ~
frac(2 ~ rho ~ (x-mu[x]) ~ (y-mu[y]), sigma[x] ~ sigma[y]),
")"),
"]")
)
)
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.