How do you verify independence of a pair of uniform random variables in MATLAB? If $A$ is a subset of $R$ and $X$ is a random variable.  I have two variables $X_1$ and $X_2$.   $I$ being $1$ if $X$ in subset $A$, and $0$ if not in $A$.  Let $U$~$U(0;1)$ and determine if this pair is independent.  Verify your claim using simulation in Matlab.
$$ X_1 = I_U \epsilon\left[\left.0,\frac{1}{3}\right.\right), X_2 = I_U\epsilon\left[\left.\frac{1}{3},\frac{2}{3}\right.\right)$$
I determined that this pair is not independent because  $P[X_1=1]=P[U∈[0,\frac{1}{3})]=\frac{1}{3}$ and similarly for $X_2$. However, $P[X_1=1,X_2=1]=0$.  Now I am quite unfamiliar with MATLAB.  To verify, do I call on $rand$ many times and tally how many times the values fall within the bounds of each random variable and make a histogram? Or do I make a plot of various values ranging from $(0,1)$ and show how the two random variables act at each value? Any type of suggestions will help get me started please!  I do have about 5 different pairs of random variables, some being independent and some not.  I must run and verify each pair. 
 A: Comment (continued):  Here are three examples in which both analytic methods and technology may be useful. Both involve taking $m = 20,000$ samples of size $n = 5$ from a particular population and trying to determine whether the $m$ sample means $\bar X$ are (or are not) stochastically independent of $m$ sample standard deviations $S.$ (I say 'stochastically' independent because the equation for $S$ involves $\bar X,$ so they can't be 'functionally' independent.)
Example 1 involves sampling from a normal population. Here there is a theorem that guarantees independence of $\bar X$ and $S.$ 
Example 2 involves sampling from exponential data. Here a plot of $S$ against $\bar X$ suggests marked correlation, and association (non-independence) is easy to prove. 
Example 3 involves samples from $\mathsf{Beta}(.1, .1)$. Here $\bar X$ and $S$ are uncorrelated, but highly dependent.

Respective sample correlations are $r \approx 0, r > 0, r \approx 0.$
In the center plot it is clear that the sample mean
can be near .5 and the sample SD can be near 3, but
the two events are disjoint. So that $\bar X$ and $S$ cannot be independent. (It is possible to derive the
equation of the diagonal line at the top edge of the data cloud.) 
At right, it is obvious by symmetry that
the population correlation must be 0, but the plot
makes it clear that $\bar X$ and $S$ cannot be
independent. (The data lie near the corners, edges,
and faces of the 5-dimensional unit hypercube. When the data are reduced to the two-dimensional $(\bar X, S),$
one can still discern images of some corners and edges.)
set.seed(421)     # for reproducibility
m = 20000;  n = 5
par(mfrow=c(1,3))
x = rnorm(m*n); DTA = matrix(x, nrow=m)
 a = rowMeans(DTA);  s = apply(DTA, 1, sd)
 plot(a, s, pch=".", main="Standard Normal")
 cor(a,s)
 [1] -0.001354763    # consistent with 0
x = rexp(m*n); DTA = matrix(x, nrow=m)
 a = rowMeans(DTA);  s = apply(DTA, 1, sd)
 plot(a, s, pch=".", main="Standard Exponential")
 cor(a,s)
 [1] 0.7695967 
x = rbeta(m*n, .1,.1); DTA = matrix(x, nrow=m)
 a = rowMeans(DTA);  s = apply(DTA, 1, sd)
 plot(a, s, pch=".", main="Standard Normal")
 cor(a,s)
 [1] -0.008673277   # consistent with 0
par(mfrow=c(1,1))

Addendum: Slight variation of your example in which discrete points have differing probabilities. Use 'jittering' (small random offsets from discrete points) to aid visualization of probabilities at points.
set.seed(2019)
m = 10^4    # for good graph, don't use too many
            # for accurate est of r, use very many
x = runif(m);  x1 = (x<.3);  x2 = (x>=.8)
cor(x1, x2)
[1] -0.327649
jit1 = runif(m, -.25, .25); jit2 = runif(m, -.25, .25)
 plot(x1+jit1, x2+jit2, pch=".")


