Quantile-Quantile Plot with Unknown Distribution? 
N=2762
I've been exploring a data set that seems to give rise to this kind of plot rather frequently. Would you say this is one population with a different than normal population? Or are two populations confounding the normal distribution?
It used matplotlib and scipy.stats:
      (x,y), (slope, inter, cor) = stats.probplot(data, dist='norm')

      osmf = x.take([0, -1])  # endpoints
      osrf = slope * osmf + inter
      self.ax.plot(x, y, ',', osmf, osrf, '-', **self.kwargs)

 A: There are a variety of different possibilities. For example, a chi-square distribution with degrees of freedom in the range of 30-40 would give rise to such a qq-plot. In R:
x <- rchisq(10000, df=35)
qqnorm(x)
qqline(x)

looks like this:

A mixture of two normals with different means doesn't apply though.
x <- c(rnorm(10000/2, mean=0), rnorm(10000/2, mean=2))
qqnorm(x)
qqline(x)

looks like this:

Note how the points cross the line, which is a different pattern than the one you observe.
A: Your dataset clearly is not normal.  (With this much data, any goodness of fit test will tell you that.)  But you can read much more than that from the normal probability plot:


*

*The generally smooth curvature does not hint at a mixture structure.

*The upper tail is too stretched out (values too high compared to the reference distribution).

*The lower tail is too compressed (values also too high).
This suggests that a mild Box-Cox transformation will produce nearly-normal, or at least symmetric, data.   To find it, consider some key values on this plot: the median, found above the x-value of 0, is about 0.90; +2 standard deviations is about 0.99; and -2 standard deviations is about 0.825.  The nonlinearity is apparent from the simple calculations 0.99 - 0.90 = 0.09 whereas 0.90 - 0.825 = 0.075: the rise from the median to the upper tail is greater than the rise from the lower tail to the median.  We can equalize the slopes by trying out some simple re-expressions of these three values only.  For example, taking the reciprocals of the three key data values (Box-Cox power of -1) gives
1/0.825 = 1.21
1/0.90  = 1.11; 1.21 - 1.11 = 0.10 (new slope is 0.050 per SD)
1/0.99  = 1.01; 1.11 - 1.01 = 0.10 (0.050 per SD)

Because the slopes of the re-expressed values are now equal, we know the plot of reciprocals of the data will be approximately linear between -2 and +2 SDs.  As a check, let's pick more points further out into the tails and see what the reciprocal does to them.  I estimate that the value in the plot at -3 SD from the mean is around 0.79 and the value +3 SD from the mean is 1.05.  The two slopes in question equal 0.053 and 0.052 per SD: close enough to each other and to the slopes found between -2 and +2 SD.
My estimates--based on the plot as shown on a monitor--are crude, so you will want to repeat these (simple, quick) calculations with the actual data.  Nevertheless, there is considerable evidence that your data when suitably re-expressed with a simple transformation will be close to normally distributed.
A: You may want to take a look at the   Anderson-Darling test for normality which empirically tests whether or not your data comes from a given distribution. @chl recommends looking at the scipy toolkit, specifically anderson() in morestats.py for an implementation.
