bivariate normal

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
[1] 0.5
attr(,"error")
[1] 2e-16
attr(,"msg")
[1] "Normal Completion"

And a single quadrant should give 0.25.

 > prob <- pmvnorm(lower=c(-Inf,-Inf), upper=c( 0, 0), mean, corr)
 > prob
[1] 0.25
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "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
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
Alameda, CA, USA