Fitting one set of data with two functions

Question

I have a set of data (42 points):

{{1.1*10^6, 1/21}, {600., 3/7}, {80., 5/7}, {100., 29/42}, {600., 3/
  7}, {600., 3/7}, {70., 31/42}, {500., 19/42}, {25., 37/42}, {25., 
  37/42}, {300., 4/7}, {25., 37/42}, {15000., 5/42}, {150., 2/
  3}, {1400., 13/42}, {10., 1}, {60., 16/21}, {200., 13/21}, {10000., 
  4/21}, {10000., 4/21}, {1000., 5/14}, {50., 17/21}, {1500., 2/
  7}, {200., 13/21}, {300., 4/7}, {5000., 3/14}, {10., 1}, {80000., 1/
  14}, {20., 13/14}, {1575., 11/42}, {400., 1/2}, {1000., 5/
  14}, {9.1671*10^6, 1/42}, {60000., 2/21}, {20., 13/14}, {50., 17/
  21}, {400., 1/2}, {150., 2/3}, {10., 1}, {10000., 4/21}, {350., 11/
  21}, {2200., 5/21}}

On a log-log plot they look like this:

The x coordinate represents area (of forest fire, in m^2), and the y coordinate represents the number of fires that are equal or greater than the value on x-coordinate, divided by total numbers of fires (empirical cumulative distribution). Some of the fires have the same burned area and that is the reason why doubled data are present.

The first part of the data fit well with a so-called "stretched exponential":

$$y=\exp(-(x/A)^\beta),$$

where $A$ and $\beta$ are constants that I must find.

The second part of the data fit better to a Power law

$$y=Cx^{-\alpha},$$

where $C$ and $\alpha$ are also constants that I must find.

My question is: how to make Mathematica find a best fit for these data and for both functions at the same time? I want to get an x-value X* that is on the "borderline" between these two functions, all constants that are present in the functions, and a curve that is smooth but consists of two parts: a left part ("before" X*) should be a part of the stretched exponential (curve in this graph) and a right part ("after" X*) should be a plot of the Power law (straight line in LogLogPlot).

How to do all of that?

---

P.S. Short explanation of why I need this:

In a forest fire (and mathematical - physical similar problems, like spread of oil on sea for example), if the fire behavior is like a stretched exponential (where are the "A" values, which represents area) then the fire is "under control" until (approximately) and area of that size. If it is similar to a Power law, which has no scale, than small fire behavior is similar to a big one (and vice versa) and fire behavior is independent of human activities, it depends only on "outside" factors like (big) rain or natural borders (like water, big slope, bald rock etc.).

Also, if the functions continue smoothly from one to each other, the ordinate of X* (the Y* value associated to X*) is the percentage of total numbers of fires whose behavior is "scale free"; X* is then the burned area after which the fire is "totally out of control."

Thanks for any help or suggestion!

Answer 1

One can define a function with the characteristics that you want. The curves meet at $x_0$ and the requirement of a smooth transition is characterized by having the first derivatives of the two curves being equal.

To make things easier we deal with the logs of the two functions. We need the parameters to satisfy the following two equalities:

$$-\left(\frac{x_0}{a}\right)^{\beta }=\log (c)-\alpha \log (x_0)$$

and

$$\frac{\partial \left(-\left(\frac{x_0}{a}\right)^{\beta }\right)}{\partial x_0}=\frac{\partial (\log (c)-\alpha \log (x_0))}{\partial x_0}$$

To meet those requirements forces the above restrictions on the parameters. Below I use $x_0$, $a$, and $\beta$ with $c$ and $\alpha$ being functions of $x_0$, $a$, and $\beta$. This results in the following:

$$c=e^{\left(\frac{x_0}{a}\right)^{\beta } (\beta \log (x_0)-1)}$$

$$\alpha=\beta \left(\frac{x_0}{a}\right)^{\beta }$$

I used Mathematica (as you used that package as a tag). Here is the code for finding $c$ and $\alpha$ in terms of the other parameters:

{c, α} = {c, α} /. Solve[{-(x0/a)^β == Log[c] - α  Log[x0],
     D[-(x0/a)^β, x0] == D[Log[c] - α  Log[x0], x0]}, {c, α}][[1]]

Below is a complete set of code for solving for the parameter restrictions and fitting the piecewise function.

(* Available data *)
data = {{1.1*10^6, 1/21}, {600., 3/7}, {80., 5/7}, {100., 
    29/42}, {600., 3/7}, {600., 3/7}, {70., 31/42}, {500., 
    19/42}, {25., 37/42}, {25., 37/42}, {300., 4/7}, {25., 
    37/42}, {15000., 5/42}, {150., 2/3}, {1400., 13/42}, {10., 
    1}, {60., 16/21}, {200., 13/21}, {10000., 4/21}, {10000., 
    4/21}, {1000., 5/14}, {50., 17/21}, {1500., 2/7}, {200., 
    13/21}, {300., 4/7}, {5000., 3/14}, {10., 1}, {80000., 
    1/14}, {20., 13/14}, {1575., 11/42}, {400., 1/2}, {1000., 
    5/14}, {9.1671*10^6, 1/42}, {60000., 2/21}, {20., 13/14}, {50., 
    17/21}, {400., 1/2}, {150., 2/3}, {10., 1}, {10000., 4/21}, {350.,
     11/21}, {2200., 5/21}};

(* Define logs of functions *)
logf = -(x/a)^β;
logg = Log[c x^(-α)];

(* Solve for where curves meet at x0 such that the first derivatives are equal *)
solc = Solve[logf == logg, c, Reals][[1]] /. x -> x0;
solα = Solve[(D[logf, x] /. x -> x0) == (D[logg, x] /. x -> x0 /. solc), α][[1]];

(* Define piecewise equation to fit by plugging in the required values of c and α *)
eqn = Piecewise[{{logf, x <= x0}, {logg /. solc /. solα, x > x0}}];

(* Take the logs of the dependent variable *)
data2 = Transpose[{data[[All, 1]], Log[data[[All, 2]]]}];

(* Fit the model *)
nlm = NonlinearModelFit[data2, eqn, {{x0, 200}, {a, 500}, {β, 1}}, x];
nlm["BestFitParameters"]
(* {x0 -> 144.7452999922721,a -> 466.25976118614983,β -> 0.7237450948802254} *)
α /. solα /. nlm["BestFitParameters"]
(* 0.3103888712712745 *)
c /. solc /. solα /. nlm["BestFitParameters"]
(* 3.0505203130634624 *)

(* Some fit and residual plots *)
Show[ListLogLogPlot[data],
 LogLogPlot[
  Exp[nlm[x]], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}]]

Histogram[nlm["FitResiduals"]]

ListPlot[Transpose[{nlm["PredictedResponse"], nlm["FitResiduals"]}],
 AxesOrigin -> {-4, 0}]

QuantilePlot[nlm["FitResiduals"]]

The residual plots don't look great but you'll need to decide about that.

Of course, standard errors for the parameters and confidence bands can be constructed (although confidence bands for piecewise regressions do have complications).