# fahrenheit degrees
far_deg <- 100
# convert to celsius
(far_deg - 32) * (5/9)[1] 37.77778Practice: Functions (part 1)
Stat 133
Learning Objectives
- Define a function that takes arguments
- Return a value from a function
- Test a function
- Set default values for function arguments
The goal of this module is to give you practice writing very simple functions. The following exercises are mostly designed for those of you without any programming experience. If you have previous programming experience, you can move to the next module.
1 Function Basics
To define a new function in R you use the function function(). Actually, you need to specify a name for the function, and then assign function() to the chosen name. You also need to define optional arguments (i.e. inputs). And of course, you must write the code (i.e. the body) so the function does something when you use it:
# anatomy of a function
some_name <- function(arguments) {
# body of the function
}2 From Fahrenheit to Celsius
Let’s consider a typical programming example that involves converting fahrenheit degrees into celsius degrees. The conversion formula is \((F - 32) \times 5/9 = C\). Here’s some R code to convert 100 fahrenheit degrees into Celsius degrees:
What if you want to convert 90 fahrenheit degrees in Celsius degrees? One option would be to rewrite the previous lines as:
# fahrenheit degrees
far_deg <- 90
# convert to celsius
(far_deg - 32) * (5/9)[1] 32.22222However, retyping many lines of code can be very boring, tedious, and inefficient. To make your code reusable in a more efficient manner, you will have to write functions.
2.1 Writing a simple function
So, how do you create a function? The first step is to write code and make sure that it works. In this case we already have the code that converts a number in Fahrenheit units into Celsius.
The next step is to encapsulate the code in the form of a function. You have to choose a name, some argument(s), and determine the output. Here’s one example with a function fahrenheit_to_celsius()
fahrenheit_to_celsius <- function(x) {
(x - 32) * (5/9)
}
fahrenheit_to_celsius(100)[1] 37.77778If you want to get the conversion of 90 fahrenheit degrees, you just simply execute it again by changing its argument:
fahrenheit_to_celsius(90)[1] 32.22222And because we are using arithmetic operators (i.e. multiplication, subtraction, division), the function is also vectorized:
fahrenheit_to_celsius(c(90, 100, 110))[1] 32.22222 37.77778 43.33333Sometimes it is recommended to add a default value to one (or more) of the arguments. In this case, we can give a default value of x = 1. When the user executes the function without any input, fahrenheit_to_celsius returns the value of 1 fahrenheit degree to Celsius degrees:
fahrenheit_to_celsius <- function(x = 1) {
(x - 32) * (5/9)
}
# default execution
fahrenheit_to_celsius()[1] -17.222223 Unit Conversion Functions
3.1 miles to kilometers
Write a function miles_to_kms() that converts miles into kilometers: 1 mile is equal to 1.6 kilometers. Give the argument a default value of 1.
Show answer
miles_to_kms = function(x = 1) {
x * 1.6
}3.2 gallons to liters
Write a function gallon_to_liters() that converts gallons to liters: 1 gallon is equal to 3.78541 liters:
Show answer
gallon_to_liters = function(x = 1) {
x * 3.78541
}3.3 seconds to years
According to Wikipedia, in 2015 the life expectancy of a person born in the US was 79 years. Consider the following question: Can a newborn baby in USA expect to live for one billion (\(10^9\)) seconds?
To answer this question, write a function seconds2years() that takes a number in seconds and returns the equivalent number of years. Assume a year with 365 days. Test the function with seconds2years(1000000000)
Show answer
seconds2years = function(x = 1) {
x / (365 * 24 * 60 * 60)
}
seconds2years(1000000000)4 Simple Math Functions
Consider the following mathematical functions:
\[ f(x) = x^2, \quad g(x) = 2x + 5 \]
4.1 Function f()
Write a function f() based on the above equation.
Show answer
f = function(x = 1) {
x^2
}Test your function with:
f(2)
f(-5)4.2 Function g()
Write a function g() based on the above equation.
Show answer
g = function(x = 1) {
2*x + 5
}Test your function with:
g(0)
g(-5/2)4.3 Composite function fog()
Write code to create the composite function fog(): \(f \circ g(x)\)
Show answer
fog = function(x = 1) {
f(g(x))
}Test your composite function with:
fog(2)
fog(-5)4.4 Composite function gof()
Write code to create the composite function gof(): \(g \circ f(x)\)
Show answer
gof = function(x = 1) {
g(f(x))
}Test your composite function with:
gof(0)
gof(-5/2)5 Polynomials
In this problem we want to see whether the graph of a given polynomial will cross or touch the x axis in a given interval.
Let’s begin with the polynomial: \(f(x) = x^2 (x -1)\). The first thing to do is write a function for the polynomial, for instance:
poly1 <- function(x) {
(x^2) * (x - 1)
}Once you have a function for the polynomial, you can create a set of pairs of points \(x\) and \(y = f(x)\), and then use them to graph the polynomial
# set of points
x <- seq(-4, 4, length.out = 30)
y <- poly1(x)
# graph polynomial
plot(x, y, type = 'l', lwd = 3, col = "#FB7215", las = 1)
abline(h = 0, v = 0, col = '#888888aa', lwd = 1.5)
title(main = expression(paste(f(x), ' = ', x^2, (x - 1))))
5.1 Some polynomials
Write functions and graph the following polynomials in the x-axis interval -4 to 4:
- \(f(x) = x^3\)
Show answer
poly2 <- function(x) {
x^3
}
x <- seq(-4, 4, length.out = 30)
y <- poly2(x)
# graph polynomial
plot(x, y, type = 'l', lwd = 3, col = "#FB7215", las = 1)
abline(h = 0, v = 0, col = '#888888aa', lwd = 1.5)
title(main = expression(paste(f(x), ' = ', x^3)))- \(f(x) = (x^2 - 1)(x + 3)^3\)
Show answer
poly3 <- function(x) {
(x^2 - 1) * (x + 3)^3
}
x <- seq(-4, 4, length.out = 30)
y <- poly3(x)
# graph polynomial
plot(x, y, type = 'l', lwd = 3, col = "#FB7215", las = 1)
abline(h = 0, v = 0, col = '#888888aa', lwd = 1.5)
title(main = expression(paste(f(x), ' = ', (x^2 - 1), (x + 3)^3)))- \(f(x) = (x^2 - 1)(x^2 - 9)\)
Show answer
poly4 <- function(x) {
(x^2 - 1) * (x^2 - 9)
}
x <- seq(-4, 4, length.out = 30)
y <- poly4(x)
# graph polynomial
plot(x, y, type = 'l', lwd = 3, col = "#FB7215", las = 1)
abline(h = 0, v = 0, col = '#888888aa', lwd = 1.5)
title(main = expression(paste(f(x), ' = ', (x^2 - 1), (x^2 - 9))))6 Pythagoras
The Pythagoras formula is used to compute the length of the hypotenuse, \(h\), of a right triangle with legs of length \(a\) and \(b\).
\[ h = \sqrt{(a^2 + b^2)} \]
- Write a function
pythagoras()that takes two argumentsaandb, and returns the length of the hypotenuse. Test yourpythagoras()with two leg values:pythagoras(3, 4)
Show answer
pythagoras <- function(a, b) {
sqrt(a^2 + b^2)
}- Rewrite your function
pythagoras()in a way that it will accept one or two leg lengths. If only one leg length is specified, this will be used for both legs. In other words, you should be able to invoke the function with just one leg value:pythagoras(5)
Show answer
pythagoras <- function(a, b = a) {
sqrt(a^2 + b^2)
}7 Quartiles
You can use the summary() function with a numeric vector to obtain descriptive statistics:
set.seed(345)
x1 <- rnorm(100)
summary(x1) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.58206 -0.77032 -0.25206 -0.07596 0.64208 2.81828 7.1 First Quartile
We can write a function to compute the first quartile of a numeric vector; this can be easily done with the advantage of the quantile() function:
quartile1 <- function(x) {
quantile(x, probs = 0.25)
}
# test it
quartile1(mtcars$wt) 25%
2.58125 7.2 Argument na.rm
Many functions that work on vectors have a special argument: na.rm. This parameter is a logical value to indicate whether NAs should be removed or not. Because the quantile() function does come with the na.rm argument, we can take advantage of it and pass it to our quartile1() function:
quartile1 <- function(x, na.rm = FALSE) {
quantile(x, probs = 0.25, na.rm = na.rm)
}Let’s get the weight of cars and add some missing values:
weight_na = mtcars$wt
weight_na[c(1, 10, 20)] = NAIf you apply quartile1() on weight_na using the default call, you will get an error:
quartile1(weight_na)Error in quantile.default(x, probs = 0.25, na.rm = na.rm): missing values and NaN's not allowed if 'na.rm' is FALSETo remove missing values, you can use na.rm = TRUE:
quartile1(weight_na, na.rm = TRUE) 25%
2.77 7.3 Your Turn: function quartiles()
As you know, there are three quartile values—lower quartile, median, and upper quartile—to divide the data set into four ranges, each containing 25% of the data points.
Write a function quartiles() that computes the three quartiles of a numeric vector. BTW: take advantage of quantile()’s vectorization behavior. Include an argument na.rm to decide whether missing values should be removed. Give a default value of na.rm = FALSE.
Test quartiles() on vectors mtcars$wt and weight_na.
Show answer
quartiles <- function(x = 1, na.rm = FALSE) {
quantile(x, probs = c(0.25, 0.50, 0.75), na.rm = na.rm)
}8 Gaussian Function
The Gaussian (Normal) function, given in the equation below, is one of the most widely used functions in science and statistics:
\[ f(x) = \frac{1}{\sqrt{2 \pi} \ s} \ \exp \left \{ -\frac{1}{2} \left (\frac{x - m}{s} \right)^2 \right \} \]
The parameters \(s\) and \(m\) are real numbers, where \(s\) must be greater than zero.
Make a function gaussian() that takes three arguments: x, m, and s. Evaluate the function with \(m = 0\), \(s = 2\), and \(x = 1\).
Show answer
gaussian <- function(x = 1, m = 0, s = 1) {
constant = 1 / sqrt(2*pi)
exponent = exp(-0.5 * ((x-m)/s)^2)
output = constant * (1/s) * exponent
output
}Test your gaussian() function and compare it with the R function dnorm()
# compare with dnorm()
dnorm(x = 1, mean = 0, sd = 2)[1] 0.1760327Now try gaussian() with a vector seq(-4.5, 4.5, by = 0.1), and pass the values to plot() to get a normal curve. Here’s code with values obtained from dnorm()
# you should get a plot like this one
x_values <- seq(-4.5, 4.5, by = 0.1)
y_values <- dnorm(x_values, mean = 0, sd = 2)
plot(x_values, y_values, las = 1, type = "l", lwd = 2)
Your turn:
Show answer
x_values <- seq(-4.5, 4.5, by = 0.1)
y_values <- gaussian(x_values, m = 0, s = 2)
plot(x_values, y_values, las = 1, type = "l", lwd = 2)