  To extend the Bose-Einstein (BE) distribution to fractional order, we turn
our attention to the differential equation, $df/dx =-f-f^2$. It is satisfied
with the stationary solution, $f(x)=1/(e^{x+\mu}-1)$, of the Kompaneets
equation, where $\mu$ is the constant chemical potential. Setting $R=1/f$, we
obtain a linear differential equation for $R$. Then, the Caputo fractional
derivative of order $p$ ($p>0$) is introduced in place of the derivative of
$x$, and fractional BE distribution is obtained, where function ${\rm e}^x$ is
replaced by the Mittag-Leffler (ML) function $E_p(x^p)$. Using the integral
representation of the ML function, we obtain a new formula. Based on the
analysis of the NASA COBE monopole data, an identity $p\simeq e^{-\mu}$ is
found.
