Hello,
I try to create a persp-plot for the bivariate normal distribution.
There are a few solutions online ([1],[2]), but I did not want to
hard-code the formula like [2] did or use a kernel density estimate like
[1].
My plan was to use dmvnorm, and use outer to calculate the z-dimension
for each point. Here is what I tried:
library(mvtnorm)
x<-seq(0,1,len=40)
y<-x
g<-function(a,b) {
dmvnorm(x=c(a,b))
}
z<-outer(x,y,g)
Which results in the following error:
Error in dim(robj) <- c(dX, dY) :
dims [product 1600] do not match the length of object [1]
What causes this error? If I set g to be a simple function like:
g<-function(a,b){
a*b
}
Everything works fine. Why does it fail when I pass the parameters to
dmvnorm?
Thank you for your time,
Severin
[1]https://stat.ethz.ch/pipermail/r-help/2003-September/038314.html
[2]http://www.stat.ucl.ac.be/ISpersonnel/lecoutre/stats/fichiers/_gallery.pdf
On Mar 28, 2010, at 3:19 AM, Severin Kacianka wrote:> Hello, > > I try to create a persp-plot for the bivariate normal distribution. > There are a few solutions online ([1],[2]), but I did not want to > hard-code the formula like [2] did or use a kernel density estimate > like [1]. > > My plan was to use dmvnorm, and use outer to calculate the z- > dimension for each point. Here is what I tried: > library(mvtnorm) > x<-seq(0,1,len=40) > y<-x > g<-function(a,b) { > dmvnorm(x=c(a,b))The construction "c(a,b)" would in the later function call turn those two 40 element vectors into an 80 element vector. See if you get better success with cbind.> } > z<-outer(x,y,g) > > Which results in the following error: > Error in dim(robj) <- c(dX, dY) : > dims [product 1600] do not match the length of object [1] > > What causes this error? If I set g to be a simple function like: > g<-function(a,b){ > a*b > } > Everything works fine. Why does it fail when I pass the parameters > to dmvnorm? > > Thank you for your time, > Severin > > [1]https://stat.ethz.ch/pipermail/r-help/2003-September/038314.html > [2]http://www.stat.ucl.ac.be/ISpersonnel/lecoutre/stats/fichiers/_gallery.pdf >David Winsemius, MD West Hartford, CT
Hello David and Dennis,
thank you two for your replies. You not only helped me solve my problem,
but also helper me to understand my problem better.
Here is my -now working- code:
x<-y<-seq(-4,4,len=40)
g<-function(a,b) {
dmvnorm(x=cbind(a,b),sigma=matrix(c(4,2,2,3), ncol = 2))
}
z<-outer(x,y,g)
persp(x,y,z)
Thank you a lot!
Severin
On Mar 28, 2010, at 9:36 AM, Severin Kacianka wrote:> Hello David and Dennis, > > thank you two for your replies. You not only helped me solve my > problem, but also helper me to understand my problem better. > Here is my -now working- code: > > x<-y<-seq(-4,4,len=40) > > g<-function(a,b) { > dmvnorm(x=cbind(a,b),sigma=matrix(c(4,2,2,3), ncol = 2)) > } > z<-outer(x,y,g) > persp(x,y,z)From that point, you will get a better idea of what you are looking at with some additional views: > persp(x,y,z) > persp(x,y,z, theta=-45) > persp(x,y,z, theta=-60) > persp(x,y,z, theta=-90) > persp(x,y,z, theta=-135)>-- David Winsemius, MD West Hartford, CT