---
title: 'Writing Functions'
author: "James M. Flegal"
output:
ioslides_presentation:
smaller: yes
beamer_presentation: default
---
## Agenda
- Defining functions: Tying related commands into bundles
- Interfaces: Controlling what the function can see and do
- Example: Parameter estimation code
- Multiple functions
- [Recursion](https://en.wikipedia.org/wiki/Recursion_(computer_science)): Making hard problems simpler
## Why Functions?
- Data structures tie related values into one object
- Functions tie related commands into one object
- In both cases: easier to understand, easier to work with, easier to build into larger things
## Example cubic function
```{r}
cube <- function(x) x ^ 3
cube
cube(3)
cube(1:10)
```
##
```{r}
cube(matrix(1:8, 2, 4))
matrix(cube(1:8), 2, 4)
# cube(array(1:24, c(2, 3, 4))) # cube each element in an array
mode(cube)
```
## Example
```{r}
# "Robust" loss function, for outlier-resistant regression
# Inputs: vector of numbers (x)
# Outputs: vector with x^2 for small entries, 2|x|-1 for large ones
psi.1 <- function(x) {
psi <- ifelse(x^2 > 1, 2*abs(x)-1, x^2)
return(psi)
}
```
Our functions get used just like the built-in ones:
```{r}
z <- c(-0.5,-5,0.9,9)
psi.1(z)
```
##
Go back to the declaration and look at the parts:
```{r}
# "Robust" loss function, for outlier-resistant regression
# Inputs: vector of numbers (x)
# Outputs: vector with x^2 for small entries, |x| for large ones
psi.1 <- function(x) {
psi <- ifelse(x^2 > 1, 2*abs(x)-1, x^2)
return(psi)
}
```
**Interfaces**: the **inputs** or **arguments**; the **outputs** or **return value**
Calls other functions `ifelse()`, `abs()`, operators `^` and `>`, and could also call other functions we've written
`return()` says what the output is; alternately, return the last evaluation
**Comments**: Not required by R, but a good idea
## What should be a function?
- Things you're going to re-run, especially if it will be re-run with changes
- Chunks of code you keep highlighting and hitting return on
- Chunks of code which are small parts of bigger analyses
- Chunks which are very similar to other chunks
## Named and default arguments
```{r}
# "Robust" loss function, for outlier-resistant regression
# Inputs: vector of numbers (x), scale for crossover (c)
# Outputs: vector with x^2 for small entries, 2c|x|-c^2 for large ones
psi.2 <- function(x,c=1) {
psi <- ifelse(x^2 > c^2, 2*c*abs(x)-c^2, x^2)
return(psi)
}
```
```{r}
identical(psi.1(z), psi.2(z,c=1))
```
##
Default values get used if names are missing:
```{r}
identical(psi.2(z,c=1), psi.2(z))
```
Named arguments can go in any order when explicitly tagged:
```{r}
identical(psi.2(x=z,c=2), psi.2(c=2,x=z))
```
## Checking Arguments
_Problem_: Odd behavior when arguments aren't as we expect
```{r}
psi.2(x=z,c=c(1,1,1,10))
psi.2(x=z,c=-1)
```
##
_Solution_: Put little sanity checks into the code
```{r}
# "Robust" loss function, for outlier-resistant regression
# Inputs: vector of numbers (x), scale for crossover (c)
# Outputs: vector with x^2 for small entries, 2c|x|-c^2 for large ones
psi.3 <- function(x,c=1) {
# Scale should be a single positive number
stopifnot(length(c) == 1,c>0)
psi <- ifelse(x^2 > c^2, 2*c*abs(x)-c^2, x^2)
return(psi)
}
```
Arguments to `stopifnot()` are a series of expressions which should all
be TRUE; execution halts, with error message, at _first_ FALSE
(try it!)
## What the function can see and do
- Each function has its own environment
- Names here over-ride names in the global environment
- Internal environment starts with the named arguments
- Assignments inside the function only change the internal environment (There _are_ ways around this, but they are difficult and best avoided)
- Names undefined in the function are looked for in the environment the function gets called from _not_ the environment of definition
## Internal environment examples
```{r}
x <- 7
y <- c("A","C","G","T","U")
adder <- function(y) { x<- x+y; return(x) }
adder(1)
x
y
```
##
```{r}
circle.area <- function(r) { return(pi*r^2) }
circle.area(c(1,2,3))
truepi <- pi
pi <- 3
circle.area(c(1,2,3))
pi <- truepi # Restore sanity
circle.area(c(1,2,3))
```
## Respect the interfaces
- Interfaces mark out a controlled inner environment for our code
- Interact with the rest of the system only at the interface
- Advice: arguments explicitly give the function all the information
+ Reduces risk of confusion and error
+ Exception: true universals like $\pi$
- Likewise, output should only be through the return value
## Fitting a Model
Fact: bigger cities tend to produce more economically per capita
A proposed statistical model (Geoffrey West et al.):
\[
Y = y_0 N^{a} + \mathrm{noise}
\]
where $Y$ is the per-capita "gross metropolitan product" of a city, $N$ is
its population, and $y_0$ and $a$ are parameters
##
```{r}
gmp <- read.table("gmp.dat")
gmp$pop <- gmp$gmp/gmp$pcgmp
plot(pcgmp~pop, data=gmp, log="x", xlab="Population", ylab="Per-Capita Economic
Output ($/person-year)", main="US Metropolitan Areas, 2006")
curve(6611*x^(1/8),add=TRUE,col="blue")
```
## Fitting a Model
\[
Y = y_0 N^{a} + \mathrm{noise}
\]
Take $y_0=6611$ for now and estimate $a$ by minimizing the mean squared error
Approximate the derivative of error w.r.t $a$ and move against it
\[
\begin{eqnarray*}
MSE(a) & \equiv & \frac{1}{n}\sum_{i=1}^{n}{(Y_i - y_0 N_i^a)^2}\\
MSE^{\prime}(a) & \approx & \frac{MSE(a+h) - MSE(a)}{h}\\
a_{t+1} - a_{t} & \propto & -MSE^{\prime}(a)
\end{eqnarray*}
\]
##
An actual first attempt at code:
```{r}
maximum.iterations <- 100
deriv.step <- 1/1000
step.scale <- 1e-12
stopping.deriv <- 1/100
iteration <- 0
deriv <- Inf
a <- 0.15
while ((iteration < maximum.iterations) && (deriv > stopping.deriv)) {
iteration <- iteration + 1
mse.1 <- mean((gmp$pcgmp - 6611*gmp$pop^a)^2)
mse.2 <- mean((gmp$pcgmp - 6611*gmp$pop^(a+deriv.step))^2)
deriv <- (mse.2 - mse.1)/deriv.step
a <- a - step.scale*deriv
}
list(a=a,iterations=iteration,converged=(iteration < maximum.iterations))
```
## What's wrong with this?
- Not _encapsulated_: Re-run by cutting and pasting code --- but how
much of it? Also, hard to make part of something larger
- _Inflexible_: To change initial guess at $a$, have to edit, cut,
paste, and re-run
- _Error-prone_: To change the data set, have to edit, cut, paste,
re-run, and hope that all the edits are consistent
- _Hard to fix_: should stop when _absolute value_ of derivative is
small, but this stops when large and negative. Imagine having five copies of
this and needing to fix same bug on each.
Will turn this into a function and then improve it
##
First attempt, with logic fix:
```{r}
estimate.scaling.exponent.1 <- function(a) {
maximum.iterations <- 100
deriv.step <- 1/1000
step.scale <- 1e-12
stopping.deriv <- 1/100
iteration <- 0
deriv <- Inf
while ((iteration < maximum.iterations) && (abs(deriv) > stopping.deriv)) {
iteration <- iteration + 1
mse.1 <- mean((gmp$pcgmp - 6611*gmp$pop^a)^2)
mse.2 <- mean((gmp$pcgmp - 6611*gmp$pop^(a+deriv.step))^2)
deriv <- (mse.2 - mse.1)/deriv.step
a <- a - step.scale*deriv
}
fit <- list(a=a,iterations=iteration,
converged=(iteration < maximum.iterations))
return(fit)
}
```
##
_Problem_: All those magic numbers!
_Solution_: Make them defaults
```{r}
estimate.scaling.exponent.2 <- function(a, y0=6611,
maximum.iterations=100, deriv.step = .001,
step.scale = 1e-12, stopping.deriv = .01) {
iteration <- 0
deriv <- Inf
while ((iteration < maximum.iterations) && (abs(deriv) > stopping.deriv)) {
iteration <- iteration + 1
mse.1 <- mean((gmp$pcgmp - y0*gmp$pop^a)^2)
mse.2 <- mean((gmp$pcgmp - y0*gmp$pop^(a+deriv.step))^2)
deriv <- (mse.2 - mse.1)/deriv.step
a <- a - step.scale*deriv
}
fit <- list(a=a,iterations=iteration,
converged=(iteration < maximum.iterations))
return(fit)
}
```
##
_Problem:_ Why type out the same calculation of the MSE twice?
_Solution:_ Declare a function
```{r}
estimate.scaling.exponent.3 <- function(a, y0=6611,
maximum.iterations=100, deriv.step = .001,
step.scale = 1e-12, stopping.deriv = .01) {
iteration <- 0
deriv <- Inf
mse <- function(a) { mean((gmp$pcgmp - y0*gmp$pop^a)^2) }
while ((iteration < maximum.iterations) && (abs(deriv) > stopping.deriv)) {
iteration <- iteration + 1
deriv <- (mse(a+deriv.step) - mse(a))/deriv.step
a <- a - step.scale*deriv
}
fit <- list(a=a,iterations=iteration,
converged=(iteration < maximum.iterations))
return(fit)
}
```
`mse()` declared inside the function, so it can see `y0`, but it's not added to the global environment
##
_Problem:_ Locked in to using specific columns of `gmp`; shouldn't have to re-write just to compare two data sets
_Solution:_ More arguments, with defaults
```{r}
estimate.scaling.exponent.4 <- function(a, y0=6611,
response=gmp$pcgmp, predictor = gmp$pop,
maximum.iterations=100, deriv.step = .001,
step.scale = 1e-12, stopping.deriv = .01) {
iteration <- 0
deriv <- Inf
mse <- function(a) { mean((response - y0*predictor^a)^2) }
while ((iteration < maximum.iterations) && (abs(deriv) > stopping.deriv)) {
iteration <- iteration + 1
deriv <- (mse(a+deriv.step) - mse(a))/deriv.step
a <- a - step.scale*deriv
}
fit <- list(a=a,iterations=iteration,
converged=(iteration < maximum.iterations))
return(fit)
}
```
##
Respecting the interfaces: We could turn the `while()` loop into
a `for()` loop, and nothing outside the function would care
```{r}
estimate.scaling.exponent.5 <- function(a, y0=6611,
response=gmp$pcgmp, predictor = gmp$pop,
maximum.iterations=100, deriv.step = .001,
step.scale = 1e-12, stopping.deriv = .01) {
mse <- function(a) { mean((response - y0*predictor^a)^2) }
for (iteration in 1:maximum.iterations) {
deriv <- (mse(a+deriv.step) - mse(a))/deriv.step
a <- a - step.scale*deriv
if (abs(deriv) <= stopping.deriv) { break() }
}
fit <- list(a=a,iterations=iteration,
converged=(iteration < maximum.iterations))
return(fit)
}
```
## What have we done?
The final code is shorter, clearer, more flexible, and more re-usable
_Exercise:_ Run the code with the default values to get
an estimate of $a$; plot the curve along with the data points
_Exercise:_ Randomly remove one data point --- how much does the estimate change?
_Exercise:_ Run the code from multiple starting points --- how different
are the estimates of $a$?
## How We Extend Functions
- Multiple functions: Doing different things to the same object
- Sub-functions: Breaking up big jobs into small ones
## Why Multiple Functions?
Meta-problems:
- You've got more than one problem
- Your problem is too hard to solve in one step
- You keep solving the same problems
Meta-solutions:
- Write multiple functions, which rely on each other
- Split your problem, and write functions for the pieces
- Solve the recurring problems once, and re-use the solutions
## Writing Multiple Related Functions
Statisticians want to do lots of things with their models: estimate, predict,
visualize, test, compare, simulate, uncertainty, ...
Write multiple functions to do these things
Make the model one object; assume it has certain components
## Consistent Interfaces
- Functions for the same kind of object should use the same arguments, and presume the same structure
- Functions for the same kind of task should use the same arguments, and
return the same sort of value (to the extent possible)
## Keep related things together
- Put all the related functions in a single file
- Source them together
- Use comments to note ***dependencies***
## Power-Law Scaling
Remember the model:
$$Y = y_0 N^{a} + \mathrm{noise}$$
$$(\mathrm{output}\ \mathrm{per}\ \mathrm{person}) = $$
$$ (\mathrm{baseline}) (\mathrm{population})^{\mathrm{scaling}\ \mathrm{exponent}} + \mathrm{noise}$$
Estimated parameters $a$, $y_0$ by minimizing the mean squared error
Exercise: Modify the estimation code from last time so it returns a list, with components `a` and `y0`
## Predicting from a Fitted Model
Predict values from the power-law model:
```
# Predict response values from a power-law scaling model
# Inputs: fitted power-law model (object), vector of values at which to make
# predictions at (newdata)
# Outputs: vector of predicted response values
predict.plm <- function(object, newdata) {
# Check that object has the right components
stopifnot("a" %in% names(object), "y0" %in% names(object))
a <- object$a
y0 <- object$y0
# Sanity check the inputs
stopifnot(is.numeric(a),length(a)==1)
stopifnot(is.numeric(y0),length(y0)==1)
stopifnot(is.numeric(newdata))
return(y0*newdata^a) # Actual calculation and return
}
```
## Predicting from a Fitted Model
```
# Plot fitted curve from power law model over specified range
# Inputs: list containing parameters (plm), start and end of range (from, to)
# Outputs: TRUE, silently, if successful
# Side-effect: Makes the plot
plot.plm.1 <- function(plm,from,to) {
# Take sanity-checking of parameters as read
y0 <- plm$y0 # Extract parameters
a <- plm$a
f <- function(x) { return(y0*x^a) }
curve(f(x),from=from,to=to)
# Return with no visible value on the terminal
invisible(TRUE)
}
```
## Predicting from a Fitted Model
When one function calls another, use \texttt{...} as a meta-argument, to pass along unspecified inputs to the called function:
```
plot.plm.2 <- function(plm,...) {
y0 <- plm$y0
a <- plm$a
f <- function(x) { return(y0*x^a) }
# from and to are possible arguments to curve()
curve(f(x), ...)
invisible(TRUE)
}
```
## Sub-Functions
Solve big problems by dividing them into a few sub-problems
- Easier to understand, get the big picture at a glance
- Easier to fix, improve and modify
- Easier to design
- Easier to re-use solutions to recurring sub-problems
Rule of thumb: A function longer than a page is probably too long
## Sub-Functions or Separate Functions?
Defining a function inside another function
- Pros: Simpler code, access to local variables, doesn't clutter workspace
- Cons: Gets re-declared each time, can't access in global environment (or in other functions)
- Alternative: Declare the function in the same file, source them together
Rule of thumb: If you find yourself writing the same code in multiple places, make it a separate function
## Plotting a Power-Law Model
Our old plotting function calculated the fitted values
But so does our prediction function
```
plot.plm.3 <- function(plm,from,to,n=101,...) {
x <- seq(from=from,to=to,length.out=n)
y <- predict.plm(object=plm,newdata=x)
plot(x,y,...)
invisible(TRUE)
}
```
## Recursion
Reduce the problem to an easier one of the same form:
```
my.factorial <- function(n) {
if (n == 1) {
return(1)
} else {
return(n*my.factorial(n-1))
}
}
```
## Recursion
Or multiple calls ([Fibonacci numbers](https://en.wikipedia.org/wiki/Fibonacci_number)):
```
fib <- function(n) {
if ( (n==1) || (n==0) ) {
return(1)
} else {
return (fib(n-1) + fib(n-2))
}
}
```
Exercise: Convince yourself that any loop can be replaced by
recursion; can you always replace recursion with a loop?
## Summary
- **Functions** bundle related commands together into objects: easier to
re-run, easier to re-use, easier to combine, easier to modify, less risk of error, easier to think about
- **Interfaces** control what the function can see (arguments,
environment) and change (its internals, its return value)
- **Calling** functions we define works just like calling built-in
functions: named arguments, defaults
- **Multiple functions** let us do multiple related jobs, either on the same object or on similar ones
- **Sub-functions** let us break big problems into smaller ones, and re-use the solutions to the smaller ones
- [Recursion](https://en.wikipedia.org/wiki/Recursion_(computer_science)) is a powerful way of making hard problems simpler