title: Add fitted (gamma) distribution density curve to histogram plot from library(MASS) fitdistr
body: <p>I would like to plot the gamma density function derived from a set of observations over the histogram of the observed data. I am able to generate the histogram and parameter estimates for the gamma fit. This is being done for multiple subsets of data from a master data set. How can I plot the gamma density function on each of the histograms created in this loop?</p>
  
  <p>I currently have:</p>
  
  <pre><code>library(MASS)
  
  species &lt;- c('acesac', 'acesac', 'acesac', 'acesac', 'acesac', 'acesac',
   'polbif', 'polbif', 'polbif', 'polbif', 'polbif', 'polbif')
  tmean &lt;- c(2,3,5,6,6,7,5,6,6,6,8,9) 
  Data &lt;- data.frame(species, tmean) 
  
  for (i in unique(Data$species)){
    subdata &lt;- subset(Data, species ==i)
    hist(subdata$tmean, main = i)
    dist &lt;- fitdistr(subdata$tmean, "gamma")
  }
  </code></pre>
  
  <p>I'm thinking that I should use <code>lines()</code>, however, not sure how to specify this?</p>