1

I've got data that, when plotted on a log-log scale, looks like two straight lines that gently bend and blend into each other. See the picture below.

I don't know how to find a function to fit that. Now, I know that I should use some kind of broken power law, and I know how to make one that asymptotically converges to a constant, but I don't know how to combine two power laws. I've found something that looks just right as far as its shape is concerned, and even some Python code to generate it, but I can't figure out how to tweak it to fit my data without breaking the continuity between the two power laws:

def chabrier03individual(m):
        k = 0.158 * exp(-(-log(0.08))**2/(2 * 0.69**2))
        return numpy.where(m <= 1,
                0.158*(1./m) * exp(-(log(m)-log(0.08))**2/(2 * 0.69**2)), k*m**-2.3)

Note that here, log means $log_{10}$. Am I even on the right track?

Two power laws blending into each other

Alexko
  • 177

3 Answers3

3

I've always liked the very simple

$$y=\frac{1}{\frac{1}{x^a}+\frac{1}{x^{-b}}}$$

where both $a$ and $b$ are positive.

When $x$ is small, the $x^a$-term dominates and you get a slope of $a$ on the log-log plot, whereas when $x$ is large the $x^{-b}$-term dominates with a slope of $-b$ on the log-log plot.

Of course this can be modified with additional constant multipliers in the denominator.

Cye Waldman
  • 7,524
  • That is nicely simple indeed! It might even be simple enough for me, thanks! I'll see how good a fit I can get. – Alexko Sep 12 '17 at 21:18
  • How delightfully simple! Just a quick fit by hand: https://s26.postimg.org/mwy25dtx5/moutain_Fit.png ; I'll see whether SciPy can do significantly better. Thanks again! – Alexko Sep 12 '17 at 21:41
2

After many years this question was posted, I came across this problem and found this thread. However, my problem is more general: how do I smoothly merge two or more power laws?

The goal is to find a function that has a derivative ($\frac{d\log y}{d\log x}$) that looks like a steps function with the discontinuity smoothly joined. The most natural candidate is the logistic function: $$ f(x) = \frac{L}{1 + \exp(-k(x - x_0))}. $$

Given sequences of power-law indices $m_1, m_2, ..., m_n$ and turn-over points (in ascending order) $x_{1,2}, x_{2,3}, ..., x_{n-1,n}$, choosing $k = 1$, we can write the derivative as $$ \frac{d\log y}{d\log x} = m_1 + \frac{m_2 - m_1}{1 + \exp(-\log x + \log x_{1,2})} + ... + \frac{m_n - m_{n-1}}{1 + \exp(-\log x + \log x_{n-1,n})}. $$

It turns out that this differential equation can be rewritten to a very simple form: $$ \frac{d\log y}{dx} = \frac{m_1}{x} + \frac{m_2 - m_1}{x + x_{1,2}} + ... + \frac{m_n - m_{n-1}}{x + x_{n-1,n}}. $$

After integration, it reduces to a very simple form that one should have thought of at the beginning: $$ y = A x^{m_1} (x + x_{1,2})^{m_2 - m_1} ... (x + x_{n-1,n})^{m_n - m_{n-1}} $$

The reason it works is that, for some $x$ that falls between two turn-over points $x_{i-1,i}$ and $x_{i,i+1}$, $x$ is much larger than $x_{j-1,j}$ for all $j \leq i$ while $x$ is much smaller than $x_{j,j+1}$ for all $j \geq i$. Therefore, all former terms reduces to $x^{m_i}$ and all latter terms reduces to a constant. And there you have it!

1

What you provided in the previous question looked like a the cumulative data function (CDF), or survival function (SF), and this the probability density function (PDF). So this is a broken power law with a single break. You can either use the negative derivative of the SF from the previous question, which gives: \begin{align} f(x) & \equiv - N\frac{\partial B(x, s, \beta, x_n)}{\partial x} \\ & = - N\beta \left[1 + \left(\frac{x}{x_n}\right)^{|\beta|s}\right]^{\operatorname{sgn}(\beta)/s - 1} \left(\frac{x^{|\beta|s-1}}{x_n^{|\beta|s}}\right)^{|\beta|s}, \end{align} or, if you'd prefer to fit the pdf with its own smoothly broken power law that has a single break, you'd go for: $$ f(x) = N x^n \left[1+\left(\frac{x}{x_n}\right)^{|\beta|s}\right]^{\operatorname{sgn}(\beta)/s}, $$ which gives one more parameter than fitting the SF with a broken power law.

Sean Lake
  • 1,737
  • That did the trick (the second one), thank you very much! I'll try to optimize both this and Cye's suggestion, and see what ends up fitting best.

    Picture of the fit

    – Alexko Sep 13 '17 at 09:19
  • His is a special case of the above with: $s=1$, $a=n$, $b = n - \beta$, and $\beta < 0$. That is, $n$ is the small $x$ power law index, $\beta$ is the change in the power law index across $x=x_n$ ($x=1$ in Cye's answer), and $s> 0$ is the sharpness (common cases: $s=1$ for simplicity of formula, $s=2$ for a Pythoagorean theorem-like break). – Sean Lake Sep 13 '17 at 17:35
  • The form I gave the power law was specifically designed to accommodate the addition of more breaks by multiplying by another factor with the same form as the one in brackets (including the exponent of the brackets). – Sean Lake Sep 13 '17 at 17:36
  • Yes, I have to admit that it took my a while to get a handle on it, but now that I have, I quite like the flexibility of this form, and its relative simplicity. The ability to add breaks at will is really neat. I suspect I'll come across this kind of data again, and if I do, I'll be glad to have this in my mathematical toolbox. Thanks again. – Alexko Sep 13 '17 at 21:35