Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

Given

$$f_{(U_1,U_2)}(u_1,u_2)=\begin{cases} 1/2& -u_1<u_2<u_1 \text{ and } u_1 - 2 < u_2 < 2 - u_1 \text{ and } 0 < u_1 <2\\ 0& \text{otherwise}\end{cases}$$

I found that $$f_{(U_1)}(u_1) = \begin{cases} u_1 & 0<u_1<1\\2-u_1 &1<u_1<2 \\0 & \text{otherwise}\end{cases}$$

And $$f_{(U_2)}(u_2) = \begin{cases} 1+u_2 & -1<u_2<0\\1-u_2 &0<u_2<1 \\0 & \text{otherwise}\end{cases}$$

I want to justify whether $U_1$ and $U_2$ are independent or not by using R.. I am wondering what is the command is like...Could anyone please give me some directions as for how to check independence for functions in R?

share|improve this question

2 Answers 2

up vote 5 down vote accepted

For independence you want to show that the joint density factorizes, i.e. $$ f_{12}(u_1,u_2) = f_1(u_1)f_2(u_2), \;\;\; \forall u_1,u_2 \;\;\; (\dagger)$$

where for convenience I've written $f_{12}=f_{(U1,U2)}, f_1=f_{(U1)}, f_2=f_{(U2)}$.

You can't use a tool like R to prove that $U_1$ and $U_2$ are independent. This is because probability densities are continuous and can potentially have an infinite support, whereas computers can only deal with individual numbers and can't handle infinity. In your case, your densities have a finite support, i.e. they're only non-zero near the origin, so this makes things easier. But a computer can't check every single possible value of $u_1$ and $u_2$.

What you can do is to use R to check that plots, etc, are consistent with independence. As a general rule, it's a good idea to start by plotting the data (leaving other things aside, it gives you a visual check you haven't made a mistake in translating from math to code). First of all define the range of $u_1$ and $u_2$ and then use grid to get R to produce all the points $(u_1,u_2)$ that you're going to want to plot probability densities for.

u1=seq(-3,3,0.05);
u2=u1;
grid=expand.grid(u1=u1,u2=u2);

Now define your joint density $f_{12}$ and put its results alongside our points

f12 = function(u1,u2) { ifelse (-u1 < u2 & u2  < u1 & u1 - 2 < u2 & u2 < 2-u1 & 0 < u1 & u1 < 2, 0.5, 0) }
grid$f12 = with(grid,f12(u1,u2))

Note the use of R's ifelse function rather than an if statement with an accompanying else. The former works with vectors, which is what we want.

We can plot this using the lattice 3d plotting package (ggplot2 is another extremely good graphics package).

library(lattice)
wireframe(f12~u1*u2,data=grid,shade=TRUE,main="f12", scales=list(arrows=FALSE))

enter image description here Now define $f_1$ and $f_2$ and plot them to check they make sense. Note that because computers work with a set of discrete points, I had to change a < into a <= in both cases.

f1=function(u1) { ifelse(0 < u1 & u1 <= 1, u1, ifelse(1 < u1 & u1 < 2, 2-u1, 0)) }
f2=function(u2) { ifelse(-1<u2 & u2<=0, 1+u2, ifelse(0<u2 & u2<1, 1-u2, 0)) }

To plot these on top of each other, use the par command and R's regular 1d plot command.

par(mfrow=c(2,1))
plot(f1(u1)~u1, type="l")
plot(f2(u2)~u2, type="l")

enter image description here

Now define a function f1f2 to represent the product $f_1(u_1)f_2(u_2)$ on the right hand side of $(\dagger)$, put it alongside our points $(u_1,u_2)$ and plot it.

f1f2=function(u1,u2) { f1(u1)*f2(u2) }
grid$f1f2=with(grid,f1f2(u1,u2))
wireframe(f1f2~u1*u2,data=grid,shade=TRUE,main="f12", scales=list(arrows=FALSE))

enter image description here

You don't need to go any further. Clearly the product of the densities is different from the joint distribution, so $U_1$ and $U_2$ are not independent. In other circumstances if you still thought the distributions were independent, you might want to plot their differences, or simply sum their (absolute) differences, as in sum(abs(grid$f12-grid$f1f2)).

As a final note, you can actually tell from the very first plot of $f_{12}$ that $u_1$ and $u_2$ are not independent. If they were independent, then the 3d plot would have borders parallel to the $x$ and $y$ axes. As it is, you can see that depending on the values that $U_1$ takes, certain values in the support of $U_2$ are possible or impossible. So clearly there is dependence here. Contrast this with the plot of $f_1f_2$ whose borders are aligned with the axes (a necessary but not sufficient condition for independence). This is illustrated in the plot below, where I plot all the $u_2$s for which $f_{12}\neq0$ against all the $u_1$s for which $f_{12}\neq0$.

plot(grid$u2[grid$f12>0] ~ grid$u1[grid$f12>0], ylab=expression(u[2]), xlab=expression(u[1]), main=expression(paste(u[2], plain(" vs "), u[1], plain(" for "), f[12]>0)))

enter image description here

share|improve this answer
    
sum(abs(grid$f12-grid$f1f2)) will calculate the difference between $f_{U_1}*f_{U_2} - f_{U_1,U_2}(u_1, u_2)$ right...if so and if they are independent, sum(abs(grid$f12-grid$f1f2)) will give 0. If they are dependent, we will a number between 0 and 1 right since their total volumn adds up to 1, right? –  afsdf dfsaf Mar 2 at 14:19
    
@afsdfdfsaf yes, if they're independent you'll get zero (however you can't safely say if the abs sum is zero then they're independent for reasons given at the start of my answer). Re if they're dependent, it will be non-zero, yes, but it could be a lot bigger than 1, because we're not taking into account the area "dxdy" as we sum. –  TooTone Mar 2 at 14:36
    
quick note about checking if difference between 2 things is zero: computers store things as bits an are unable to store a general rational number exactly. So you should check if difference between two numbers is small (e.g. 10^(-8)) for example try b <- 0.8 - 0.2; b - 0.6 in R –  Benjamin Mar 2 at 14:42
1  
you know I'm not really sure what you mean. But also, this discussion feels like it's getting beyond the original question and it's difficult to address what you mean in an extended set of comments. Here's what I suggest: (i) decide whether you want to upvote / downvote / accept my answer; (ii) try to take stock of what you've learnt from asking your question (which is a good question), and the answers you have received (iii) formulate a new question based on taking stock, and put that question here is it's stats related or on stackoverflow if it is more related to programming in R. hth –  TooTone Mar 4 at 13:29
1  
I just want to show a 2D-plot with $U_1$ in the x-axis and $U_2$ in the y-axis...then plot the support for $f_{U_1, U_2}(U_1,U_2)$ one one graph and another 2D-plot with $U_1$ in the x-axis and $U_2$ in the y-axis...then plot for the support for $f_{U_1}(U_1)f_{U_2}(U_2)$...tell me what part you do not understand...the reason I asked that ... it is easy to from the support to show that they are dependent –  afsdf dfsaf Mar 4 at 15:41

Had you been using data to show that two continuous random variables are independent (i.e., distributions factorize as above), Hoeffding's $D$ test is one to use. It is implemented in the R Hmisc package's hoeffd function. $D$ is a measure of the difference between the product of the two vectors of ranks and the scaled bivariate rank. Nonparametric statistics texts discuss it in much more detail.

After installing the Hmisc package run

require(Hmisc)
?hoeffd    # to see documentation, which has this example below
example(hoeffd)   # run examples
# Hoeffding's test can detect even one-to-many dependency
set.seed(1)
x <- seq(-10,10,length=200)
y <- x*sign(runif(200,-1,1))
plot(x,y)   # a large X
hoeffd(x,y)

D
     x    y
x 1.00 0.06
y 0.06 1.00

avg|F(x,y)-G(x)H(y)|
       x      y
x 0.0000 0.0407
y 0.0407 0.0000

max|F(x,y)-G(x)H(y)|
       x      y
x 0.0000 0.0763
y 0.0763 0.0000

n= 200 

P
  x  y 
x     0
y  0  

The $P$-value for testing $H_0$: $X$ and $Y$ are independent prints as zero. It is $< 10^{-4}$.

share|improve this answer
    
I haven't heard of it...could you show me how to do it? –  afsdf dfsaf Mar 2 at 13:23
    
when I run it, it says the following > require(Hmisc) Loading required package: Hmisc Warning message: In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, : there is no package called ‘Hmisc’ > ?hoeffd No documentation for ‘hoeffd’ in specified packages and libraries: you could try ‘??hoeffd’ > example(hoeffd) Warning message: In example(hoeffd) : no help found for ‘hoeffd’ > set.seed(1) > x <- seq(-10,10,length=200) > y <- x*sign(runif(200,-1,1)) > plot(x,y) # a large X > hoeffd(x,y) Error: could not find function "hoeffd" How do I fix those issues –  afsdf dfsaf Mar 2 at 14:05
4  
It is best to get more familiar with the workings of R before proceeding, e.g., how to easily download and install packages. RStudio.org makes this especially easy. –  Frank Harrell Mar 2 at 14:33

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.