How can I fit a set of points to the form $y = a(b-x)^c$?

Question

I have a list of points $(x, y)$, and I expect them to be related as $$y = a(b-x)^c,$$ where $a$, $b$, and $c$ are unknown constants. I also know the approximate values of constants, which can be used as initial guess.

How can I fit the curve to my data?

I tried using NonlinearModelFit, but it returns errors. Here, the actual values (which I am trying to estimate) are $a = 1, b = 1, c = 0.5$.

data = N[Table[{n/50, Sqrt[1 - n/50] + RandomReal[]/10}, {n, 1, 49}]]
(* {{0.02, 0.989949}, {0.04, 0.979796}, {0.06, 0.969536}, {0.08, 
  0.959166}, {0.1, 0.948683}, {0.12, 0.938083}, {0.14, 
  0.927362}, {0.16, 0.916515}, {0.18, 0.905539}, {0.2, 
  0.894427}, {0.22, 0.883176}, {0.24, 0.87178}, {0.26, 
  0.860233}, {0.28, 0.848528}, {0.3, 0.83666}, {0.32, 
  0.824621}, {0.34, 0.812404}, {0.36, 0.8}, {0.38, 0.787401}, {0.4, 
  0.774597}, {0.42, 0.761577}, {0.44, 0.748331}, {0.46, 
  0.734847}, {0.48, 0.72111}, {0.5, 0.707107}, {0.52, 0.69282}, {0.54,
   0.678233}, {0.56, 0.663325}, {0.58, 0.648074}, {0.6, 
  0.632456}, {0.62, 0.616441}, {0.64, 0.6}, {0.66, 0.583095}, {0.68, 
  0.565685}, {0.7, 0.547723}, {0.72, 0.52915}, {0.74, 
  0.509902}, {0.76, 0.489898}, {0.78, 0.469042}, {0.8, 
  0.447214}, {0.82, 0.424264}, {0.84, 0.4}, {0.86, 0.374166}, {0.88, 
  0.34641}, {0.9, 0.316228}, {0.92, 0.282843}, {0.94, 
  0.244949}, {0.96, 0.2}, {0.98, 0.141421}} *)

fit = NonlinearModelFit[data, {a (b - q)^c}, {{a, 0.8}, {b, 0.9}, {c, 0.6}}, q, 
  Method -> NMinimize]

NonlinearModelFit::nrnum: The function value 1.50013 -0.120551 I is not a real number at {a,b,c} = {0.8,0.9,0.6}.

NMinimize::nnum: The function value Experimental`NumericalFunction[{Hold[Block[{q={0.02,0.04,0.06,0.08,0.1,0.12,0.14,0.16,0.18,0.2,<<39>>},Optimization`FindFit`y$1610={0.989949,0.979796,0.969536,0.959166,0.948683,0.938083,0.927362,0.916515,0.905539,0.894427,<<39>>}},1/2 (Times[<<2>>]-Optimization`FindFit`y\$1610).(Times[<<2>>]-Optimization`FindFit`y\$1610)]],Block},{0,{{{},1,0,Hold[0.8],0,0},{{},1,1,Hold[0.9],0,0},{{},1,2,Hold[0.6],0,0}}},<<3>>,{None,None,None}][{0.8,0.9,0.6}] is not a number at {a,b,c} = {0.8,0.9,0.6}.

Answer 1

Method-1

• Remove the restriction of a and b.Also remove the Method -> "NMinimize.

data = N[Table[{n/50, Sqrt[1 - n/50] + RandomReal[]/10}, {n, 1, 49}]];
nlm = NonlinearModelFit[data, {a (b - x)^c}, {a, b, {c, .6}}, x]
Plot[nlm[x], {x, 0, 1}, Epilog -> {Red, Point[data]}, 
 AspectRatio -> Automatic]

Method-2

• If we want the c near .5, the initial value of c may set to be (c - .6) etc.

data = N[Table[{n/50, Sqrt[1 - n/50] + RandomReal[]/10}, {n, 1, 49}]];
nlm = NonlinearModelFit[data, {a (b - x)^(c - .6)}, {a, b, c}, x]

Plot[nlm[x], {x, 0, 1}, Epilog -> {Red, Point[data]}, 
 AspectRatio -> Automatic]

Method-3

• Or add a perturbation of the exponent.

Clear["Global`*"];
data = N[
   Table[{n/50, Sqrt[1 - n/50] + RandomReal[]/10}, {n, 1, 49}]];
nlm = NonlinearModelFit[data, {a (b - x)^(c + d)}, {a, b, c, d}, x]
Plot[nlm[x], {x, 0, 1}, Epilog -> {Red, Point[data]}, 
 AspectRatio -> Automatic]

Answer 2

If appropriate to better understand the functional form near x=1 just fit that portion of the data.

fit = NonlinearModelFit[data[[30 ;;]], {a (b - x)^c }, {a, b, c}, x];
Plot[fit[x], {x, 0.02, 0.98}, Epilog -> Point[data]]

However the fit is still not the best.

Altering the parameterization greatly improves the fit.

fit = NonlinearModelFit[data[[30 ;;]], {a (b - x)^c - d}, {a, b, c, d}, x];
Plot[fit[x], {x, 0.02, 0.98}, Epilog -> Point[data]]

Answer 3

For some values of a,b,c y may be complex. MMA gives you a warning but gets nevertheless a fit:

MMA version 13.3 Windows 10

fit = NonlinearModelFit[data, {a (b - x)^c}, {a, b, c}, x];
Plot[fit[x], {x, 0.02, 0.98}, Epilog -> Point[data]]

Answer 4

This gives a result similar to some of the others posted. The idea is to linearize the objective function by taking logs. This means sweeping b and the x values inside a log together so we create a function for that purpose and use it in the defining objective.

I had to use NMinimize because FindMinimum, at least with default settings, was unwilling to move from the initial values.

logx[b_, vals_] := Log[b - vals]

bestFit[data_] := Module[
  {xvals, yvals, obj, ly, y},
  {xvals, yvals} = Transpose[data];
  obj[a_Real, b_Real, c_Real] := 
   Module[{ovec = 
      Log[yvals] - (Log[a] + c*logx[b, xvals]) /. 
       Indeterminate :> Nothing, res},
    Abs[ovec . ovec]];
  NMinimize[obj[a, b, c], {{a, .7, .9}, {b, .8, .9}, {c, .6, .7}}]
  ]

data = N[Table[{n/50, Sqrt[1 - n/50] + RandomReal[]/10}, {n, 1, 49}]];

Let's try it out.

In[287]:= {min, vals} = bestFit[data]

(* Out[287]= {0.0842183, {a -> 1.0391, b -> 1.02703, c -> 0.515608}} *)

Plot[a*(b - x)^c /. vals, {x, 0.02, 0.98}, Epilog -> Point[data]]

Answer 5

The cause for the errors is that you give a starting value for b that is too small. You need to give a starting value for b larger than Max[data[[All,1]] otherwise you end up with complex numbers. It also doesn't hurt to put in that explicit condition.

SeedRandom[12345];
data = N[Table[{n/50, Sqrt[1 - n/50] + RandomReal[]/10}, {n, 1, 49}]]

fit = NonlinearModelFit[data, {a (b - q)^c, b > Max[data[[All, 1]]]}, 
  {{a, 0.8}, {b, Max[data[[All, 1]]]}, {c, 0.6}}, q]

Answer 6

A minor variant:

d = N[Table[{n/50, Sqrt[1 - n/50] + RandomReal[]/10}, {n, 1, 49}]];
td = {#1, Log@#2} & @@@ d;
n = NonlinearModelFit[td, a + c Log[b - x], {a, b, c}, x];
Plot[Exp[n[x]], {x, 0.01, 1}, Epilog -> Point[d]]
n["ParameterTable"]
model = Exp[a] (b-x)^c /. n["BestFitParameters"]