# [R] bivariate normal

David Winsemius dwinsemius at comcast.net
Sat Jul 28 01:38:25 CEST 2012

```On Jul 27, 2012, at 10:45 AM, Anders Holm wrote:

> Dear list members
>
>
>
> I need a function that calculates the bivariate normal distribution
> for each observation. It is part of a likelihood function and I have
> 1000's of cases. As I understand it I cannot use packages like
> "mvtnorm" because it requres a covariance matrix of the same
> dimension as the number of observations. Basically what I need is a
> function that takes as arguments a n*2 matrix of bivariate values
> given a common mean and covariance matrix, where n is the number of
> cases and which returns a n*1 vector of the probabilities of the
> bivariate normal distribution of the n*2 vector of values.
>

I do not think you are reading the documentation for the d- and p-
functions in mvtnorm correctly:

Here is the 2-dimensional case using pmvnorn (since you asked for
probabilities)  and the result is as expected:

> n <- 2
> mean <- rep(0, 2)
> corr <- diag(2)
> corr[lower.tri(corr)] <- 0
> corr[upper.tri(corr)] <- 0

Half of the xy plane of a 2d-mvnorm with 0 covariances should give 0.5

> prob <- pmvnorm( lower= c(-Inf,-Inf), upper= c( 0, Inf), mean, corr)
> prob
 0.5
attr(,"error")
 2e-16
attr(,"msg")
 "Normal Completion"

And a single quadrant should give 0.25.

> prob <- pmvnorm(lower=c(-Inf,-Inf), upper=c( 0, 0), mean, corr)
> prob
 0.25
attr(,"error")
 1e-15
attr(,"msg")
 "Normal Completion"

I guessing that you want to provide the data to pmvnorm with apply()
using -Inf as the 'lower' limits and your datapoints as the 'upper'
values:

Perhaps:

apply(n_x_2_values,
>
>
>
> all the best
>
> Anders Holm
>
> Department of Education
>
> Aarhus University
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help