MAT 4375
Multivariate statistical methods
Summary of R commands (Chapter 4)
Univariate Normal Distribution
1. To generate a random sample of size n from a univariate normal distribution N (µ, σ 2 ) with mean µ and standard deviation σ (say n = 100,
µ = 0 and σ = 1):
rnorm(100,mean=0,sd=1)
or simply
rnorm(100,0,1)
2. To compute the value of the N (µ, σ 2 ) density function:
(x − µ)2
1
exp −
f (x) = √
2σ 2
2πσ 2
at point x: (say µ = 0, σ = 1 and x = 1.5)
dnorm(1.5,mean=0,sd=1)
or simply
dnorm(1.5,0,1)
3. To plot the graph of the density function of N (µ, σ 2 ) between the points
a and b (say we use again µ = 0, σ = 1 and the points a = −3 and
b = 3):
curve(dnorm(x,mean=0,sd=1),-3,3)
or simply
curve(dnorm(x,0,1),-3,3)
1
If we omit writing the values a and b, by default the graph will be
between a = 0 and b = 1.
4. To compute a normal probability, i.e to find the value P (X ≤ x) for
given x, where X ∼ N (µ, σ 2 ): (say µ = 0, σ = 1 and x = 1.96)
pnorm(1.96,mean=0,sd=1)
or simply
pnorm(1.96,0,1)
5. To compute a normal quantile, i.e. to find the value x such that
P (X ≤ x) = p for a given value p ∈ (0, 1), where X ∼ N (µ, σ 2 ):
(say µ = 0, σ = 1 and p = 0.975):
qnorm(0.975,mean=0,sd=1)
or simply
qnorm(0.975,0,1)
6. To do a normal QQ plot of univariate data (saved in a variable x),
there are two methods:
Method 1:
qqnorm(x)
abline(mean(x),sd(x))
# Do the QQ plot
# Plot the line of best fit
Method 2: (assume that n is the sample size)
p=ppoints(n)
# Compute the probabilities p_j=(j-3/8)/(n+1/4)
z=qnorm(p,0,1)
# Compute the normal quantiles corresponding to p_j
y=sort(x)
# Arrange the data in increasing order
plot(z,y)
# Plot the pairs (z_j,y_j) for j=1,...,n
abline(mean(x),sd(x)) #Plot the line of best fit
2
Multivariate Normal distribution
1. To generate a random sample of size n from a multivariate normal
distribution Np (µ, Σ) with mean vector µ and covariance matrix Σ.
Say p = 2, n = 100, µ = (0, 0)0 and
1 0.5
Σ=
0.5 1
We need to install the package MASS:
install.packages("MASS")
library(MASS)
mvrnorm(100,mu=c(0,0),Sigma=matrix(c(1,0.5,0.5,1),2))
In the last line, we could simply type:
mvrnorm(100,c(0,0),matrix(c(1,0.5,0.5,1),2))
2. To plot the Np (µ, Σ) density function:
f (x) =
1
(2π)p/2 |Σ|1/2
1
exp{− (x − µ)0 Σ−1 (x − µ)}
2
In the example below, we take p = 2, µ = (0, 0)0 and Σ as above.
#Specifies the parameters of the bivariate normal
mu1 = 0 # expected value of x
mu2 = 0.5 # expected value of y
sigma1 = 0.5 # variance of x
sigma2 = 2 # variance of y
rho = 0.5 # corr(x, y)
#Gives the range for x-axis and y-axis
xm =-3
xp = 3
ym = -3
yp = 3
3
#Specifies a sequence of values (x,y) for which we will do the plot
x = seq(xm, xp, length= as.integer((xp + abs(xm)) * 10))
y = seq(ym, yp, length= as.integer((yp + abs(ym)) * 10))
#Defines the density function
bivariate = function(x,y){
term1 = 1 / (2 * pi * sigma1 * sigma2 * sqrt(1 - rho^2))
term2 = (x - mu1)^2 / sigma1^2
term3 = -(2 * rho * (x - mu1)*(y - mu2))/(sigma1 * sigma2)
term4 = (y - mu2)^2 / sigma2^2
z = term2 + term3 + term4
term5 = term1 * exp((-z / (2 *(1 - rho^2))))
return (term5)
}
# Computes the density values
z =outer(x,y,bivariate)
# Plot
persp(x, y, z)
3. To compute a quantile corresponding to the χ2k distribution, i.e. to find
a value x such that P (X ≤ x) = p for a given value p ∈ (0, 1), where
X ∼ χ2k : (here k = 7 and p = 0.95)
qchisq(0.95,df=7)
or simply
qchisq(0.95,7)
4. Consider a variable x (which contains observations x1 , . . . , xn ) and a a
fixed value. To compute how many observations xj are smaller than a,
we first create a boolean variable b which takes values TRUE (i.e. 1)
if xj ≤ a and FALSE (i.e. 0) if xj > a, and then sum all the values of
b: (here we take a = 1.5)
b=x<1.5
sum(b)
# Create the boolean variable
# Count how many observations are smaller than 1.5
4
5. Suppose we have data consisting of n p-dimensional vectors x1 , . . . , xn
with sample mean x and sample covariance matrix S. We will assume
here that data is saved in a data frame called “table” whose columns
are called V1 , . . . , Vp . For simplicity we work only with the first two
variables in the table, so that p = 2.
a) To compute the squared Mahalanobis distances
d2j = (xj − x)S−1 (xj − x),
x1=table$V1
#
x2=table$V2
#
x=cbind(x1,x2)
#
m=colMeans(x)
#
S=cov(x)
#
d=mahalanobis(x,m,S)
print(d)
j = 1, . . . , n
Extract variable x1 from table
Extract variable x2 from table
Create a new variable x from x1 and x2
Compute the column means of x
Compute the covariance matrix of x
# Compute the squared Mahalanobis distances
# Print the values of d
b) To compute the number and the proportion of observations which
are smaller than χ22 (0.50) we use the following commands:
chisq.quantile=qchisq(0.50,df=2)
#Calculate the chi-square quantile
print(chisq.quantile)
number.obs=sum(d<chisq.quantile) #Calculate the number of observations
print(number.obs)
proportion=number.obs/n # Calculate the proportion of observations
print(proportion)
b) To see if the data in columns V1 and V2 has a bivariate normal
distribution, we use the following commands:
p=ppoints(n)
# Compute the probabilities p_j=(j-3/8)/(n+1/4)
q=qchisq(p,df=2) # Compute the chi-square quantiles corresponding to p_j
y=sort(d)
# Arrange d in increasing order
plot(q,d)
# Plot the pairs (q_j,d_j) for j=1,...,n
abline(0,1)
# Plot the line of best fit
5