2
$\begingroup$

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:

enter image description here

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!

$\endgroup$
  • 1
    $\begingroup$ There clearly is a challenging set of statistical questions here (+1). $\endgroup$ – whuber Aug 11 '15 at 19:31
0
$\begingroup$

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 meet those requirements forces restrictions on the parameters. Below I use $x_0$, $a$, and $\beta$ with $c$ and $\alpha$ being functions of $x_0$, $a$, and $\beta$. And while this can be performed in R or SAS (or many other packages), I used Mathematica (as you used that package as a tag).

(* 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)^\[Beta];
logg = Log[c x^(-\[Alpha])];

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

(* Define piecewise equation to fit by plugging in the required values of c and \[Alpha] *)
eqn = Piecewise[{{logf, x <= x0}, {logg /. solc /. sol\[Alpha], 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}, {\[Beta], 1}}, x];
nlm["BestFitParameters"]
(* {x0 -> 144.7452999922721,a -> 466.25976118614983,\[Beta] -> 0.7237450948802254} *)
\[Alpha] /. sol\[Alpha] /. nlm["BestFitParameters"]
(* 0.3103888712712745 *)
c /. solc /. sol\[Alpha] /. 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]]]}]]

Data and fit

Histogram[nlm["FitResiduals"]]

Histogram of residuals

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

Predicted vs residuals

QuantilePlot[nlm["FitResiduals"]]

qq plot of residuals

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.