R for Statistics and Math_Part 3

1. Lợi thế của R

R có một lợi thế so với các ngôn ngữ khác đó là khả năng biểu diễn trực quan dữ liệu rất mạnh, hiện tại ít có ngôn ngữ lập trình nào được tích hợp nhiều package hỗ trợ biểu đồ như R, kể đến trong số đó là \(ggplot2\). Đây là package hỗ trợ vẽ biểu đồ trong R. Một điểm mạnh nữa của R là khả năng hỗ trợ thống kê và Toán học. Bài viết thứ 3 này sẽ cung cấp một số nội dung liên quan đến học phần Toán Cao Cấp (Advanced Math), một môn học cơ bản và rất quan trọng đối với các bạn sinh viên năm thứ nhất các trường Đại học về kỹ thuật, công nghệ, kinh tế….

2. Sử dụng R cho học phần Toán Cao Cấp

2.1. Một số phép toán cơ bản

Chúng ta cùng điểm qua các phép toán thông dụng, đó là, phép cộng \(+\), phép trừ \(-\), phép nhân \(*\) và phép chia \(/\). Ngoài ra R có hỗ trợ tất cả các hàm thông dụng khác như hàm logarit, hàm \(sin\), \(cos\)….. Xét một số ví dụ sau:

-Tính: \(a=(4^5-1/6)\)(\(\sqrt[3]{e}\)+\(\pi\))$

-Tính: \(\log_{3}(4)\)

-Tính: \(arcsin(1/2)\)

-Cho: \(f(x)= (sin(x)+x^{2})/(e^{x}+1)\). Hãy tính \(f(\pi/6)\)

a = (4^5-1/6)*(exp(1/3)+pi)
log(4)/log(3) #Sử dụng công thức đổi cơ số

## [1] 1.26186

asin(1/2)

## [1] 0.5235988

(sin(pi/6)+(pi/6)^2)/(exp(pi/6)+1)

## [1] 0.2879945

2.2. Ma trận và định thức

Như đã trình bày ở phần trước, ma trận và định thức là hai khái niệm quan trọng trong đại số tuyến tính. Tại phần này chúng ta sẽ cùng thực hiện một số lệnh cụ thể và ứng dụng với R cho tính toán ma trận. Phân rã \(QR\) - sử dụng để tính hạng của ma trận

Trong đại số tuyến tính, phân rã QR, còn được gọi là phân tích nhân tố \(QR\) hoặc phân tích nhân tố \(QU\) là phân rã ma trận \(A\) thành tích \(A = QR\) của ma trận trực giao \(Q\) và ma trận tam giác trên R. Phân rã QR thường được sử dụng để giải quyết vấn đề bình phương tối thiểu tuyến tính và là cơ sở cho một thuật toán \(eigenvalue\) cụ thể, thuật toán \(QR\).

#Tạo ra 2 ma trận B và C
y1 <- c(11,4,20, 4, 9, 8,3,6,9)
B <- matrix(y1, nrow=3)
y2 <- c(0,4,17,-2,5,8,3.5,8,-9.2)
C <- matrix(y2, nrow=3)
#Cộng, trừ hai ma trận
B+C

##      [,1] [,2] [,3]
## [1,]   11    2  6.5
## [2,]    8   14 14.0
## [3,]   37   16 -0.2

B-C

##      [,1] [,2] [,3]
## [1,]   11    6 -0.5
## [2,]    0    4 -2.0
## [3,]    3    0 18.2

#Phép nhân hai ma trận
B%*%C

##      [,1] [,2] [,3]
## [1,]   67   22 42.9
## [2,]  138   85 30.8
## [3,]  185   72 51.2

#Ma trận chuyển vị của ma trận B
t(B)

##      [,1] [,2] [,3]
## [1,]   11    4   20
## [2,]    4    9    8
## [3,]    3    6    9

#Tìm hạng của ma trận B dựa vào phân rã qr(matrix_name)$rank
qr(B)$rank

## [1] 3

#Tìm ma trận nghịch đảo của ma trận B
solve(B)

##            [,1]        [,2]        [,3]
## [1,]  0.1294118 -0.04705882 -0.01176471
## [2,]  0.3294118  0.15294118 -0.21176471
## [3,] -0.5803922 -0.03137255  0.32549020

#Tìm định thức của B(B là ma trận vuông)
det(B)

## [1] 255

2.3. Hệ phương trình tuyến tính

Ví dụ 1: Tìm nghiệm riêng của hệ phương trình tuyến tính \[\begin{cases} x_{1} + x_{2} + x_{3} = 6\\ x_{1} - x_{2} = -1 \\ x_{1} + x_{2} +2x_{3} = 9 \end{cases}\] Ví dụ 2: Tìm nghiệm riêng của hệ phương trình tuyến tính \[\begin{cases} x_{1} + x_{2} - x{3} = 0\\ 3x_{1} - x_{3} = 3 \end{cases}\] Ví dụ 3: Tìm nghiệm riêng của phương trình: \[\begin{cases} x_{1} + x_{2} - x_{3} = 3 \end{cases}\]

#Giải hệ phương trình tuyến tính, tìm nghiệm riêng Ví dụ 1:
A <- matrix(c(1,1,1,1,-1,1,1,0,2), nrow=3)
y <- matrix(c(6,-1,9),nrow =3)
lm.fit(A, y)$coefficients #Chúng ta thu được nghiệm duy nhất x=(1,2,3)

## x1 x2 x3 
##  1  2  3

#Giải hệ phương trình tuyến tính Ví dụ 2:
A1 <- matrix(c(1,3,1,0,-1,-1), nrow=2)
y1 <- matrix(c(0,3), nrow=2)
lm.fit(A1, y1)$coefficients#Trong trường hợp này,  nghiệm của hệ có 1 tham số, nghiệm riêng là: x=(1,-1,0)

## x1 x2 x3 
##  1 -1 NA

#Giải hệ phương trình ví dụ 3:
A2 <- matrix(c(1,1,-1), nrow = 1)
y2 <- matrix(3, nrow = 1)
lm.fit(A2, y2)$coefficients

## x1 x2 x3 
##  3 NA NA

#Trường hợp này nghiệm của hệ có 2 tham số, nghiệm riêng của hệ là:
#x=(3,0,0)

2.4. Vẽ đồ thị hàm số trong R

Vẽ đồ thị các hàm số một biến hoặc nhiều biến số trong Toán học là phương pháp biểu diễn hình học các giá trị của hàm số. Chúng ta cùng điểm qua một số hàm số thông dụng và cách vẽ các hàm số này trong R

2.4.1. Vẽ đồ thị hàm số một biến (2D)

Giả sử chúng ta cần vẽ đồ thị hàm số: \(y=x^{2}+1\) trên đoạn\([-10;10]\) Thực hiện như sau:

x <- seq(-10,10, length.out=200) # Tạo ra 200 điểm trong khoảng [-10;10]
y <- (x^2 + 1)
plot(y, type = "l", col = "red", main = "The graph of f(x)", xlab = "x value", ylab="f(x) value")

#Đồ thị trong R có chất lượng cao, chúng ta hoàn toàn có thể save lại graph này và đưa vào các bản báo cáo.

2.4.2. Vẽ đồ thị hàm 2 biến số (3D)

Ví dụ: Vẽ đồ thị hàm số \(z=sin(y)*y\) với \(x\in[0;2\pi], y\in[0;5]\) Tuy nhiên, trước tiên chúng ta cần biết cách sử dụng hàm \(outer()\) trong R. \(outer()\) là một hàm (function) hay đôi khi gọi là phương thức, sử dụng cho các phép tính toán trên 2 véc tơ hoặc 2 mảng một chiều. Phương thức \(persp()\) trong R - Perspective plots, vẽ các khối graphs đa chiều, sử dụng gồm 3 tham số truyền vào là các véc tơ và 1 mảng gồm các véc tơ. Trong đó các tham số bổ sung như :

• \(theta\) và \(phi\) cho phép xác định góc nhìn với \(theta\) dành cho hướng phương vị và \(phi\) là góc nhìn. Thông thường chúng ta chọn theta và phi bằng 30.
• \(expand\) là tham số phụ cho phép kéo hoặc dãn biểu đồ theo \(scale\). Thường cho \(expand =0.5\)
• \(col\) cho phép xác định màu của đồ thị
• \(ticktype\) có hai giá trị:
  • Nếu \(ticktype="simple"\) thì vẽ một mũi tên song song với các trục
  • Nếu \(ticktype="detailed"\) thì hiển thị các giá trị trên các trục
• \(shade\) cho phép tạo độ bóng và độ râm cho đồ thị, với các giá trị nằm trong khoảng [0.5, 0.75] cho phép tạo hiệu ứng giống như ánh sáng ban ngày. Chúng ta cùng xem một số ví dụ sau với hàm \(outer()\) và hàm \(persp()\) Thực hiện như sau:

# Sử dụng hàm outer() với một véc tơ số và 1 giá trị
x1 <- 1:5               # Create x vector
y1 <- 3                 # Create y value
outer(x1, y1, "+") -> output1
output1

##      [,1]
## [1,]    4
## [2,]    5
## [3,]    6
## [4,]    7
## [5,]    8

# Sử dụng hàm outer() cho 2 véc tơ số
x2 <- 1:5   # Create x vector (same as in Example 1)
y2 <- 3:7   # Create y vector
names(x2) <- x2         # Prepare row names of output matrix
names(y2) <- y2         # Prepare column names of output matrix
output2 <- outer(x2, y2, "+")       # Apply R outer function
output2       # Print output to RStudio console

##   3 4  5  6  7
## 1 4 5  6  7  8
## 2 5 6  7  8  9
## 3 6 7  8  9 10
## 4 7 8  9 10 11
## 5 8 9 10 11 12

# Sử dụng hàm outer() cho các dạng dữ liệu khác nhau

x3 <- letters[1:5]        # Create x vector with characters
y3 <- LETTERS[1:5]        # Create y vector with upper case
names(x3) <- x3         # Prepare row names of output matrix
names(y3) <- y3       # Prepare column names of output matrix
output3 <- outer(x3, y3, "paste") # Apply R outer function to characters
output3                            # Print output to RStudio console

##   A     B     C     D     E    
## a "a A" "a B" "a C" "a D" "a E"
## b "b A" "b B" "b C" "b D" "b E"
## c "c A" "c B" "c C" "c D" "c E"
## d "d A" "d B" "d C" "d D" "d E"
## e "e A" "e B" "e C" "e D" "e E"

# Sử dụng hàm outer() cho các function khác (custom function)
x4 <- 1:5                   # Create x vector (same as in Example 1 & 2)
y4 <- 3:7                 # Create y vector (same as in Example 2)
custom_fun <- function(x, y) {  # Create custom function in R
  z <- (x + y)^2
  return(z)
}
names(x4) <- x4           # Prepare row names of output matrix
names(y4) <- y4         # Prepare column names of output matrix

output4 <- outer(x4, y4, custom_fun)
# Apply R outer function to custom function
output4     # Print output to RStudio console

##    3  4   5   6   7
## 1 16 25  36  49  64
## 2 25 36  49  64  81
## 3 36 49  64  81 100
## 4 49 64  81 100 121
## 5 64 81 100 121 144

########################################################
# Để sử dụng comment nhiều code trong R, bôi đen code cần comment và gõ Ctrl+Shift+C
x <- seq(0, 2*pi, length = 30)
y <- seq(0,5, length = 30)
f <- function(x,y){
  r <- sin(x)*y}
z <- outer(x, y, f)
persp(x,y,z, theta = 30, phi =30, expand=0.5, col="lightblue",
      ltheta = 120, shade = 0.5, ticktype = "simple",
      xlab = "Ox axis", ylab = "Oy axis", main = "Plot 3D in R")

# Ví dụ khác dùng persp ()
x <- seq(-10, 10, length=40)
y <- seq(-10, 10, length=40)

z <- matrix (
  outer(x, y,
        function(x, y) {
          14 - 
          sqrt( (x+cos(y)/10) ^2 + (y - sin(x)/5) ^2) +
          cos((x-2)/3)*sin(x) + 
          cos((y-2)/4)*cos(y) +
        # 
          sqrt( (x+3) ^ 2 + (y/2-2) ^2)
        }
  ),
  nrow=length(x)
)
persp(x, y, z)