What are we modelling when a gamma distribution has non-integer shape parameter

Question

I wish to receive a clear and concise answer as to what is being modeled for a gamma distribution with non-integer shape parameter, and a more detailed derivation of its distribution function for all positive numbers (rather than simply integers).

I have derived the gamma distribution as the sum of $n \in \{0, 1, 2, ...\}$ exponentially distributed random variables. I also have shown via contradiction (or rather, my lecturer has), that a situation in which $k-1$ events occur at some given rate $\lambda_1$, and one event occurs at a rate $\lambda_2 = \frac{\lambda_1}{2}$, cannot be modelled by a gamma distribution (no matter the rate parameter, as a compound Poisson process does not follow a gamma distribution, although the two distributions are very close).

Thus, if this the two distributions are not the same, how would I derive the gamma pdf for a non-integer shape parameter?

I have read the answer here: Is there another interpretation for a Gamma distribution with non-integer shape parameter?

However this is with a rate parameter of one, and thus not general enough. All the textbooks I have merely introduce or define the distribution without any explanation or derivation, which I find ridiculous. We may be able to use this definition to prove that it is a valid pmf, however it does not prove that it has the property of modelling the things we wish to model. This is why I am searching for a derivation based on a non-integer shape parameter.

Answer 1

A part of the answer is that there is a limited number of model distributions and lots of phenomena in physics, biology, etc. that we want to model. Gamma distribution can approximate many distributions that behave as a power law for small values of variable and sharply decay in the infinity (see here for some examples).

Mathematically Gamma distribution may also arise as a stationary solution of Fokker-Planck equation: $$ \partial_t P(x,t)=-\partial_xJ(x,t),\\J(x,t)=\mu(x)P(x,t)- \partial_x [D(x)P(x)] $$ with some choices of parameters, like $D(x)=D_0x$ and $\mu(x)=-\mu_1x+\mu_0$, which gives $$ J(x,t)=0\Rightarrow P_0(x)\propto x^{\frac{\mu_0}{D_0}-1}e^{-\frac{\mu_1}{D_0}x} $$ Such equations arise, e.g., in mathematical population genetics.