Re: [R] (no subject)

From: Thomas Lumley <tlumley_at_u.washington.edu>
Date: Thu 06 Oct 2005 - 00:45:03 EST

On Wed, 5 Oct 2005, allan_sta_staff_sci_main_uct@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@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley@u.washington.edu	University of Washington, Seattle

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Oct 06 00:53:50 2005

This archive was generated by hypermail 2.1.8 : Sun 23 Oct 2005 - 18:21:06 EST