Tutorial_1-R_solution.Rmd

---
title: "Applied Statistics with R (STAT 215)"
subtitle: "Tutorial 1: The R Language"
author: 
  - "Ziling Ma"
  - "ziling.ma@kaust.edu.sa"
date: "September 29, 2024"
output:
  prettydoc::html_pretty:
    theme: leonids
    highlight: vignette
    toc: true
    number_sections: true
---

<!-- <style type="text/css"> -->
<!-- body{ -->
<!--   font-size: 10pt; -->
<!-- } -->
<!-- </style> -->

# Assigning Objects {#Ex1}
a. Create an object that stores the value $32 × \frac{41}{8}$.

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
a <- 32 * 41 / 8
print(a)
# or you can simply use 'a'
```
</details>

b. Overwrite your object in (a) by itself divided by $2.33$. Print the result to the console.

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
# Notice here we would like to overwrite 
a <- a / 2.33
# check value of a 
a
```
</details>

c. Create a new object with the value $−8.2 × 10^{−3}$.

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
b <- -8.2 * 10^(-3)
```
</details>


d. Print directly to the console the result of multiplying (b) by (c).

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
(a * b)
```
</details>
<a href="#Ex1">Back to</a> <a href="#top">top</a>



# Vectors {#Ex2}
a. Create and store a sequence of values from 5 to −11 that progresses in steps of 0.3.


<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
sq <- seq(5,-11,by=-0.3)
sq
```
</details>

b. Overwrite the object from (a) using the same sequence with the order reversed.

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
sq <- rev(sq)
sq
```
</details>

c. Repeat the vector c(-1,3,-5,7,-9) twice, with each element repeated 10 times, and store the result. Display the result sorted from largest to smallest.


<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
# Create the vector
vec <- c(-1, 3, -5, 7, -9)

# Repeat the vector twice, with each element repeated 10 times
repeated_vec <- rep(vec, each = 10, times = 2)

# Sort the result from largest to smallest
sorted_vec <- sort(repeated_vec, decreasing = TRUE)

# Print the sorted result
sorted_vec
```
</details>

d. Extract the first and last elements of your vector from (b), storing them as a new object.


<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
# We name the vector 'sq'

# Extract the first element of 'sq', notice in R (also in Matlab), we start from 1 instead of 0 like in Python
sq[1]
# Extract the last element, same as in Python you use -1, in Matlab you use 'end'
sq[length(sq)]
```
</details>
<a href="#Ex1">Back to</a> <a href="#top">top</a>


# Matrix operations {#Ex3}
In `R`, use just one line of code to create the following mtrices and assign them to the variables `W`, and `X`, respectively.
\[
W=\begin{pmatrix}
1 & 4 & 7 \\
2 & 5 & 8 \\
3 & 6 & 9
\end{pmatrix};
X = \begin{pmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
7 & 8 & 9
\end{pmatrix};
\]

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
(W <- matrix(data = 1:9, nrow = 3, ncol = 3, byrow = FALSE))
```
```{r}
(X <- matrix(data = 1:9, nrow = 3, ncol = 3, byrow = TRUE))
```

</details>
Also, compute $W\odot X$ and $W\cdot X$, such that $\odot$ represents the element-wise multiplication and $\cdot$ represents matrix multiplication.
<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
W * X
```
```{r}
W %*% X
```
</details>

Given a matrix Y, where
```{r}
 Y <- matrix(c(2,3,1,4),nrow = 2)
 Y
```
How can you compute the inverse matrix of Y?

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
solve(Y)
```
</details>
<a href="#Ex1">Back to</a> <a href="#top">top</a>


# Matrix and functions {#Ex4}
Create a matrix $A\in\mathbb R^{3\times3}$, such that
\[
a_{ij} = i + j,\quad\forall i,j.
\]

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
A <- matrix(nrow = 3, ncol = 3)
for(i in 1:3) {
  for(j in 1:3) {
    A[i,j] <- i + j
  }
}
A
```
</details>

Also, create $B\in \mathbb R^{3\times3}$, such that
\[
b_{ij}=\begin{cases}
\text{TRUE}, & a_{ij} \text{ is even}, \\
\text{FALSE}, &\text{otherwise}.
\end{cases}
\]

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
fB <- function(b) return(!as.logical(b %% 2))
(B <- matrix(fB(A), nrow = 3, ncol = 3))
```
</details>

c. Create $C\in \mathbb R^{3\times3}$, such that
\[
c_{ij}=\begin{cases}
1, & a_{ij} \in \{1,2,3\}, \\
0, &\text{otherwise}.
\end{cases}
\]

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
fC <- function(c) return(as.numeric(c %in% c(1,2,3)))
(C <- apply(A, MARGIN = c(1,2), fC))
```
</details>

d. Confirm the dimensions of the matrix from (a) are $2 × 3$ if you remove any one row.

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
A_modified <- A[-1, ]
A_modified
dim(A_modified)
```
</details>
<a href="#Ex1">Back to</a> <a href="#top">top</a>

# Recursive functions 1 {#Ex5}
Write a function for computing $f(n)=n!,\ (n\in \mathbb N)$. Print the 10 first numbers from the sequence $f$.
Also, for $n=10$, compare your result with the `gamma()` function from `R` (recall that, for any positive integer $n$, $\Gamma(n) = (n-1)!$).

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
factorial <- function(n){
  if((n == 0) || (n == 1))
    return(1)
  else
    return (n * factorial(n-1))
}
sapply(0:9, FUN = factorial)
cat("Gamma(10) =", gamma(10), "; 9! =", factorial(9))
```
  
</details>
<a href="#Ex1">Back to</a> <a href="#top">top</a>


# Graph function {#Ex6}
Use the functions `outer` y `persp` to graph the following functions\n\n")
$$f(x, y) = 1 - exp(-1 / (x^2 + y^2)), \text{for}\: x, y ∈ [-5, 5]$$
$$f(x, y) = sin(x)cos(y), \text{for}\: x, y ∈ [-\pi/2, \pi/2]$$

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
# Define the range for x and y
x <- seq(-5, 5, length = 50)
y <- seq(-5, 5, length = 50)

# Define the function
f1 <- function(x, y) {
  1 - exp(-1 / (x^2 + y^2))
}

# Create a grid using outer
z1 <- outer(x, y, f1)

# Plot using persp
persp(x, y, z1, col = "salmon2", 
      xlab = "X", ylab = "Y", zlab = "f(x, y)",
      main = "Graph of f(x, y) = 1 - exp(-1/(x^2 + y^2))")

# Notice, you can adjust the theta and phi values in persp() if you want a different viewing angle, for example
persp(x, y, z1, theta = 30, phi = 30,col = "salmon2", 
      xlab = "X", ylab = "Y", zlab = "f(x, y)",
      main = "Graph of f(x, y) = 1 - exp(-1/(x^2 + y^2))")

```
  
</details>

<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
# Define the range for x and y
x2 <- seq(-pi/2, pi/2, length = 100)
y2 <- seq(-pi/2, pi/2, length = 100)

# Define the function
f2 <- function(x, y) {
  sin(x) * cos(y)
}

# Create a grid using outer
z2 <- outer(x2, y2, f2)

# Plot using persp
persp(x2, y2, z2, theta = 30, phi = 30, col = "lightgreen", 
      xlab = "X", ylab = "Y", zlab = "f(x, y)",
      main = "Graph of f(x, y) = sin(x) * cos(y)")

```
</details>

<a href="#Ex1">Back to</a> <a href="#top">top</a>




# Linear regression {#Ex7}
Suppose that we model the height of male teenagers as
\[
\mathrm{height}_i = \beta_0 + \beta_1 \cdot \mathrm{age}_i + \epsilon_i,\ \text{s.t. } \epsilon_i \sim \mathcal N\left(0, \sigma^2\right).
\]
Now, for a random sample of size $n=20$, we want to estimate $\boldsymbol \beta = (\beta_0,\beta_1)^\intercal$ as
\[
\boldsymbol \beta = (X^\intercal X)^{-1} X^{\intercal} \boldsymbol y,
\]
where
\[
X = \begin{pmatrix}
1 & \textrm{age}_1 \\
\vdots & \vdots \\
1 & \textrm{age}_{20}
\end{pmatrix},\quad
\boldsymbol y = \begin{pmatrix}
\textrm{height}_1 \\
\vdots\\
\textrm{height}_{20}
\end{pmatrix}.
\]
Compute $\hat {\boldsymbol \beta}$ and plot $\hat {\boldsymbol y} =X\hat {\boldsymbol \beta}$. For this exercise, generate the data set using the following code:
```{r}
set.seed(999)
beta.0 <- 140
beta.1 <- 2
ages <- sample(x = 12:18, size = 20, replace = TRUE)
heights <- round(x = beta.0 + beta.1 * ages + rnorm(n = 10, mean = 0, sd = 2), 
                 digits = 2)
data <- data.frame(heights = heights, ages = ages)
```
<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
X <- cbind(rep(x = 1, times = 20), data$ages)
y <- as.matrix(data$heights)
(beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y)
y_hat <- X %*% beta_hat
```
```{r}
plot(x = data$ages, y = data$heights, xlab = "Age", ylab = "Height (cm)")
par(new = TRUE)
plot(x = data$ages, y = y_hat, type = "b" , col = "red", axes = FALSE, xlab = "", ylab = "")
summary(lm(heights ~ ages, data))
```
</details>
<a href="#Ex1">Back to</a> <a href="#top">top</a>


# `R Markdown` and `R Sweave`

* Very elegent and powerful tools to use with `Rstudio`.
* Work very well for homework, reports, presentations, etc.
* `R Markdown` is mainly for generating `.html` files though it can also generate `.pdf` files and many more (this file is generated by `R Markdown`).
* `R Sweave` is like $\LaTeX$ but with `R` chunks within. It generates standard `.pdf` files.

# A fun exercise if you would like to try 

## Recursive functions
Create 3 `N`-element row vectors `p`, `q`, and `r`:
\[
\begin{cases}
\mathtt{p}: &\text{prime number sequence},\\
\mathtt{q}: &\text{Fibonacci sequence},\\
\mathtt{r}: &r_n=\lfloor\sqrt{n\pi}\rfloor.
\end{cases}
\]
Do NOT compute the terms manually (write functions!).
<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
prime_seq <- function(n, i = 1, primes = c()){
  if(length(primes) < n) {
    if (i == 2 || sum(as.numeric(i %% (1:i) == 0)) == 2){
      prime_seq(n = n, i = i + 1, primes = c(primes, i))
    }else{
      prime_seq(n = n, i = i + 1, primes = primes)
    }
  }else{
    return(primes)
  }
}

fibonacci_seq <- function (n){
  fib = c()
  for (i in 1:n){
    if (i == 1){
      fib <- c(fib, 0)
    }else if(i == 2){
      fib <- c(fib, 1)
    }else{
      fib <- c(fib, fib[i - 1] + fib[i - 2])
    }
  }
  return(fib)
}
r_seq <- function(n) return(floor(sqrt((1:n) * pi)))

N <- 5
(p <- prime_seq(n = N))
(q <- fibonacci_seq(n = N))
(r <- r_seq(N))
```
</details>
Also, create matrices:
\[
A=\begin{pmatrix}
\mathtt p \\ \mathtt q \\ \mathtt r
\end{pmatrix},
B=\begin{pmatrix}
\mathtt p^\intercal & \mathtt q^\intercal & \mathtt r^\intercal
\end{pmatrix}
\]
<details>
  <summary><span style="color:blue">Click for answer</span></summary>
```{r}
(A <- rbind(p, q, r))
(B <- cbind(p, q, r))
```
</details>
<a href="#Ex1">Back to</a> <a href="#top">top</a>