Topic #10 – Self-Made Functions

1. Topic #10 – Self-Made Functions • Self-Made Functions • Extremely useful tool: • Can be put in separate script and sourced in different code • Prevents duplication of code • Allows easy modification later • Saves lots of time! • Necessary for nonlinear optimization • Requires a function to be minimized

2. Topic #10 – Self-Made Functions Self-Made Functions Simple syntax Function=function(parameters){ ...expressions... object to return }

3. Topic #10 – Self-Made Functions Self-Made Functions Example #1 Creating a function for use in nonlinear optimization

4. Topic #10 – Self-Made Functions Self-Made Functions Example #1 Y ~ X1^b + X2^b Y=c(10, 15, 14, 12, 9, 10, 12) X1=c(1, 1, 1, 2, 2, 2, 2) X2=c(5, 4, 1, 2, 3, 4, 2) RSS=function(a, b){ Y.Est=X1^a + X2^b RSS=sum( (Y-Y.Est)^2 ) RSS }

5. Topic #10 – Self-Made Functions Self-Made Functions Example #2 Modifying a call-function for a specific use * Typing a call-function without parentheses displays its code* help.search(“gibbs”) library(“LearnBayes”) gibbs

6. Topic #10 – Self-Made Functions Gibbs function > gibbs function (logpost, start, m, scale, data) { p = length(start)‏ vth = array(0, dim = c(m, p))‏ f0 = logpost(start, data)‏ arate = array(0, dim = c(1, p))‏ th0 = start th1 = th0 for (i in 1:m) { for (j in 1:p) { th1[j] = th0[j] + rnorm(1) * scale[j] f1 = logpost(th1, data)‏ u = runif(1) < exp(f1 - f0)‏ th0[j] = th1[j] * (u == 1) + th0[j] * (u == 0)‏ f0 = f1 * (u == 1) + f0 * (u == 0)‏ vth[i, j] = th0[j] arate[j] = arate[j] + u } } arate = arate/m stuff = list(par = vth, accept = arate)‏ return(stuff)‏ }

7. Topic #10 – Self-Made Functions Modified Gibbs function > gibbs.modified=function (logpost, start, m, scale) { p = length(start)‏ vth = array(0, dim = c(m, p))‏ f0 = logpost(start)‏ arate = array(0, dim = c(1, p))‏ th0 = start th1 = th0 for (i in 1:m) { for (j in 1:p) { th1[j] = th0[j] + rnorm(1) * scale[j] f1 = logpost(th1, data)‏ u = runif(1) < exp(f1 - f0)‏ th0[j] = th1[j] * (u == 1) + th0[j] * (u == 0)‏ f0 = f1 * (u == 1) + f0 * (u == 0)‏ vth[i, j] = th0[j] arate[j] = arate[j] + u } } arate = arate/m stuff = list(par = vth, accept = arate)‏ return(stuff)‏ }