How does a distribution of change in species translate to a distribution of changes in function?

changeInFunctionBEF.R

#Michaelis-Menten function with the max function set to 1
#takes h as an argument, the half-saturation
mm <- function(x, h) x/(x+h)

#simple power function
pow <- function(x, p) x^p

#x<-1:100
#plot(x, mm(x), type="l")
#plot(x, pow(x), type="l")



plotBEFDist <- function(init = 30, changeMean = 0, changeSD = init/2, 
                        h=2, p=0.5, 
                        nsim=500, breaks=200, ...){
  
  #get the distribution of species change based on 
  #inputs
  change <- rnorm(nsim, changeMean, changeSD)
  
  #get rid of impossibly low changes
  change <- change[which(!(-1*change>init))]
  
  #calculate the change in function resulting from the distribution 
  #based on a MM curve
  func_change <- (mm(init+change, h)-mm(init, h))/mm(init, h)
  
  #calculate the change in function resulting from the dist.
  #based on a power function
  func_change_power <- (pow(init+change, p)-pow(init, p))/pow(init, p)
  
  #plot the resulting distributions
  par(mfrow=c(1,2))

  hist(func_change, breaks=breaks, 
       main=paste("Michaelis-Menten BEF Relationship"), 
       xlab="Porportion Change in Function", ...)

  hist(func_change_power, breaks=breaks,
       main="Power Function BEF Relationship", 
       xlab="Porportion Change in Function", ...)

  par(mfrow=c(1,1))
  
}

#try it with defaults
plotBEFDist()

#change the half-saturation in the MM curve
#and the coefficient in the power function
plotBEFDist(h=3, p=0.05)