Nonlinearmodelfit with conditionals - how to.

Questions
,
1

Hey guys I am trying to adapt a SAS code to mathematica, but I am experiencing some problems in doing that. I tried to run a model, but wasn’t sure how to insert some conditionals.

SAS code:

proc nlin method = newton;
parms a=1 b=10^6 d=2.8;
If x > a then do;
Model S=0;
end;
else
if x < b then do;
model S =1;
end;
else
model S = (x^(d - 3) - b^(d - 3))/(a^(d - 3) - b^(d - 3));
end;
run;

My first mathematica model, without using the conditionals (x > a;x < b)

data = {{0.1 , 1},{0.5 , 0.99245283},{1 , 0.981132075},{2 , 0.957629013},{3 , 0.94647978},{5 , 0.892872862},{8 , 0.825362586},{10 , 0.794884307},{15 , 0.737820069},{30 , 0.688456089},{50 , 0.710852149},{100 , 0.686350949},{300 , 0.630188679},{500 , 0.622243133},{1500 , 0.61197978}};

model = (x^(d - 3) - b^(d - 3))/(a^(d - 3) - b^(d - 3));

nlm = NonlinearModelFit[data,model, {{a, 1}, {b, 10^6}, {d, 2.8}}, x,Method -> "Newton"]
2

The following gives a result with a warning:

Clear[x, a, b, d];
data = {{0.1, 1}, {0.5, 0.99245283}, {1, 0.981132075}, {2, 
    0.957629013}, {3, 0.94647978}, {5, 0.892872862}, {8, 
    0.825362586}, {10, 0.794884307}, {15, 0.737820069}, {30, 
    0.688456089}, {50, 0.710852149}, {100, 0.686350949}, {300, 
    0.630188679}, {500, 0.622243133}, {1500, 0.61197978}};
model = Piecewise[{{0, x > a}, {1, 
     x < b}, {(x^(d - 3) - b^(d - 3))/(a^(d - 3) - b^(d - 3)), True}}];
nlm = NonlinearModelFit[data, model, {{a, 1}, {b, 10^6}, {d, 2.8}}, x]
3

Are you sure your conditions are x>a or x<b, and not a>b or x<a? That would make more sense.

4

So. The case is that S is the relative saturation, b is the moment in which the saturation is 0 and a is the moment in which the saturation is 1, so that if x(that is the potential in a specific water content)>b saturation is equal to 0, and if the potential in a specific water content (x) is minor than a, so that the saturation is 1.

I think your code has work succesfully, I only changed the way you wrote the conditionals.

5

I am not sure if a explained it in an understandable way. Well I tried to use this model by using LevenbergMarquardt method and it worked, wasn’t the same when tried Newton method. Some specific reason why it didn’t work?

Thank you anyway.

6

I am no expert. It may have something to do with the fact that your model is not differentiable. Sorry I can’t help you any further. Here is a code that gives an output, even though with a warning:

data = {{0.1, 1}, {0.5, 0.99245283}, {1, 0.981132075}, {2, 
    0.957629013}, {3, 0.94647978}, {5, 0.892872862}, {8, 
    0.825362586}, {10, 0.794884307}, {15, 0.737820069}, {30, 
    0.688456089}, {50, 0.710852149}, {100, 0.686350949}, {300, 
    0.630188679}, {500, 0.622243133}, {1500, 0.61197978}};
model = Piecewise[{{1, x < a}, {0, 
     x > b}, {(x^(d - 3) - b^(d - 3))/(a^(d - 3) - b^(d - 3)), True}}];
nlm = NonlinearModelFit[data, model, {{a, 1}, {b, 10^6}, {d, 2.8}}, x]
7

I see. Well I was trying to use Piecewise function, but as I started working with mathematica last week, I don’t know how to use yet. Thank you again..

8

Your SAS code has several syntax errors. “^” should be “**”, “b=10^6” should be “b=10E6”, and the last “end;” generates an error as it is not needed. Those errors seem a bit more than cut-and-paste issues. Does the code really run in SAS?

9

No. Actually not. I got this code in a scientific paper, so that I am trying to contact the author to see what’s wrong with the code. I am not sure.

10

Following the good suggestion from Gianluca Gorni you’ll find the following result:

piecewiseequation.png
piecewiseequation.png645×100 2.27 KB

Note that the estimate of b is huge and the coefficient -1.02625 x 10^(-11) is essentially zero. That leaves an equation that is linear in the log of x. Given that and that the plot of the fit and the data looks like the following:

dataandfit.png
dataandfit.png564×363 4.49 KB

this screams for taking the logs of both the dependent and independent variables. From the initial description of the issue of wanting a piecewise fit, one possibility is to have the log of the prediction be zero for x < a and b + (d-3)\log(x) otherwise. Because you probably want the fit to be continuous, you’ll want b + (d-3)\log(x) = 0 when x=a. That means b=-(d-3)\log(a). The associated Mathematica code is

model = Piecewise[{{0, x < a}, {-(d - 3) Log[a] + (d - 3) Log[x], True}}];
logY = Log[data[[All, 2]]];
nlm = NonlinearModelFit[Transpose[{data[[All, 1]], logY}], {model, a > 0.1}, {{a, 0.9}, {d, 2.9}}, x];
nlm["BestFitParameters"]
(* {a -> 0.624483, d -> 2.92726} *)
nlm["BestFit"]
Show[ListLogLogPlot[data, PlotRange -> {{0.1, 1600}, {0.5, 1.1}}, PlotStyle -> {PointSize[0.02], Red}, 
  PlotRangeClipping -> False, Frame -> True],
  LogLogPlot[Exp[nlm[x]], {x, 0.1, 1600}]]

piecewiselogequation.png
piecewiselogequation.png569×437 5.46 KB

The fit is not great but the residuals are a bit better behaved than not taking logs.