- Multiple functions: Doing different things to the same object
- Sub-functions: Breaking up big jobs into small ones
- Example: Back to resource allocation
Reading for Friday: 1.3, 7.3–7.5, 7.11, 7.13 of Matloff (skipping "extended examples")
17
Reading for Friday: 1.3, 7.3–7.5, 7.11, 7.13 of Matloff (skipping "extended examples")
Functions tie together related commands:
my.clever.function <- function(an.argument,another.argument) {
# many lines of clever calculations
return(important.result)
}
Inputs/arguments and outputs/return values define the interface
A user only cares about turning inputs into outputs correctly
Meta-problems:
Meta-solutions:
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
(to the extent possible)
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
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
}
# 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)
}
When one function calls another, use 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)
}
Solve big problems by dividing them into a few sub-problems
Rule of thumb: A function longer than a page is probably too long
Defining a function inside another function
Rule of thumb: If you find yourself writing the same code in multiple places, make it a separate function
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)
}
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))
}
}
or multiple calls:
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?
planner <- function(output,factory,available,slack,tweak=0.1) {
needed <- plan.needs(output,factory)
if (all(needed <= available) && all(available-needed <= slack)) {
return(list(output=output,needed=needed))
}
else {
output <- adjust.plan(output,needed,available,tweak)
return(planner(output,factory,available,slack))
}
}
plan.needs <- function(output,factory) { factory %*% output }
adjust.plan <- function(output,needed,available,tweak) {
if (all(needed >= available)) { return(output*(1-tweak)) }
if (all(needed < available)) { return((1+tweak)) }
return(output*runif(n=length(output),min=1-tweak,max=1+tweak))
}
Next time: Designing functions from the top down