# Bivariate normal probability density function # AUTHOR: M.V.Wickerhauser # DATE: 6 February 2011 # 9 February 2012 corrected for(j) loop limit # # Inputs: # # x: Vector of x-values # y: Vector of y-values # mu: Mean of the pdf [default (0,0)] # Sigma: Covariance of the pdf [default Id]. # This must be a positive definite (hence nonsingular) # symmetric 2x2 matrix. # # Output: # z: Matrix of `length(x)' rows and `length(y)' columns # containing f(x[i],y[j]) at index i,j, where # f is the bivariate normal pdf with mean `mu' # and covariance `Sigma'. # bvnpdf <- function( x, y, mu=c(0,0), Sigma=diag(c(1,1)) ) { A <- 0.5*solve(Sigma); # matrix for the quadratic form xx <- x-mu[1]; yy <- y-mu[2]; # mean-subtraction c <- 1/(2 * pi * sqrt(det(Sigma))); # normalization z<- matrix( 0, nrow=length(x), ncol<-length(y) ); # initialization for( i in 1:length(x) ) { for( j in 1:length(y) ) { xy <- rbind( xx[i], yy[j] ) z[i,j] <- c*exp(- t(xy) %*% A %*% xy ) } } return(z); }