2014-02-18 70 views
3

我想在R中創建一個圖形。它包含矢量變量(x,y)的二元正態分佈的等值線圖以及邊緣f(x),f (y)基條件分佈f(y | x)和通過條件值X = x的線(它將是一個簡單的abline(v = x))。 我已經得到了輪廓和abline:具有邊際和條件密度的雙變量法線

http://i62.tinypic.com/noyfls.png

,但我不知道該怎麼繼續。

這裏是我迄今爲止所使用的代碼:

bivariate.normal <- function(x, mu, Sigma) { 
    exp(-.5 * t(x-mu) %*% solve(Sigma) %*% (x-mu))/sqrt(2 * pi * det(Sigma)) 
} 

mu <- c(0,0) 
Sigma <- matrix(c(1,.8,.8,1), nrow=2) 
x1 <- seq(-3, 3, length.out=50) 
x2 <- seq(-3, 3, length.out=50) 

z <- outer(x1, x2, FUN=function(x1, x2, ...){ 
      apply(cbind(x1,x2), 1, bivariate.normal, ...) 
      }, mu=mu, Sigma=Sigma) 

contour(x1, x2, z, col="blue", drawlabels=FALSE, nlevels=4, 
     xlab=expression(x[1]), ylab=expression(x[2]), lwd=1) 
abline(v=.7, col=1, lwd=2, lty=2) 
text(2, -2, labels=expression(x[1]==0.7)) 
+1

請提供[重複的例子(http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example),其中包括迄今爲止使用的代碼。 – Thomas

回答

2

這本來是有益的,如果你提供了計算邊緣分佈的功能。我可能已經得到了邊緣分佈函數錯了,但我認爲這得到你想要的東西:

par(lwd=2,mgp=c(1,1,0)) 
# Modified to extract diagonal. 
bivariate.normal <- function(x, mu, Sigma) 
    exp(-.5 * diag(t(x-mu) %*% solve(Sigma) %*% (x-mu)))/sqrt(2 * pi * det(Sigma)) 

mu <- c(0,0) 
Sigma <- matrix(c(1,.8,.8,1), nrow=2) 
x1 <- seq(-3, 3, length.out=50) 
x2 <- seq(-3, 3, length.out=50) 

plot(1:10,axes=FALSE,frame.plot=TRUE,lwd=1) 

# z can now be calculated much easier. 
z<-bivariate.normal(t(expand.grid(x1,x2)),mu,Sigma) 
dim(z)<-c(length(x1),length(x2)) 
contour(x1, x2, z, col="#4545FF", drawlabels=FALSE, nlevels=4, 
     xlab=expression(x[1]), ylab=expression(x[2]), lwd=2,xlim=range(x1),ylim=range(x2),frame.plot=TRUE,axes=FALSE,xaxs = "i", yaxs = "i") 
axis(1,labels=FALSE,lwd.ticks=2) 
axis(2,labels=FALSE,lwd.ticks=2) 
abline(v=.7, col=1, lwd=2, lty=2) 
text(2, -2, labels=expression(x[1]==0.7)) 

# Dotted 
f<-function(x1,x2) bivariate.normal(t(cbind(x1,x2)),mu,Sigma) 
x.s<-seq(from=min(x1),to=max(x1),by=0.1) 
vals<-f(x1=0.7,x2=x.s) 
lines(vals-abs(min(x1)),x.s,lty=2,lwd=2) 

# Marginal probability distribution: http://mpdc.mae.cornell.edu/Courses/MAE714/biv-normal.pdf 
# Please check this, I'm not sure it is correct. 
marginal.x1<-function(x) exp((-(x-mu[1])^2)/2*(Sigma[1,2]^2))/(Sigma[1,2]*sqrt(2*pi)) 
marginal.x2<-function(x) exp((-(x-mu[1])^2)/2*(Sigma[2,1]^2))/(Sigma[2,1]*sqrt(2*pi)) 

# Left side solid 
vals<-marginal.x2(x.s) 
lines(vals-abs(min(x1)),x.s,lty=1,lwd=2) 

# Bottom side solid 
vals<-marginal.x1(x.s) 
lines(x.s,vals-abs(min(x2)),lty=1,lwd=2) 

enter image description here

+0

謝謝,但我沒有創建邊際分佈的功能。你在上面看到的圖(我認爲是在Matlab中創建的)是我試圖在R中重新創建的圖,而我沒有它的代碼。 – Umberto

+0

我不是指代碼,而是創建邊際分佈的數學函數。我在代碼中包含了一個鏈接,顯示了二元正態分佈的邊際分佈的派生,您可以在那裏查看它。我想我假設你理解如何計算邊際,你只是不知道如何編碼。如果你喜歡這個答案,你可以點擊它左邊的複選框。 – nograpes

+0

是的,謝謝! – Umberto