Discrete approximation to bivariate normal distribution
I would like to do a discrete approximation to the bivariate normal distribution. That is, I would like to compute a matrix where each entry is the probability of hitting one of the small squares in the bottom image.
Here's what I've done so far.
library(mvtnorm)
library(graphics)
euclide = function(x,y){sqrt(x^2+y^2)}
maxdist = 40
sigma = diag(2)
m = matrix(0,ncol=maxdist*2 + 1, nrow=maxdist*2 + 1)
for (row in -maxdist:maxdist){
for (col in -maxdist:maxdist){
if ( euclide(abs(row), abs(col)) < maxdist ){
lower = c(row-0.5, col-0.5)
upper = c(row+0.5, col+0.5)
p = pmvnorm(lower = lower , upper = upper, mean = c(0,0), sigma = sigma)
} else {
p = 0
}
m[row + maxdist + 1,col + maxdist + 1] = p
}
}
m = m[rowSums(m)!=0,colSums(m)!=0]
contour(m, levels = exp(-20:0), xlim=c(0.3,0.7), ylim=c(0.3,0.7))
It works fine. This is quite slow (for large ones maxdist
), although I hope to improve my computational time. But this is not my main problem ...
The main problem is that with my method I cannot change the number of small squares close to the center to get closer to the mean. I can only add squares to the surrounding. In other words, I would like to be able to set the variance of both axes of the bivariate normal distribution.
source to share
Here's a simple implementation. As @DanielJohnson says you can just use the cdf form one dimensional, but it should be the same as using pmvnorm
as shown below. The using version pnorm
is much faster.
## Choose the matrix dimensions
yticks <- xticks <- seq(-3, 3, length=100)
side <- diff(yticks[1:2]) # side length of squares
sigma <- diag(2) # standard devs. for f2
mu <- c(0,0) # means
## Using pnorm
f <- Vectorize(function(x, y, side, mu1, mu2, s1, s2)
diff(pnorm(x+c(-1,1)*side/2, mu1, s1)) * diff(pnorm(y+c(-1,1)*side/2, mu2, s2)),
vec=c("x", "y"))
## Using pmvnorm
f2 <- Vectorize(function(x, y, side, mu, sigma)
pmvnorm(lower=c(x,y)-side/2, upper=c(x,y)+side/2, mean=mu, sigma=sigma),
vec=c("x", "y"))
## get prob. of squares, mu are means, s are standards devs.
mat <- outer(xticks, yticks, f, side=side, mu1=0, mu2=0, s1=1,s2=1)
mat2 <- outer(xticks, yticks, f2, side=side, mu=mu, sigma=sigma)
## test equality
all(abs(mat2-mat) < 1e-11) # TRUE
all.equal(mat2, mat) # TRUE
## See how it looks
library(lattice)
persp(mat, col="lightblue", theta=35, phi=35, shade=0.1)
source to share
I'm not an R person, but I'm pretty sure there is a CDF for normal distribution. If you want to literally plot the probability matrix of falling per square, we can use this CDF to get the answer. Since the bivariate normal distribution has independent marginal distributions, the question here is simply to ask two questions for each square described by the [x_left, x_right] and [y_left, y_right] axis positions:
- What is the probability that a 1D normal random variable is in the interval [x_left, x_right]?
- What is the probability that another independent 1D normal random variable is in the interval [y_left, y_right]?
Since the two are independent, the total probability for a square is:
P = (CDF(x_right) - CDF(x_left))*(CDF(y_right) - CDF(y_left))
This is the exact answer, so computing time shouldn't be an issue!
EDIT: I must also say that you can select a grid with a lot of ticks on each axis close to zero to get the resolution you want. The above probability formula for each square is kept.
source to share