# R code accompanying lecture 5 for 36-350, Fall 2012 # See lecture for context # Example of predicting from a model # Predict response values from a power-law scaling model # Presumes: power-law model objects are lists containing a scaling factor # called "y0" and a scaling exponent called "a" # Inputs: fitted power-law model (object), vector of values at which to make # predictions at (newdata) # Argument names follow conventions for other predict() functions # 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 } # Examples of plotting # 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 # 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) } # With other arguments for curve() # illustrates the use of "..." to pass along extra arguments # Plot fitted curve from power law model over specified range # Inputs: list containing parameters (plm), start and end of range (from, to), # other graphical arguments for curve (...) # Outputs: TRUE, silently, if successful # Side-effect: Makes the plot 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) } # Using predict.plm() # Plot fitted curve from power law model over specified range # Inputs: list containing parameters (plm), start and end of range (from, to), # number of points for plotting (n), other graphical arguments for plot (...) # Outputs: TRUE, silently, if successful # Side-effect: Makes the plot 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) } # Examples of recursion # Calculate a factorial # Input: a number (n) # Output: the factorial of n # Presumes: n is a single positive integer my.factorial <- function(n) { if (n == 1) { # Base case return(1) } else { # Recursion return(n*my.factorial(n-1)) } } # Calculate a Fibonacci number # Input: a number (n) # Output: the nth Fibonacci number # Presumes: n is a single positive integer fib <- function(n) { if ( (n==1) || (n==0) ) { # Base cases return(1) } else { # Recursion return (fib(n-1) + fib(n-2)) } } # EXERCISE: Try tracing through fib(5) # Our resource allocation example, made more modular and recursive # Adjust an output plan until it's under budget and close to full utilization # Inputs: vector of initial planned output levels (output), matrix giving # resources needed per unit of output (factory), vector of resources # available (available), vector of acceptable idle resource amounts # (slack), factor by which to adjust plans (tweak) # Output: List containing acceptable output levels (output) and resources # required (needed) 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)) } } # Calculate resources needed by a given output plan # Inputs: vector of output levels (output), matrix giving resources needed per # unit of output (factory) # Output: column matrix giving resources needed plan.needs <- function(output,factory) { factory %*% output } # Make multiplicative adjustments to an unsatisfactory output plan # Inputs: vector of planned output levels (output), resources needed by plan # (needed), vector of resources available (available), factor by which to # adjust plans # Output: vector of new output levels # Presumes: resources needed do not exactly match resources available adjust.plan <- function(output,needed,available,tweak) { # If we use too much, reduce all planned output levels by the same factor if (all(needed >= available)) { return(output*(1-tweak)) } # If we use too little, increase all planned output levels if (all(needed < available)) { return((1+tweak)) } # If we use too much of some things and too little of others, randomly # change all output levels # See help(runif) return(output*runif(n=length(output),min=1-tweak,max=1+tweak)) # N.B., should never call this when needed == available }