Variance

Summary: The mean value locates the center of the probability mass distribution induced by X on the real line. In this unit, we examine how expectation may be used for further characterization of the distribution for X. In particular, we deal with the variance and its square root the standard deviation. We identify some important properties of variance and identify the concept covariance. The variance is calculated for several distributions, comparing analytical and Matlab results.

In the treatment of the mathematical expection of a real random variable X, we note that the mean value locates the center of the probability mass distribution induced by X on the real line. In this unit, we examine how expectation may be used for further characterization of the distribution for X. In particular, we deal with the concept of variance and its square root the standard deviation. In subsequent units, we show how it may be used to characterize the distribution for a pair {X,Y{ considered jointly with the concepts covariance, and linear regression

Variance

Location of the center of mass for a distribution is important, but provides limited information. Two markedly different random variables may have the same mean value. It would be helpful to have a measure of the spread of the probability mass about the mean. Among the possibilities, the variance and its square root, the standard deviation, have been found particularly useful.

Definition. The variance of a random variable X is the mean square of its variation about the mean value:

Var[X]=σ_X²=E[(X−μ_X)²]whereμ_X=E[X](1)

The standard deviation for X is the positive square root σ_X of the variance.

Remarks

• If X(ω) is the observed value of X, its variation from the mean is X(ω)−μ_X. The variance is the probability weighted average of the square of these variations.
• The square of the error treats positive and negative variations alike, and it weights large variations more heavily than smaller ones.
• As in the case of mean value, the variance is a property of the distribution, rather than of the random variable.
• We show below that the standard deviation is a “natural” measure of the variation from the mean.
• In the treatment of mathematical expectation, we show that
  E[(X−c)²]isaminimumiffc=E[X],inwhichcaseE[(X−E[X])²]=E[X²]−E²[X](2)
  This shows that the mean value is the constant which best approximates the random variable, in the mean square sense.

Basic patterns for variance

Since variance is the expectation of a function of the random variable X, we utilize properties of expectation in computations. In addition, we find it expedient to identify several patterns for variance which are frequently useful in performing calculations. For one thing, while the variance is defined as E[(X−μ_X)²], this is usually not the most convenient form for computation. The result quoted above gives an alternate expression.

• (V1): Calculating formula. Var[X]=E[X²]−E²[X].
• (V2): Shift property. Var[X+b]=Var[X]. Adding a constant b to X shifts the distribution (hence its center of mass) by that amount. The variation of the shifted distribution about the shifted center of mass is the same as the variation of the original, unshifted distribution about the original center of mass.
• (V3): Change of scale. Var[aX]=a²Var[X]. Multiplication of X by constant a changes the scale by a factor |a|. The squares of the variations are multiplied by a². So also is the mean of the squares of the variations.
• (V4): Linear combinations
  1. Var[aX±bY]=a²Var[X]+b²Var[Y]±2ab(E,[,X,Y,],−,E,[,X,],E,[,Y,])
  2. More generally,
     Var[

     n

     ∑

     k=1

      ,a_k,X_k]=

     n

     ∑

     k=1

      a_k²Var[X_k]+2∑a_ia_j(E,[X_iX_j],−,E,[X_i],E,[X_j])(3)
  The term c_ij=E[X_iX_j]−E[X_i]E[X_j] is the covariance of the pair {X_i,X_j{, whose role we study in the unit on that topic. If the c_ij are all zero, we say the class is uncorrelated.

Remarks

• If the pair {X,Y{ is independent, it is uncorrelated. The converse is not true, as examples in the next section show.
• If the a_i=±1 and all pairs are uncorrelated, then
  Var[

  n

  ∑

  k=1

   ,a_i,X_i]=

  n

  ∑

  k=1

   Var[X_i](4)
  The variance add even if the coefficients are negative.

We calculate variances for some common distributions. Some details are omitted—usually details of algebraic manipulation or the straightforward evaluation of integrals. In some cases we use well known sums of infinite series or values of definite integrals. A number of pertinent facts are summarized in Appendix B. Some Mathematical Aids. The results below are included in the table in Appendix C.

Variances of some discrete distributions

1. Indicator function X=I_EP(E)=p,q=1−p E[X]=p
   E[X²]−E²[X]=E[I_E²]−p²=E[I_E]−p²=p−p²=p(1−p)=pq(5)
2. Simple random variableX=

   n

   ∑

   i=1

    t_iI_Ai (primitive form) P(A_i)=p_i.
   Var[X]=

   n

   ∑

   i=1

    t_i²p_iq_i−2∑t_it_jp_ip_j,sinceE[I_AiI_Aj]=0i≠j(6)
3. Binomial(n,p). X=

   n

   ∑

   i=1

    I_Eiwith{I_Ei:1≤i≤n{iidP(E_i)=p
   Var[X]=

   n

   ∑

   i=1

    Var[I_Ei]=

   n

   ∑

   i=1

    pq=npq(7)
4. Geometric(p). P(X=k)=pq^k∀k≥0E[X]=q/p 
   We use a trick: E[X²]=E[X(X−1)]+E[X]
   E[X²]=p

   ∞

   ∑

   k=0

    k(k−1)q^k+q/p=pq²

   ∞

   ∑

   k=2

    k(k−1)q^k−2+q/p=pq²

   2

   (1−q)3

    +q/p=2

   q2

   p2

    +q/p(8)
   Var[X]=2

   q2

   p2

    +q/p−(q/p)²=q/p²(9)
5. Poisson(μ)P(X=k)=e^−μ

   μk

   k!

    ∀k≥0 
   Using E[X²]=E[X(X−1)]+E[X], we have
   E[X²]=e^−μ

   ∞

   ∑

   k=2

    k(k−1)

   μk

   k!

    +μ=e^−μμ²

   ∞

   ∑

   k=2

    

   μk−2

   (k−2)!

    +μ=μ²+μ(10)
   Thus, Var[X]=μ²+μ−μ²=μ. Note that both the mean and the variance have common value μ.

Some absolutely continuous distributions

1. Uniform on (a,b)f_X(t)=

   1

   b−a

    a<t<bE[X]=

   a+b

   2

   E[X²]=

   1

   b−a

    ∫_a^bt²dt=

   b3−a3

   3(b−a)

    soVar[X]=

   b3−a3

   3(b−a)

    −

   (a+b)2

   4

    =

   (b−a)2

   12

   (11)
2. Symmetric triangular(a,b) Because of the shift property (V2), we may center the distribution at the origin. Then the distribution is symmetric triangular (−c,c), wherec=(b−a)/2. Because of the symmetry
   Var[X]=E[X²]=∫_−c^ct²f_X(t)dt=2∫₀^ct²f_X(t)dt(12)
   Now, in this case,
   f_X(t)=

   c−t

   c2

    0≤t≤csothatE[X²]=

   2

   c2

    ∫₀^c(ct²−t³)dt=

   c2

   6

    =

   (b−a)2

   24

   (13)
3. Exponential (λ) f_X(t)=λe^−λt,t≥0E[X]=1/λ
   E[X²]=∫₀^∞λt²e^−λtdt=

   2

   λ2

    sothatVar[X]=1/λ²(14)
4. Gamma(α,λ)f_X(t)=

   1

   Γ(α)

    λ^αt^α−1e^−λtt≥0E[X]=

   α

   λ

   E[X²]=

   1

   Γ(α)

    ∫₀^∞λ^αt^α+1e^−λtdt=

   Γ(α+2)

   λ2Γ(α)

    =

   α(α+1)

   λ2

   (15)
   Hence Var[X]=α/λ².
5. Normal(μ,σ²)E[X]=μ
   Consider Y∼N(0,1),E[Y]=0,Var[Y]=

   2

   \2π

    ∫₀^∞t²e^−t2/2dt=1.
   X=σY+μimpliesVar[X]=σ²Var[Y]=σ²(16)

Extensions of some previous examples

In the unit on expectations, we calculate the mean for a variety of cases. We revisit some of those examples and calculate the variances.

EXAMPLE 1: Expected winnings (Example 8 from "Mathematical Expectation: Simple Random Variables")

A bettor places three bets at $2.00 each. The first pays $10.00 with probability 0.15, the second $8.00 with probability 0.20, and the third $20.00 with probability 0.10.

SOLUTION

The net gain may be expressed

X=10I_A+8I_B+20I_C−6,withP(A)=0.15,P(B)=0.20,P(C)=0.10(17)

We may reasonbly suppose the class {A,B,C{ is independent (this assumption is not necessary in computing the mean). Then

Var[X]=10²P(A)[1−P(A)]+8²P(B)[1−P(B)]+20²P(C)[1−P(C)](18)

Calculation is straightforward. We may use MATLAB to perform the arithmetic.

c = [10 8 20];
p = 0.01*[15 20 10];
q = 1 - p;
VX = sum(c.^2.*p.*q)
VX =  58.9900

EXAMPLE 2: A function of X (Example 9 from "Mathematical Expectation: Simple Random Variables")

Suppose X in a primitive form is

X=−3I_C1−I_C2+2I_C3−3I_C4+4I_C5−I_C6+I_C7+2I_C8+3I_C9+2I_C10(19)

with probabilities P(C_i)=0.08,0.11,0.06,0.13,0.05,0.08,0.12,0.07,0.14,0.16.

Let g(t)=t²+2t. Determine E[g(X)] and Var[g(X)]

c = [-3 -1 2 -3 4 -1 1 2 3 2];            % Original coefficients
pc = 0.01*[8 11 6 13 5 8 12 7 14 16];     % Probabilities for C_j
G = c.^2 + 2*c                            % g(c_j)
EG = G*pc'                                % Direct calculation E[g(X)]
EG =  6.4200
VG = (G.^2)*pc' - EG^2                  % Direct calculation Var[g(X)]
VG = 40.8036
[Z,PZ] = csort(G,pc);                   % Distribution for Z = g(X)
EZ = Z*PZ'                              % E[Z]
EZ =  6.4200
VZ = (Z.^2)*PZ' - EZ^2                  % Var[Z]
VZ = 40.8036

EXAMPLE 3: Z=g(X,Y) (Example 10 from "Mathematical Expectation: Simple Random Variables")

We use the same joint distribution as for Example 10 from "Mathematical Expectation: Simple Random Variables" and let g(t,u)=t²+2tu−3u. To set up for calculations, we use jcalc.

jdemo1                      % Call for data
jcalc                       % Set up
Enter JOINT PROBABILITIES (as on the plane)  P
Enter row matrix of VALUES of X  X
Enter row matrix of VALUES of Y  Y
 Use array operations on matrices X, Y, PX, PY, t, u, and P
G = t.^2 + 2*t.*u - 3*u;    % Calculation of matrix of [g(t_i, u_j)]
EG = total(G.*P)            % Direct calculation of E[g(X,Y)]
EG =   3.2529
VG = total(G.^2.*P) - EG^2  % Direct calculation of Var[g(X,Y)]
VG =  80.2133
[Z,PZ] = csort(G,P);        % Determination of distribution for Z
EZ = Z*PZ'                  % E[Z] from distribution
EZ =   3.2529
VZ = (Z.^2)*PZ' - EZ^2      % Var[Z] from distribution
VZ =  80.2133

EXAMPLE 4: A function with compound definition (Example 12 from "Mathematical Expectation; General Random Variables")

Suppose X∼ exponential (0.3). Let

Z={

X2	forX≤4
16	forX>4

=I_[0,4](X)X²+I_(4,∞](X)16(20)

Determine E[Z] and Var[Z].

ANALYTIC SOLUTION

E[g(X)]=∫g(t)f_X(t)dt=∫₀^∞I_[0,4](t)t²0.3e^−0.3tdt+16E[I_(4,∞](X)](21)

=∫₀⁴t²0.3e^−0.3tdt+16P(X>4)≈7.4972(byMaple)(22)

Z²=I_[0,4](X)X⁴+I_(4,∞](X)256(23)

E[Z²]=∫₀^∞I_[0,4](t)t⁴0.3e^−0.3tdt+256E[I_(4,∞](X)]=∫₀⁴t⁴0.3e^−0.3tdt+256e^−1.2≈100.0562(24)

Var[Z]=E[Z²]−E²[Z]≈43.8486(byMaple)(25)

APPROXIMATION

To obtain a simple aproximation, we must approximate by a bounded random variable. Since P(X>50)=e⁻¹⁵≈3·10⁻⁷ we may safely truncate X at 50.

tappr
Enter matrix [a b] of x-range endpoints  [0 50]
Enter number of x approximation points  1000
Enter density as a function of t  0.3*exp(-0.3*t)
Use row matrices X and PX as in the simple case
M = X <= 4;
G = M.*X.^2 + 16*(1 - M);  % g(X)
EG = G*PX'                 % E[g(X)]
EG =  7.4972
VG = (G.^2)*PX' - EG^2     % Var[g(X)]
VG = 43.8472               % Theoretical = 43.8486
[Z,PZ] = csort(G,PX);      % Distribution for Z = g(X)
EZ = Z*PZ'                 % E[Z] from distribution
EZ =  7.4972
VZ = (Z.^2)*PZ' - EZ^2     % Var[Z]
VZ = 43.8472

EXAMPLE 5: Stocking for random demand (Example 13 from "Mathematical Expectation; General Random Variables")

The manager of a department store is planning for the holiday season. A certain item costs c dollars per unit and sells for p dollars per unit. If the demand exceeds the amount m ordered, additional units can be special ordered for s dollars per unit (s>c). If demand is less than the amount ordered, the remaining stock can be returned (or otherwise disposed of) at r dollars per unit (r<c). Demand D for the season is assumed to be a random variable with Poisson (μ) distribution. Supposeμ=50,c=30,p=50,s=40,r=20. What amount m should the manager order to maximize the expected profit?

PROBLEM FORMULATION

Suppose D is the demand and X is the profit. Then

• For D≤m,X=D(p−c)−(m−D)(c−r)=D(p−r)+m(r−c)
• For D>m,X=m(p−c)+(D−m)(p−s)=D(p−s)+m(s−c)

It is convenient to write the expression for X in terms of I_M, where M=(−∞,m]. Thus

X=I_M(D)[D(p−r)+m(r−c)]+[1−I_M(D)][D(p−s)+m(s−c)](26)

=D(p−s)+m(s−c)+I_M(D)[D(p−r)+m(r−c)−D(p−s)−m(s−c)](27)

=D(p−s)+m(s−c)+I_M(D)(s−r)[D−m](28)

Then

E[X]=(p−s)E[D]+m(s−c)+(s−r)E[I_M(D)D]−(s−r)mE[I_M(D)](29)

We use the discrete approximation.

APPROXIMATION

>> mu = 50;
>> n = 100;
>> t = 0:n;
>> pD = ipoisson(mu,t);         % Approximate distribution for D
>> c  = 30;
>> p  = 50;
>> s  = 40;
>> r  = 20;
>> m  = 45:55;
>> for i = 1:length(m)          % Step by step calculation for various m
    M = t<=m(i);
    G(i,:) = (p-s)*t + m(i)*(s-c) + (s-r)*M.*(t - m(i));
end
>> EG = G*pD';
>> VG = (G.^2)*pD' - EG.^2;
>> SG =sqrt(VG);
>> disp([EG';VG';SG']')
   1.0e+04 *
    0.0931    1.1561    0.0108
    0.0936    1.3117    0.0115
    0.0939    1.4869    0.0122
    0.0942    1.6799    0.0130
    0.0943    1.8880    0.0137
    0.0944    2.1075    0.0145
    0.0943    2.3343    0.0153
    0.0941    2.5637    0.0160
    0.0938    2.7908    0.0167
    0.0934    3.0112    0.0174
    0.0929    3.2206    0.0179

EXAMPLE 6: A jointly distributed pair (Example 14 from "Mathematical Expectation; General Random Variables")

Suppose the pair {X,Y{ has joint density f_XY(t,u)=3u on the triangular region bounded by u=0, u=1+t, u=1−t. Let Z=g(X,Y)=X²+2XY.

Determine E[Z] and Var[Z].

ANALYTIC SOLUTION

E[Z]=∫∫(t²+2tu)f_XY(t,u)dudt=3∫₋₁⁰∫₀^1+tu(t²+2tu)dudt+3∫₀¹∫₀^1−tu(t²+2tu)dudt=1/10(30)

E[Z²]=3∫₋₁⁰∫₀^1+tu(t²+2tu)²dudt+3∫₀¹∫₀^1−tu(t²+2tu)²dudt=3/35(31)

Var[Z]=E[Z²]−E²[Z]=53/700≈0.0757(32)

APPROXIMATION

tuappr
Enter matrix [a b] of X-range endpoints  [-1 1]
Enter matrix [c d] of Y-range endpoints  [0 1]
Enter number of X approximation points  400
Enter number of Y approximation points  200
Enter expression for joint density  3*u.*(u<=min(1+t,1-t))
Use array operations on X, Y, PX, PY, t, u, and P
G = t.^2 + 2*t.*u;          % g(X,Y) = X^2 + 2XY
EG = total(G.*P)            % E[g(X,Y)]
EG =   0.1006               % Theoretical value = 1/10
VG = total(G.^2.*P) - EG^2
VG =   0.0765               % Theoretical value 53/700 = 0.0757
[Z,PZ] = csort(G,P);        % Distribution for Z
EZ = Z*PZ'                  % E[Z] from distribution
EZ =  0.1006
VZ = Z.^2*PZ' - EZ^2
VZ =  0.0765

EXAMPLE 7: A function with compound definition (Example 15 from "Mathematical Expectation; General Random Variables")

The pair {X,Y{ has joint density f_XY(t,u)=1/2 on the square region bounded by u=1+t, u=1−t, u=3−t, and u=t−1.

W={

X	formax{X,Y{≤1
2Y	formax{X,Y{>1

=I_Q(X,Y)X+I_Q^c(X,Y)2Y(33)

where Q={(t,u):max{t,u{≤1{={(t,u):t≤1,u≤1{.

Determine E[W] and Var[W].

ANALYTIC SOLUTION

The intersection of the region Q and the square is the set for which 0≤t≤1 and 1−t≤u≤1. Reference to Figure 11.3.2 shows three regions of integration.

E[W]=

1

2

 ∫₀¹∫_1−t¹tdudt+

1

2

 ∫₀¹∫₁^1+t2ududt+

1

2

 ∫₁²∫_t−1^3−t2ududt=11/6≈1.8333(34)

E[W²]=

1

2

 ∫₀¹∫_1−t¹t²dudt+

1

2

 ∫₀¹∫₁^1+t4u²dudt+

1

2

 ∫₁²∫_t−1^3−t4u²dudt=103/24(35)

Var[W]=103/24−(11/6)²=67/72≈0.9306(36)

tuappr
Enter matrix [a b] of X-range endpoints  [0 2]
Enter matrix [c d] of Y-range endpoints  [0 2]
Enter number of X approximation points  200
Enter number of Y approximation points  200
Enter expression for joint density  ((u<=min(t+1,3-t))& ...
      (u$gt;=max(1-t,t-1)))/2
Use array operations on X, Y, PX, PY, t, u, and P
M = max(t,u)<=1;
G = t.*M + 2*u.*(1 - M);     % Z = g(X,Y)
EG = total(G.*P)              % E[g(X,Y)]
EG =  1.8340                  % Theoretical 11/6 = 1.8333
VG = total(G.^2.*P) - EG^2
VG =  0.9368                  % Theoretical 67/72 = 0.9306
[Z,PZ] = csort(G,P);          % Distribution for Z
EZ = Z*PZ'                    % E[Z] from distribution
EZ =  1.8340
VZ = (Z.^2)*PZ' - EZ^2
VZ =  0.9368

EXAMPLE 8: A function with compound definition

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

Z=I_Q(X,Y)X+I_Q^c(X,Y)forQ={(t,u):u+t≤1{(38)

The value t₀ where the line u=1−t and the curve u=t² meet satisfies t₀²=1−t₀.

E[Z]=3∫₀^t₀t∫₀^t2dudt+3∫_t0¹t∫₀^1−tdudt+3∫_t0¹∫_1−t^t2dudt=

3

4

 (5t₀−2)(39)

For E[Z²] replace t by t² in the integrands to get E[Z²]=(25t₀−1)/20.

Using t₀=(\5−1)/2≈0.6180, we get Var[Z]=(2125t₀−1309)/80≈0.0540.

APPROXIMATION

% Theoretical values
t0 = (sqrt(5) - 1)/2
t0 =  0.6180
EZ = (3/4)*(5*t0 -2)
EZ =  0.8176
EZ2 = (25*t0 - 1)/20
EZ2 = 0.7225
VZ = (2125*t0 - 1309)/80
VZ =  0.0540
tuappr
Enter matrix [a b] of X-range endpoints  [0 1]
Enter matrix [c d] of Y-range endpoints  [0 1]
Enter number of X approximation points  200
Enter number of Y approximation points  200
Enter expression for joint density  3*(u <= t.^2)
Use array operations on X, Y, t, u, and P
G = (t+u <= 1).*t + (t+u > 1);
EG = total(G.*P)
EG =  0.8169                   % Theoretical = 0.8176
VG = total(G.^2.*P) - EG^2
VG =  0.0540                   % Theoretical = 0.0540
[Z,PZ] = csort(G,P);
EZ = Z*PZ'
EZ =  0.8169
VZ = (Z.^2)*PZ' - EZ^2
VZ =  0.0540

Standard deviation and the Chebyshev inequality

In Example 5 from "Functions of a Random Variable," we show that if X∼N(μ,σ²) then Z=

X−μ

σ

 ∼N(0,1). Also, E[X]=μ and Var[X]=σ². Thus

P(

|X−μ|

σ

 ,≤,t)=P(|X−μ|≤tσ)=2Φ(t)−1(40)

For the normal distribution, the standard deviation σ seems to be a natural measure of the variation away from the mean.

For a general distribution with mean μ and variance σ², we have the

Chebyshev inequality

P(

|X−μ|

σ

 ,≥,a)≤

1

a2

 orP(|X−μ|≥aσ)≤

1

a2

(41)

In this general case, the standard deviation appears as a measure of the variation from the mean value. This inequality is useful in many theoretical applications as well as some practical ones. However, since it must hold for any distribution which has a variance, the bound is not a particularly tight. It may be instructive to compare the bound on the probability given by the Chebyshev inequality with the actual probability for the normal distribution.

t = 1:0.5:3;
p = 2*(1 - gaussian(0,1,t));
c = ones(1,length(t))./(t.^2);
r = c./p;
h = ['       t     Chebyshev   Prob     Ratio'];
m = [t;c;p;r]';
disp(h)
       t     Chebyshev   Prob     Ratio
disp(m)
    1.0000    1.0000    0.3173    3.1515
    1.5000    0.4444    0.1336    3.3263
    2.0000    0.2500    0.0455    5.4945
    2.5000    0.1600    0.0124   12.8831
    3.0000    0.1111    0.0027   41.1554

— □

DERIVATION OF THE CHEBYSHEV INEQUALITY

Let A={|X−μ|≥aσ{={(X−μ)²≥a²σ²{. Then a²σ²I_A≤(X−μ)².

Upon taking expectations of both sides and using monotonicity, we have

a²σ²P(A)≤E[(X−μ)²]=σ²(42)

from which the Chebyshev inequality follows immediately.

— □

We consider three concepts which are useful in many situations.

Definition. A random variable X is centered iff E[X]=0.

X^'=X−μisalwayscentered.(43)

Definition. A random variable X is standardized iff E[X]=0 and Var[X]=1.

X^*=

X−μ

σ

 =

X'

σ

 isstandardized(44)

Definition. A pair {X,Y{ of random variables is uncorrelated iff

E[XY]−E[X]E[Y]=0(45)

It is always possible to derive an uncorrelated pair as a function of a pair {X,Y{, both of which have finite variances. Consider

U=(X^*+Y^*)V=(X^*−Y^*),whereX^*=

X−μX

σX

 ,Y^*=

Y−μY

σY

(46)

Now E[U]=E[V]=0 and

E[UV]=E(X^*+Y^*)(X^*−Y^*)]=E[(X^*)²]−E[(Y^*)²]=1−1=0(47)

so the pair is uncorrelated.

EXAMPLE 9: Determining an uncorrelated pair

We use the distribution for Examples Example 10 from "Mathematical Expectation: Simple Random Variables" and Example 3, for which

E[XY]−E[X]E[Y]≠0(48)

jdemo1
jcalc
Enter JOINT PROBABILITIES (as on the plane)  P
Enter row matrix of VALUES of X  X
Enter row matrix of VALUES of Y  Y
 Use array operations on matrices X, Y, PX, PY, t, u, and P
EX = total(t.*P)
EX =   0.6420
EY = total(u.*P)
EY =   0.0783
EXY = total(t.*u.*P)
EXY = -0.1130
c = EXY - EX*EY
c  =  -0.1633                % {X,Y} not uncorrelated

VX = total(t.^2.*P) - EX^2
VX =  3.3016
VY = total(u.^2.*P) - EY^2
VY =  3.6566
SX = sqrt(VX)
SX =  1.8170
SY = sqrt(VY)
SY =  1.9122
x = (t - EX)/SX;            % Standardized random variables
y = (u - EY)/SY;
uu = x + y;                 % Uncorrelated random variables
vv = x - y;
EUV = total(uu.*vv.*P)      % Check for uncorrelated condition
EUV = 9.9755e-06            % Differs from zero because of roundoff