First we load the package:
library(matlib)EXAMPLE 1: Then we create the coefficient matrix and the vector for the right hand side of the system of linear equations.
A <- matrix(c(0, -1, 1, 0, 1, 1, 0, 1, 3, -4, 2, 0, -1, 0, 4, -4), 4, 4)
b <- c(1, 1, 5, -2)Here, we did not type “nrow=” and “ncol=”. If it is clear we can skip typying them. Using the showEqn() function we can see the system of linear equations.
showEqn(A, b)## 0*x1 + 1*x2 + 3*x3 - 1*x4 = 1
## -1*x1 + 1*x2 - 4*x3 + 0*x4 = 1
## 1*x1 + 0*x2 + 2*x3 + 4*x4 = 5
## 0*x1 + 1*x2 + 0*x3 - 4*x4 = -2Now using the echelon() function we can see how Gaussian Elimination works. With this example, we can type in R as:
echelon(A, b, verbose=TRUE, fractions=TRUE)##
## Initial matrix:
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 1 3 -1 1
## [2,] -1 1 -4 0 1
## [3,] 1 0 2 4 5
## [4,] 0 1 0 -4 -2
##
## row: 1
##
## exchange rows 1 and 2
## [,1] [,2] [,3] [,4] [,5]
## [1,] -1 1 -4 0 1
## [2,] 0 1 3 -1 1
## [3,] 1 0 2 4 5
## [4,] 0 1 0 -4 -2
##
## multiply row 1 by -1
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -1 4 0 -1
## [2,] 0 1 3 -1 1
## [3,] 1 0 2 4 5
## [4,] 0 1 0 -4 -2
##
## subtract row 1 from row 3
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -1 4 0 -1
## [2,] 0 1 3 -1 1
## [3,] 0 1 -2 4 6
## [4,] 0 1 0 -4 -2
##
## row: 2
##
## multiply row 2 by 1 and add to row 1
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 7 -1 0
## [2,] 0 1 3 -1 1
## [3,] 0 1 -2 4 6
## [4,] 0 1 0 -4 -2
##
## subtract row 2 from row 3
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 7 -1 0
## [2,] 0 1 3 -1 1
## [3,] 0 0 -5 5 5
## [4,] 0 1 0 -4 -2
##
## subtract row 2 from row 4
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 7 -1 0
## [2,] 0 1 3 -1 1
## [3,] 0 0 -5 5 5
## [4,] 0 0 -3 -3 -3
##
## row: 3
##
## multiply row 3 by -1/5
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 7 -1 0
## [2,] 0 1 3 -1 1
## [3,] 0 0 1 -1 -1
## [4,] 0 0 -3 -3 -3
##
## multiply row 3 by 7 and subtract from row 1
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 6 7
## [2,] 0 1 3 -1 1
## [3,] 0 0 1 -1 -1
## [4,] 0 0 -3 -3 -3
##
## multiply row 3 by 3 and subtract from row 2
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 6 7
## [2,] 0 1 0 2 4
## [3,] 0 0 1 -1 -1
## [4,] 0 0 -3 -3 -3
##
## multiply row 3 by 3 and add to row 4
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 6 7
## [2,] 0 1 0 2 4
## [3,] 0 0 1 -1 -1
## [4,] 0 0 0 -6 -6
##
## row: 4
##
## multiply row 4 by -1/6
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 6 7
## [2,] 0 1 0 2 4
## [3,] 0 0 1 -1 -1
## [4,] 0 0 0 1 1
##
## multiply row 4 by 6 and subtract from row 1
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 1
## [2,] 0 1 0 2 4
## [3,] 0 0 1 -1 -1
## [4,] 0 0 0 1 1
##
## multiply row 4 by 2 and subtract from row 2
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 1
## [2,] 0 1 0 0 2
## [3,] 0 0 1 -1 -1
## [4,] 0 0 0 1 1
##
## multiply row 4 by 1 and add to row 3
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 1
## [2,] 0 1 0 0 2
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 1Therefore the solution for this system is x1 = 1, x2 = 2, x3 = 0, x4 = 1
EXAMPLE 2:
A <- matrix(c(0, -9, 4, 4, 5, 6, 3, 8, 4, 6, 0, 4, 7, 7, 5, 12), 4, 4)
b <- c(-9, -17, -3, -12)Using the showEqn() function we can see the system of linear equations.
showEqn(A, b)## 0*x1 + 5*x2 + 4*x3 + 7*x4 = -9
## -9*x1 + 6*x2 + 6*x3 + 7*x4 = -17
## 4*x1 + 3*x2 + 0*x3 + 5*x4 = -3
## 4*x1 + 8*x2 + 4*x3 + 12*x4 = -12Now using the echelon() function we obtain the reduced echelon form of the augmented matrix. From this we have the reduced echelon form of the augmented matrix as follows:
echelon(A, b, verbose=TRUE, fractions=TRUE)##
## Initial matrix:
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 5 4 7 -9
## [2,] -9 6 6 7 -17
## [3,] 4 3 0 5 -3
## [4,] 4 8 4 12 -12
##
## row: 1
##
## exchange rows 1 and 2
## [,1] [,2] [,3] [,4] [,5]
## [1,] -9 6 6 7 -17
## [2,] 0 5 4 7 -9
## [3,] 4 3 0 5 -3
## [4,] 4 8 4 12 -12
##
## multiply row 1 by -1/9
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -2/3 -2/3 -7/9 17/9
## [2,] 0 5 4 7 -9
## [3,] 4 3 0 5 -3
## [4,] 4 8 4 12 -12
##
## multiply row 1 by 4 and subtract from row 3
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -2/3 -2/3 -7/9 17/9
## [2,] 0 5 4 7 -9
## [3,] 0 17/3 8/3 73/9 -95/9
## [4,] 4 8 4 12 -12
##
## multiply row 1 by 4 and subtract from row 4
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -2/3 -2/3 -7/9 17/9
## [2,] 0 5 4 7 -9
## [3,] 0 17/3 8/3 73/9 -95/9
## [4,] 0 32/3 20/3 136/9 -176/9
##
## row: 2
##
## exchange rows 2 and 4
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -2/3 -2/3 -7/9 17/9
## [2,] 0 32/3 20/3 136/9 -176/9
## [3,] 0 17/3 8/3 73/9 -95/9
## [4,] 0 5 4 7 -9
##
## multiply row 2 by 3/32
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 -2/3 -2/3 -7/9 17/9
## [2,] 0 1 5/8 17/12 -11/6
## [3,] 0 17/3 8/3 73/9 -95/9
## [4,] 0 5 4 7 -9
##
## multiply row 2 by 2/3 and add to row 1
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 -1/4 1/6 2/3
## [2,] 0 1 5/8 17/12 -11/6
## [3,] 0 17/3 8/3 73/9 -95/9
## [4,] 0 5 4 7 -9
##
## multiply row 2 by 17/3 and subtract from row 3
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 -1/4 1/6 2/3
## [2,] 0 1 5/8 17/12 -11/6
## [3,] 0 0 -7/8 1/12 -1/6
## [4,] 0 5 4 7 -9
##
## multiply row 2 by 5 and subtract from row 4
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 -1/4 1/6 2/3
## [2,] 0 1 5/8 17/12 -11/6
## [3,] 0 0 -7/8 1/12 -1/6
## [4,] 0 0 7/8 -1/12 1/6
##
## row: 3
##
## exchange rows 3 and 4
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 -1/4 1/6 2/3
## [2,] 0 1 5/8 17/12 -11/6
## [3,] 0 0 7/8 -1/12 1/6
## [4,] 0 0 -7/8 1/12 -1/6
##
## multiply row 3 by 8/7
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 -1/4 1/6 2/3
## [2,] 0 1 5/8 17/12 -11/6
## [3,] 0 0 1 -2/21 4/21
## [4,] 0 0 -7/8 1/12 -1/6
##
## multiply row 3 by 1/4 and add to row 1
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 1/7 5/7
## [2,] 0 1 5/8 17/12 -11/6
## [3,] 0 0 1 -2/21 4/21
## [4,] 0 0 -7/8 1/12 -1/6
##
## multiply row 3 by 5/8 and subtract from row 2
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 1/7 5/7
## [2,] 0 1 0 31/21 -41/21
## [3,] 0 0 1 -2/21 4/21
## [4,] 0 0 -7/8 1/12 -1/6
##
## multiply row 3 by 7/8 and add to row 4
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 1/7 5/7
## [2,] 0 1 0 31/21 -41/21
## [3,] 0 0 1 -2/21 4/21
## [4,] 0 0 0 0 0
##
## row: 4=====PRACTICAL APPLICATION===== For this practical application, we go back to the example of a contingency table shown in Table 1.6 in the beginning of this section. The table was constructed by the data from the 1996 General Society Survey, National Opinion Research Center. Recall that the measurements in Table 1.6 are income in terms of the US dollars and job satisfaction among black men in the US. Let us assign the numbers to the status of satisfactory. Let us assign the status “Very Dissatisfied” as 1, the status “Little Dissatisfied” as 2, the status “Moderately Satisfied” as 3 and the status “Very Satisfied” as 4. Similarly, for the income, let us assign the income < 15, 000 as 1, the income 15, 000 − 25, 000 as 2, the income 25, 000 − 40, 000 as 3 and the income > 40, 000 as 4. From the row and column sums computed from Table 1.6, we have a 4 × 4 table such that where the variable xij is the number of people who answered the questions for the income as i = 1, 2, 3, 4 and for the satisfactory status as j = 1, 2, 3, 4. As we discussed before, in order to make sure that we can protect personal information, we perturb the column sums (4, 13, 43, 63) as (4.5, 14.2, 42.4, 62.9) and the row sums (20, 22, 33, 21) as (21.5, 20.9, 33.4, 22.1). Now the question is whether there exists a 4 × 4 table satisfying such row and column sums: This can be written as a system of linear equations such that
A <- matrix(c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0,
0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,
1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
0, 0, 0, 1 ), nrow = 8, ncol = 16, byrow = TRUE)which is a 8 × 16 matrix and the vector for the right hand side of the system of linear equations is
b <- c(21.5, 20.9, 33.4, 22.1, 4.5, 14.2, 42.4, 62.9)Then we use the Solve() function. The following is the output from R:
Solve(A,b)## x1 - 1*x6 - 1*x7 - 1*x8 - 1*x10 - 1*x11 - 1*x12 - 1*x14 - 1*x15 - 1*x16 = -71.9
## x2 + x6 + x10 + x14 = 14.2
## x3 + x7 + x11 + x15 = 42.4
## x4 + x8 + x12 + x16 = 36.8
## x5 + x6 + x7 + x8 = 20.9
## x9 + x10 + x11 + x12 = 33.4
## x13 + x14 + x15 + x16 = 22.1
## 0 = 26.1When we look at the last equation, it says 0 = 26.1, which does not make sense. This is the form in Theorem 1.6. Therefore, there is no solution to the system of linear equations. Thus, we conclude that there does not exist such a table. This means that this perturbed information does not make sense for a statistical analysis for studying relationship. Thus perturbation does not work for releasing this information.