Random Vectors and MATLAB

Summary: The systematic formulation in the previous module Minterms shows that each Boolean combination, as a union of minterms, can be designated by a vector of zero-one coefficients. A coefficient one in the ith position (numbering from zero) indicates the inclusion of minterm Mi in the union. We formulate this pattern carefully below and show how MATLAB logical operations may be utilized in problem setup and solution.

m-procedures for a pair of simple random variables

We examine, first, calculations on a pair of simple random variables X,Y, considered jointly. These are, in effect, two components of a random vectorW=(X,Y), which maps from the basic space Ω to the plane. The induced distribution is on the (t,u)-plane. Values on the horizontal axis (t-axis) correspond to values of the first coordinate random variable X and values on the vertical axis (u-axis) correspond to values of Y. We extend the computational strategy used for a single random variable.

First, let us review the one-variable strategy. In this case, data consist of values t_i and corresponding probabilities P(X=t_i) arranged in matrices

X=[t₁,t₂,⋯,t_n]andPX=[P(X=t₁),P(X=t₂),⋯,P(X=t_n)]

(1)

To perform calculations on Z=g(X), we we use array operations on X to form a matrix

G=[g(t₁)g(t₂)⋯g(t_n)]

(2)

which has g(t_i) in a position corresponding to P(X=t_i) in matrix PX.

Basic problem. Determine P(g(X)∈M), where M is some prescribed set of values.

Use relational operations to determine the positions for which g(t_i)∈M. These will be in a zero-one matrix N, with ones in the desired positions.
Select the P(X=t_i) in the corresponding positions and sum. This is accomplished by one of the MATLAB operations to determine the inner product of N and PX

We extend these techniques and strategies to a pair of simple random variables, considered jointly.

The data for a pair {X,Y{ of random variables are the values of X and Y, which we may put in row matrices
X=[t₁t₂⋯t_n]andY=[u₁u₂⋯u_m]
(3)
and the joint probabilities P(X=t_i,Y=u_j) in a matrix P. We usually represent the distribution graphically by putting probability mass P(X=t_i,Y=u_j) at the point (t_i,u_j) on the plane. This joint probability may is represented by the matrix P with elements arranged corresponding to the mass points on the plane. Thus
PhaselementP(X=t_i,Y=u_j)atthe(t_i,u_j)position
(4)
To perform calculations, we form computational matrices t and u such that — t has element t_i at each (t_i,u_j) position (i.e., at each point on the ith column from the left) — u has element u_j at each (t_i,u_j) position (i.e., at each point on the jth row from the bottom) MATLAB array and logical operations on t,u,P perform the specified operations on t_i,u_j, andP(X=t_i,Y=u_j) at each (t_i,u_j) position, in a manner analogous to the operations in the single-variable case.

Formation of the t and u matrices is achieved by a basic setup m-procedure called jcalc. The data for this procedure are in three matrices: X=[t₁,t₂,⋯,t_n] is the set of values for random variable X Y=[u₁,u₂,⋯,u_m] is the set of values for random variable Y, and P=[p_ij], where p_ij=P(X=t_i,Y=u_j). We arrange the joint probabilities as on the plane, with X-values increasing to the right and Y-values increasing upward. This is different from the usual arrangement in a matrix, in which values of the second variable increase downward. The m-procedure takes care of this inversion. The m-procedure forms the matrices t and u, utilizing the MATLAB function meshgrid, and computes the marginal distributions for X and Y. In the following example, we display the various steps utilized in the setup procedure. Ordinarily, these intermediate steps would not be displayed.

EXAMPLE 1: Setup and basic calculations

>> jdemo4                         % Call for data in file jdemo4.m
>> jcalc                          % Call for setup procedure
ix of VALUES of X  X
Enter row matrix of VALUEEnter JOINT PROBABILITIES (as on the plane)  P
Enter row mat
rS of Y  Y
 Use array operations on matrices X, Y, PX, PY, t, u, and P
0.0360    0.0198    0.0297    0.0209    0.0180
    0.0372    0.055>> disp(P)                        % Optional call for display of P
   
 8    0.0837    0.0589    0.0744
    0.0516    0.0774    0.1161    0.0817    0.1032
    0.0264    0.0270    0.0405    0.0285    0.0132
                    % Optional call for display of PY
PY =  0.1356 >> PX                             % Optional call for display of PX
PX =  0.1512    0.1800    0.2700    0.1900    0.2088
>> PY        
    0.4300    0.3100    0.1244
- - - - - - - - - -               % Steps performed by jcalc
>> PX = sum(P)                    % Calculation of PX as performed by jcalc
[t,u] = meshgrid(X,fliplr(Y)); % Formation of t, u mPX =  0.1512    0.1800    0.2700    0.1900    0.2088
>> PY = fliplr(sum(P'))           % Calculation of PY (note reversal)
PY = 0.1356    0.4300    0.3100    0.1244
>>
 atrices (note reversal)
>> disp(t)                        % Display of calculating matrix t
    -3     0     1     3     5    % A row of X-values for each value of Y
    -3     0     1     3     5

     1     1     1     1     -3     0     1     3     5
    -3     0     1     3     5
>> disp(u)                        % Display of calculating matrix u
     2     2     2     2     2    % A column of Y-values (increasin
g1    % upward) for each value of X
     0     0     0     0     0

-2    -2    -2    -2    -2

Suppose we wish to determine the probability P(X²−3Y≥1). Using array operations on t and u, we obtain the matrix G=[g(t_i,u_j)].

>> G = t.^2 - 3*u                 % Formation of G = [g(t_i,u_j)] matrix
G  = 3    -6    -5     3    19
     6    -3    -2     6    22
>> M = G >=  1                9     0     1     9    25
    15     6     7    15    31
    % Positions where G >= 1
1
     1     1     1     1    M =  1     0     0     1     1
     1     0     0     1     1
     1     0     1     1    
  1
>> pM = M.*P                    % Selection of probabilities
pM =
    0.0360         0         0    0.0209    0.0180
0    0.0405    0.0285    0.0132
>> PM = total(pM) 0.0372         0         0    0.0589    0.0744
    0.0516         0    0.1161    0.0817    0.1032
    0.0264    0.02
7              % Total of selected probabilities

PM =  0.7336                    % P(g(X,Y) >= 1)

In Example 3 from "Random Vectors and Joint Distributions" we note that the joint distribution function F_XY is constant over any grid cell, including the left-hand and lower boundaries, at the value taken on at the lower left-hand corner of the cell. These lower left-hand corner values may be obtained systematically from the joint probability matrix P by a two step operation.
- Take cumulative sums upward of the columns of P.
- Take cumulative sums of the rows of the resultant matrix.
This can be done with the MATLAB function cumsum, which takes column cumulative sums downward. By flipping the matrix and transposing, we can achieve the desired results.
EXAMPLE 2: Calculation of F_XY values for Example 3 from "Random Vectors and Joint Distributions"
```
>> P = 0.1*[3 0 0; 0 6 0; 0 0 1];
>> FXY = flipud(cumsum(flipud(P)))  % Cumulative column sums upwardFXY =
    0.3000    0.6000    0.1000
         0    0.6000    0.1000
 Cumulative row sums
FXY =
   0         0    0.1000
>> FXY = cumsum(FXY')'              
% 0.3000    0.9000    1.0000
         0    0.6000    0.7000

0         0    0.1000
```
Figure 1: The joint distribution for Example 3 in "Random Vectors and Joint Distributions'.

Comparison with Example 3 from "Random Vectors and Joint Distributions" shows agreement with values obtained by hand.
The two step procedure has been incorprated into an m-procedure jddbn. As an example, return to the distribution in Example Example 1
EXAMPLE 3: Joint distribution function for Example 1
```
>> jddbn
Enter joint probability matrix (as on the plane)  PTo view joint distribution function, call for FXY
912    1.000>> disp(FXY)
    0.1512    0.3312    0.6012    0.
70
    0.1152    0.2754    0.5157    0.6848    0.8756
 0.0264    0.0534    0.0939    0.1224    0.1356
0.0780    0.1824    0.3390    0.4492    0.5656
  
 
```
These values may be put on a grid, in the same manner as in Figure 2 for Example 3 in "Random Vectors and Joint Distributions".
As in the case of canonic for a single random variable, it is often useful to have a function version of the procedure jcalc to provide the freedom to name the outputs conveniently.function[x,y,t,u,px,py,p] = jcalcf(X,Y,P) The quantities x,y,t,u,px,py, and p may be given any desired names.

Joint absolutely continuous random variables

In the single-variable case, the condition that there are no point mass concentrations on the line ensures the existence of a probability density function, useful in probability calculations. A similar situation exists for a joint distribution for two (or more) variables. For any joint mapping to the plane which assigns zero probability to each set with zero area (discrete points, line or curve segments, and countable unions of these) there is a density function.

Definition

If the joint probability distribution for the pair {X,Y{ assigns zero probability to every set of points with zero area, then there exists a joint density function f_XY with the property

P[(X,Y)∈Q]=∫∫_Qf_XY

(5)

We have three properties analogous to those for the single-variable case:

(f1)f_XY≥0(f2)∫∫_R²f_XY=1(f3)F_XY(t,u)=∫_−∞^t∫_−∞^uf_XY

(6)

At every continuity point for f_XY, the density is the second partial

f_XY(t,u)=

∂²F_XY(t,u)

∂t∂u

(7)

Now

F_X(t)=F_XY(t,∞)=∫_−∞^t∫_−∞^∞f_XY(r,s)dsdr

(8)

A similar expression holds for F_Y(u). Use of the fundamental theorem of calculus to obtain the derivatives gives the result

f_X(t)=∫_−∞^∞f_XY(t,s)dsandf_Y(u)=∫_−∞^∞f_XY(r,u)du

(9)

Marginal densities. Thus, to obtain the marginal density for the first variable, integrate out the second variable in the joint density, and similarly for the marginal for the second variable.

EXAMPLE 4: Marginal density functions

Let f_XY(t,u)=8tu0≤u≤t≤1. This region is the triangle bounded by u=0, u=t, and t=1 (see Figure 2)

f_X(t)=∫f_XY(t,u)du=8t∫₀^tudu=4t³,0≤t≤1

(10)

f_Y(u)=∫f_XY(t,u)dt=8u∫_u¹tdt=4u(1−u²),0≤u≤1

(11)

P(0.5≤X≤0.75,Y>0.5)=P[(X,Y)∈Q] where Q is the common part of the triangle with the strip between t=0.5 and t=0.75 and above the lineu=0.5. This is the small triangle bounded by u=0.5, u=t, and t=0.75. Thus

p=8∫_1/2^3/4∫_1/2^ttududt=25/256≈0.0977

(12)

**Figure 2:** Distribution for Example 4.

EXAMPLE 5: Marginal distribution with compound expression

The pair {X,Y{ has joint density f_XY(t,u)=

(t+2u) on the region bounded by t=0, t=2, u=0, and u=max{1,t{ (see Figure 3). Determine the marginal density f_X.

SOLUTION

Examination of the figure shows that we have different limits for the integral with respect to u for 0≤t≤1 and for 1<t≤2.

For 0≤t≤1
f_X(t)=

6

37

∫₀¹(t+2u)du=

6

37

(t+1)
(13)
For 1<t≤2
f_X(t)=

6

37

∫₀^t(t+2u)du=

12

37

t²
(14)

We may combine these into a single expression in a manner used extensively in subsequent treatments. Suppose M=[0,1] and N=(1,2]. Then I_M(t)=1 for t∈M(i.e., 0≤t≤1) and zero elsewhere. Likewise, I_N(t)=1 for t∈N and zero elsewhere. We can, therefore express f_X by

f_X(t)=I_M(t)

(t+1)+I_N(t)

t²

(15)

**Figure 3:** Marginal distribution for Example 5.

Discrete approximation in the continuous case

For a pair {X,Y{ with joint density f_XY, we approximate the distribution in a manner similar to that for a single random variable. We then utilize the techniques developed for a pair of simple random variables. If we have n approximating values t_i for X and m approximating values u_j for Y, we then have n·m pairs (t_i,u_j), corresponding to points on the plane. If we subdivide the horizontal axis for values of X, with constant increments dx, as in the single-variable case, and the vertical axis for values of Y, with constant increments dy, we have a grid structure consisting of rectangles of size dx·dy. We select t_i and u_j at the midpoint of its increment, so that the point (t_i,u_j) is at the midpoint of the rectangle. If we let the approximating pair be{X^*,Y^*{, we assign

p_ij=P((X^*,Y^*),=,(t_i,u_j))=P(X^*=t_i,Y^*=u_j)=P((X,Y)inijthrectangle)

(16)

As in the one-variable case, if the increments are small enough,

P((X,Y)∈ijthrectangle)≈dx·dy·f_XY(t_i,u_j)

(17)

The m-procedure tuappr calls for endpoints of intervals which include the ranges of X and Y and for the numbers of subintervals on each. It then prompts for an expression for f_XY(t,u), from which it determines the joint probability distribution. It calculates the marginal approximate distributions and sets up the calculating matrices t and u as does the m-process jcalc for simple random variables. Calculations are then carried out as for any joint simple pair.

EXAMPLE 6: Approximation to a joint continuous distribution

f_XY(t,u)=3on0≤u≤t²≤1

(18)

Determine P(X≤0.8,Y>0.1).

>> tuappr
Enter matrix [a b] of X-range endpoints  [0 1]Enter matrix [c d] of Y-range endpoints  [0 1]
er number of Y approximation points  200
EnEnter number of X approximation points  200
En
tter expression for joint density  3*(u <= t.^2)

>> M = (t <= 0.8)&(u > 0.1);
>> p = total(M.*P) Use array operations on X, Y, PX, PY, t, u, and 
P         % Evaluation of the integral with
3352455531
p =   0.3355                % Maple gives 0
.

The discrete approximation may be used to obtain approximate plots of marginal distribution and density functions.

**Figure 4:** Marginal density and distribution function for Example 7.

EXAMPLE 7: Approximate plots of marginal density and distribution functions

f_XY(t,u)=3u on the triangle bounded by u=0, u≤1+t, and u≤1−t.

>> tuappr
Enter matrix [a b] of X-range endpoints  [-1 1]Enter matrix [c d] of Y-range endpoints  [0 1]
er number of Y approximation points  200
EnEnter number of X approximation points  400
En
tter expression for joint density  3*u.*(u<=min(1+t,1-t))
PX/dx;                % Density for X  (see Use array operations on X, Y, PX, PY, t, u, and P
>> fx =
 Figure 4)
                              % Theoretical (3/2)(1 - |t|)^2>> fy = PY/dy;                % Density for Y
unction for X (Figure 4)>> FX = cumsum(PX);           % Distribution 
f
> FY = cumsum(PY);           % Distribution function for Y
>>> plot(X,fx,X,FX)            % Plotting details omitted