hi all
why does the following not work???
this was someone elses code and i couldnt explain why it doesn't work.
m=matrix(c(0,0),2,1)
v=matrix(c(1,0,0,1),2,2)
Y=function(X1,X2,mu=m,V=v)
{
X=matrix(c(X1,X2),2,1)
a=(1/((2*pi)*sqrt(det(V))))*exp((-0.5)*(t(X-mu)%*%solve(V)%*%(X-mu)))
a[1]
}
x1=seq(-1,1)
x2=x1
Z=outer(x1,x2,FUN="Y",mu=m,V=v)
persp(x1,x2,Z)
my code:
BINORMAL<-function(varx=1,vary=1,covxy=0,meanx=0,meany=0)
{
#the following function plots the density of a bi variate normal distribution
covXY<-matrix(c(varx,covxy,covxy,vary),2,2)
A<-solve(covXY)
#up<-max(meanx+4*varx^.5,meanx-4*varx^.5,meany+4*vary^.5,meany-4*vary^.5)
#x <- seq(-up,up,length=50)
#y <- x
x <- seq(meanx-3*varx^.5,meanx+3*varx^.5,length=50)
y <- seq(meany-3*vary^.5,meany+3*vary^.5,length=50)
f <- function(x,y,...)
{
detA<-det(A)
quadForm<-A[1,1]*(x-meanx)^2+2*A[1,2]*(x-meanx)*(y-meany)+A[2,2]*(y-meany)^2
K<-sqrt(detA)/(2*pi)
exp(-0.5*quadForm)*K
}
z <- outer(x, y, f)
par(mfrow=c(1,2))
persp(x, y, z,theta = 30, phi = 30,col="white",main="BI-VARIATE
NORMAL
DISTRIBUTION")
contour(x,y,z,main=paste("xy plot, corr(X,Y)=
",(covxy/(varx*vary)^.5)))
print("NOTE -sqrt(varx*vary)<=covxy<=sqrt(varx*vary)")
#print(A)
}
BINORMAL(varx=1,vary=1,covxy=0,meanx=0,meany=0)
thanx
/allan
On Wed, 5 Oct 2005, allan_sta_staff_sci_main_uct at mail.uct.ac.za wrote:> hi all > > why does the following not work??? > > this was someone elses code and i couldnt explain why it doesn't work.I think this is a case of FAQ 7.17: Why does outer() behave strangely with my function? -thomas> m=matrix(c(0,0),2,1) > v=matrix(c(1,0,0,1),2,2) > > Y=function(X1,X2,mu=m,V=v) > { > X=matrix(c(X1,X2),2,1) > a=(1/((2*pi)*sqrt(det(V))))*exp((-0.5)*(t(X-mu)%*%solve(V)%*%(X-mu))) > a[1] > } > > x1=seq(-1,1) > x2=x1 > > Z=outer(x1,x2,FUN="Y",mu=m,V=v) > > persp(x1,x2,Z) > > > > > my code: > > BINORMAL<-function(varx=1,vary=1,covxy=0,meanx=0,meany=0) > { > #the following function plots the density of a bi variate normal distribution > > covXY<-matrix(c(varx,covxy,covxy,vary),2,2) > A<-solve(covXY) > > #up<-max(meanx+4*varx^.5,meanx-4*varx^.5,meany+4*vary^.5,meany-4*vary^.5) > #x <- seq(-up,up,length=50) > #y <- x > > x <- seq(meanx-3*varx^.5,meanx+3*varx^.5,length=50) > y <- seq(meany-3*vary^.5,meany+3*vary^.5,length=50) > > f <- function(x,y,...) > { > detA<-det(A) > quadForm<-A[1,1]*(x-meanx)^2+2*A[1,2]*(x-meanx)*(y-meany)+A[2,2]*(y-meany)^2 > K<-sqrt(detA)/(2*pi) > exp(-0.5*quadForm)*K > } > > z <- outer(x, y, f) > > par(mfrow=c(1,2)) > persp(x, y, z,theta = 30, phi = 30,col="white",main="BI-VARIATE NORMAL > DISTRIBUTION") > contour(x,y,z,main=paste("xy plot, corr(X,Y)= ",(covxy/(varx*vary)^.5))) > > print("NOTE -sqrt(varx*vary)<=covxy<=sqrt(varx*vary)") > #print(A) > } > BINORMAL(varx=1,vary=1,covxy=0,meanx=0,meany=0) > > thanx > > /allan > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! R-project.org/posting-guide.html >Thomas Lumley Assoc. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle