title: How to plot a Q-Q plot of data versus a custom theoretical function
body: <p>I'd like to make a visual assessment of whether my data conform to a particular distribution function. To do this, I am using R to generate a quantile-quantile (Q-Q) plot. The distribution function is very specific and is not featured in the standard list of probability distributions, so I wrote my own R function to describe it. It is called 'DistFunc' in the code below, and consists of the ratio of two gamma functions.</p>
  
  <p>In brief, what I do in my code is read my data in from a file, 'DistributionEstimate.txt', which contains two columns. Column 1 is the x values and column 2 is the y values. The variables 'a' and 'b' are best-fit parameters which I had determined previously in another program using a least-squares fit of this distribution function to the data. I then define DistFunc and attempt to plot the Q-Q plot with the qqmath function.</p>
  
  <p>The problem occurs at this point. R proceeds to give me a lot of warnings saying that the DistFunc returns values out of range in 'gammafn', and fails to plot anything. That is fair enough, as I know that the function contains a pole close to the origin. As you can see in the code, I attempt to normalise DistFunc to attempt to convert it to a probability distribution (which, I think, is what is required in order to use qqmath?), however, this doesn't help.</p>
  
  <p>Do any of you have any idea how to overcome this problem - for example, by using a different plotting function that doesn't require normalisation, or to convert it to a pseudo-probability distribution, without affecting the outcome too seriously?</p>
  
  <p>I'd be very grateful for any helpful input!</p>
  
  <pre><code>install.packages('lattice')
  library(lattice)
  x&lt;-read.table("C:/DistributionEstimate.txt", colClasses = c(rep("NULL",1),rep("numeric",1)), header = FALSE)
  y&lt;-read.table("C:/DistributionEstimate.txt", colClasses = c(rep("numeric",1),rep("NULL",1)), header = FALSE)
  x&lt;-sapply(x, as.numeric)
  y&lt;-sapply(y, as.numeric)
  a&lt;-16359727025.407821410;
  b&lt;-198838619.13262583836;
  DistFunc &lt;- function(k,ampl=a,stretch=b) {
      fdist&lt;-ampl*gamma(k*stretch-1/2)/gamma(k*stretch+1)
      fnorm&lt;-fdist/sum(fdist)
  }
  qqmath(DistFunc(x), y, col="blue", envelope=.95, xlab="Quantiles of the best-fit model", ylab="Quantiles of the data")
  abline(0,1, col="red", lwd=2)
  grid()
  </code></pre>