######################################################################## ## Sun Sep 18 17:05:16 PDT 2011 ## ~carlm/213/ProgrammingI/demonstration.r ## ## Demonstrates grouping commands; using if/else and loops ## and functions. These are the building blocks of a complete ##programming language. ######################################################################## ## 1) everything you write in R is an "expression" -- with an ## optional assignment to an object ## expressions are mathematical statements composed of functions and ## objects. Expressions always produce something which can be ## assigned to an object. ##2) expressions can be grouped together in { } in which ##case the last expression value is returned for optional assignment ############################################# ## Example 0 ## A simple example of a "grouped expression" where in the result of ##the last expression in the group is returend. ############################################# median.uniform <- { x<- runif(1000) y<-sort(x) (y[500]+y[501])/2 } median.uniform ## Although only the last result of last expression in a "group" ## is returned such that it can be assigned, the objects created ## inside the {} become part of the workspace and *are* accessible head(x) head(y) ## grouping expressions is not all that useful alone. The above ## grouped expresion is no different from running this: x<- runif(1000) y<-sort(x) median.uniform <- (y[500]+y[501])/2 median.uniform ## And the ungrouped version uses TWO fewer keystrokes. ## Yet grouping expressions is the basis of "loops", "branches" and ## "functions" all of which are tremendously useful. ## Loops [for() and while()] allow for repeated execution of a grouped ## expression over a set of different starting values. ## Branches [if()/else] allow for conditional execution -- doing one ## thing under one condition and another under another. ## functions are grouped expressions that can be run on a set of ## inputs in a way that insulates them from the workspace. ############################################# ## QUESTION: ### write very simple for() loop that take the square of each ### element of a vector of any length. The result of the loop should ### be a vector, of the same length as the original vector. For this ### purpose you may as well use the vector "xvec" below and call your ### result "yvec" ############################################# xvec <- -(x<-runif(1)*1000):x ## HINT length(xvec) ## Your answer goes here: ## plot the result: plot(yvec~xvec) ## or to embelish a little: plot(yvec~xvec, col=rainbow(15),cex=yvec/1000) ## Q: how does your loop based solution differ from y<-x^2 ? ## A: it requires more typing ## Suppose that we wanted to take the log() of each element of xvec? ## Taking the log of a negative number is not considered mathematically cool. ## One solution is to use if/else to "branch" inside of a loop ############################################# ## QUESTION 2: Add code where it says "your code goes here" such that ## the resulting loop will take the log of positive numbers but will ## assign 0 to negative numbers xvec <- -(x<-runif(1)*1000):x for(i in 1:length(xvec)){ ##Your code goes here } ## Since R is so clever about mixing scalars and ## vectors one might ask... ## Q: Why bother with a loop? why not just do this? logxvec<-if(xvec>0){log(xvec)} else {0} ## A: It won't work ## But THIS might. logxvec<-ifelse(xvec>0,log(xvec),0) ## Read on to discover what's up with the warning message. ######################################################################## ## if(){}else{} vs ifelse ######################################################################## ## if() TAKES A SCALAR VALUED EXPRESSION as an argument -- it's only ## purpose is to decide whether or not to run the single command OR ## {grouped code} that follows it. It is best thought of as a ## programming structure only. ## ifelse() OPERATES LIKE A TYPICAL FUNCTION IN R. If you give it ## vectors it returns a vector if you give it something else it ## returns something else. ## But as we saw above, ifelse doesn't really branch. It executes the ## BOTH expressions and chooses which to assign based on the the first ## arguement. Thus the warning messages. ## the scalar version, if() is useful in situations where decisions ## are called for. the ifelse() function is useful by itself. ############################################# ## Not to beat a dead horse, but the if()/else / ifelse() difference ## is important and subtle so here is another example: ## A for loop with an if/else branch (just like the one you wrote in ## Question 2. Followed by an ifelse() expression that doesn't do the ## same thing ############################################# ## evaluate a test expression (sum(x) > 5 ) and if it's true then ## do something -- perhps do a lot of things; if not do something else x<-runif(10);y<-1:10 if(sum(x) > 5){ ### run this code if x sums to more than 5 y<-y*100 x<-sqrt(x) plot(x,y,main="x > 5... Behold",col='red') } else { ### otherwise run this grouped code x<-x*50 y<-sqrt(x) plot(x,y,main="x <= 5... Behold this instead",col='blue') } ############################################# ## ifelse() is NOT useful in the above situation. In the above case we ## need to make just ONE decision: is sum(x) GT or LT then 5. And ## based on that decision we fiddle all the values of x and draw a ## plot. ifelse is useful when we need to inspect *each* value of x and do one ## of two things depending on that ONE element of x. ############################################# x<-runif(10);y<-1:10 plot(y~x,col=ifelse(x > .5,'red','blue')) ############################################# ## while() loops: ## we have used for() loops to execute code a predetermined number of ## times. We even cleverly extended it by using for(x in ## 1:length(xvec)) BUT even here when the for() loop began, it's exit ## was predetermined -- even if not known to us directly. while() ## loops are for when no one knows. ############################################# ## Here is a for() loop that illustrates a hypothetical frog jumping ## half way to the end of a hypothetical log exactly 100 times ## Set up an empty plot so that we can watch the frog plot(c(0,100),c(0,100), type='n',xlab='jump', ylab='position') start<-1 end<-100 points(1,start,cex=5,col='red') # the original point is in red #NOTE how points() adds a #point to an already #established set of axes for (i in 1:100){ ## plot exactly 100 points and then stop start<- start + .5*(end-start) # the new start is half way to # the end points(i, start, cex=max(5/i,.2), col='blue') } abline(h=100) # a horizontal line at 100 ### while() is used to execute as long as a condition is met ############################################# ## Now the same idea BUT this time we want to see just how many times ## the frog can actually jump half way to the end of a log ... on the ## planet R. plot(c(0,100),c(0,100), type='n',xlab='jump', ylab='position') start<-1 end<-100 points(1,start,cex=5,col='purple') i<-1 while(start < end){ ## keep on plotting points until start < end i<-i+1; start<- start + .5*(end-start) points(i,start, cex=.2, col='blue') print(c(i,format(start,digits=17,scientific=FALSE))) } abline(h=100) ############################################# ## Functions -- stored expressions ############################################# ## A "function" is a stored group of expressions/commands with an ## "entrance" for values to be passed in, and an exit by which results ## are "returned" #### ### the function "function()" returns objects of type "function". Is ### that an informative sentance or what? ## this code creates a function called choose which expects a ## single argument, 'stuff'. After executing the code below, you ## will have new function in your workspace. choose<- function(stuff){ ## returns a randomly selected element of the vector 'stuff' sel<- runif(min=1,max=length(stuff)+1,n=1) sel<-trunc(sel) return(stuff[sel]) } ## see the object which you have just created: choose ## executes the choose() function choose(1:125) choose(1:25) choose(1:125) choose(c('red','yellow','blue','purple','papayawhip')) choose(c('red','yellow','blue','purple','papayawhip')) ## # # # # # # # # # # # # # # # # # # # # # # # # # ### Critical Idea: variables created inside a function will be ### destroyed when the function completes execution. Prexisting ### variables which are changed inside a function will revert to their ### previous values when the function completes execution ## # # # # # # # # # # # # # # # # # # # # # # # # # ## create the pointless function which simply adds arg to previous ## where previous must already is the workspace (is it a good idea to ## write function that rely on objects that exist in the workspace?) pointless<- function(arg){ ## illustrates "scope" inside functions previous<- previous+ arg return(previous) } previous<-13 ## previous is set to 13 BEFORE function execution pointless(5) ## pointless adds 5 to previous and returns it previous ## huh? didn't we add 5 to previous? ## A slightly more useful function footer<-function(){ mtext(side=1,line=3,cex=.8,adj=1,text=paste("Your Name",date(),sep="\n")) } ######################################################################## ## The Project: A stochastic simulation ## Consider a simplified but stochastic Malthusian world: wages are ## determined by the population at the start of each period, and ## population growth is determined by the period's wage rate + a ## stochastic (random) element. Below is a set of five STEPS with ## largely fill-in-the-blank sort of coding. Do what each says and ## then you can go home. At the bottom of this file are answers, for ## those who prefer less challenge. ######################################################################## ## STEP 1: write a function that accepts a population and returns a wage rate. wageRate<-function(pop,maxpop=10000){ if(missing(pop)) stop ("must supply pop(ulation)") ## maxpop defaults to 100K -- if it is not set when the funtion is ##executed that is the value it will have. ## YOUR CODE GOES HERE ## As population approaches maxpop, the growth rate should decline ## but should stay positive. if it exceeds maxpop, wages should fall ## abruptly and become negative. Keep the wage between -1 and 1. ######END your code return(wage) } ## test your function wageRate(100) wageRate(100000) ### ## Note the use of the default value for maxpop ### ## STEP 2: write a function that accepts wage ## and returns a stochastically determined population growth rate popGrate<-function(wage){ ##---- ## Requires a wage rate returns a stochastically ## determined growth rate which is corelated with the wage and ## agrees in sign ### ## Note the use of if(missing()), and stop() to check for bad arguments. ## if (missing(wage)) stop ("must supply wage") ## YOUR CODE GOES HERE: ## generate a random growth maybe use runif(1) or rnorm(1); ## make it a function of the wage -- higher wage --> higher growth rate. ## if the wage is negative, the growth rate must be negative ## END your code return(grate) } ## let's try it ... 10 or so times currentPop<-1000 for (w in seq( .9,-.9, length=25)){ currentPop<-currentPop*(1+popGrate(wage=w)) print(c(w,currentPop)) } ## Is the population dropping fast enough when wages are very low? ## whatever, it's just a simulation. ### STEP 3 ## Now write a program that consists of a while() loop that executes ## the simulation year by year and prints the result to the screen ## cpop<-5000 year<-0 while (year < 100){ ### YOUR CODE GOES HERE ################## print(cpop) } ## Does it work? ... great. ### # Note the use of the while loop requires that the index variable be # updated with each iteration. The for() loop takes care of that for you. ### # Adding graphics ### ## STEP 4 Add graphics ## Copy the bits that you wrote for Step 3 into the code below--which # does only one new thing cpop<-5000 year<-0 pophist<-NULL while (year < 100){ ## INSERT CODE FROM STEP 3 ############ pophist<-c(pophist,cpop) } #x11() plot(pophist,type='s',col="#FFaa11") footer() ### # To demonstrate the stochasticity of this thing, let's put the whole # thing into a single function that can be called repeatedly ### do.everything<- function(){ ## INSERT ALL CODE NEEDED TO RUN an entire simulation INSTEAD of the ## plot(pophist, ... ) above, use the lines() command that is ## provided as the last line of the function. Note that this ## function need not return anything -- we're going to use it in a ## loop wherein it will run a simulation and draw a plot again and ## again. cpop<-5000 ## YOUR CODE GOES HERE ## END YOUR CODE lines(pophist,type='s',col=colors[trunc(runif(1)*25)]) } ##### ## STEP 5: run a bunch of trial, sit back and relax. ## create a blank plot plot(c(1,100),c(5000,10000),type='n',xlab='Population', ylab='Time') footer() ## do everything 25 times for(i in 1:25){ do.everything() } ## You have finished the assignment -- if you still have energy and a ## burning desire to learn, have a go at improving the plot that is ## produced in Step 5. 1# ######################################################################### ## FRUSTRATION Prevention ######################################################################## ## ## STEP 1: write a function that accepts a population and returns a wage rate. ## wageRate<-function(pop,maxpop=10000){ ## if(missing(pop)) stop ("must supply pop(ulation)") ## ## maxpop defaults to 100K -- if it is not set when the funtion is ## ##executed that is the value it will have. ## ## YOUR CODE GOES HERE ## ## As population approaches maxpop, the growth rate should decline ## ## but should stay positive. if it exceeds maxpop, wages should fall ## ## abruptly and become negative. Keep the wage between -1 and 1. ## pct.of.max<- pop/maxpop; ## wage<- ifelse(pct.of.max > .9, -.1,pct.of.max) ## return(wage) ## } ## ## test your function ## wageRate(100) ## wageRate(100000) ## ### ## ## Note the use of the default value for maxpop ## ### ## ## STEP 2: write a function that accepts wage ## ## and returns a stochastically determined population growth rate ## popGrate<-function(wage){ ## ##---- ## ## Requires a wage rate returns a stochastically ## ## determined growth rate which is corelated with the wage and ## ## agrees in sign ## ### ## ## Note the use of if(missing()), and stop() to check for bad arguments. ## ## ## if (missing(wage)) stop ("must supply wage") ## ## YOUR CODE GOES HERE: ## ## generate a random growth maybe use runif(1) or rnorm(1); ## ## make it a function of the wage -- higher wage --> higher growth rate. ## ## if the wage is negative, the growth rate must be negative ## grate<- runif(1)*.1*wage ## return(grate) ## } ## ## let's try it ... 10 or so times ## currentPop<-1000 ## for (w in seq( .9,-.9, length=25)){ ## currentPop<-currentPop*(1+popGrate(wage=w)) ## print(c(w,currentPop)) ## } ## ## Is the population dropping fast enough when wages are very low? ## ## whatever, it's just a simulation. ## ### STEP 3 ## ## Now write a program that consists of a while() loop that executes ## ## the simulation year by year and prints the result to the screen ## ## ## cpop<-5000 ## year<-0 ## while (year < 100){ ## ### YOUR CODE GOES HEE ## cwage<-wageRate(cpop) ## cpop<-cpop *(1+ popGrate(cwage)) ## year<-year+1 ## ################## ## print(cpop) ## } ## ## Does it work? ... great. ## ### ## # Note the use of the while loop requires that the index variable be ## # updated with each iteration. The for() loop takes care of that for you. ## ### ## # Adding graphics ## ### ## ## STEP 4 Add graphics ## ## Copy the bits that you wrote for Step 3 into the code below--which ## # does only one new thing ## cpop<-5000 ## year<-0 ## pophist<-NULL ## while (year < 100){ ## ## INSER CODE FROM STEP 3 ## cwage<-wageRate(cpop) ## cpop<-cpop *(1+ popGrate(cwage)) ## year<-year+1 ## ############ ## pophist<-c(pophist,cpop) ## } ## #x11() ## plot(pophist,type='s',col="#FFaa11") ## ### ## # To demonstrate the stochasticity of this thing, let's put the whole ## # thing into a single function that can be called repeatedly ## ### ## do.everything<- function(){ ## ## INSERT ALL CODE NEEDED TO RUN an entire simulation INSTEAD of the ## ## plot(pophist, ... ) above, use the lines() command that is ## ## provided as the last line of the function. Note that this ## ## function need not return anything -- we're going to use it in a ## ## loop wherein it will run a simulation and draw a plot again and ## ## again. ## cpop<-5000 ## year<-0 ## pophist<-NULL ## colors<-rainbow(25) ## while (year < 100){ ## cwage<-wageRate(cpop) ## cpop<-cpop *(1+ popGrate(cwage)) ## year<-year+1 ## pophist<-c(pophist,cpop) ## } ## lines(pophist,type='s',col=colors[trunc(runif(1)*25)]) ## } ## ##### ## ## STEP 5: run a bunch of trial, sit back and relax. ## ## create a blank plot ## plot(c(1,100),c(5000,10000),type='n',xlab='Population', ylab='Time') ## ## do everything 25 times ## for(i in 1:25){ ## do.everything() ## }