calculus and analysis - Taking log derivatives of data in Mathematica

We have two lists (vectors) of data, y and x, we can imagine x being time steps (0,1,2,...) and y some system property computed at each value of x. I'm interested in calculating the derivative of log of y with respect to log of x, and the question is how can I preform such calculations in Mathematica? We first go to log space, e.g. logx = Map[Log,x] similarly for y. Then what method do we use for the differentiation dlog(y)/dlog(x)? and then eventually we'll plot it as a function of log(x).

Context: When we have a quantity y for which we expect to see a power law behaviour as a function of x, so $y \propto x^{\alpha},$ taking the log derivative of y w.r.t to x allows us to extract the actual $\alpha$ corresponding to our data.

Answer

In principle, I like most of all the approach proposed by @Daniel Lichtblau as a comment. I would propose Daniel to formulate it as a regular answer.

If the data is noisy, the interpolation can give rise to problems, since it goes through all points. The function, thus, strongly oscillates, and its derivatives can appear unnecessarily large and also noisy.

Instead, the strategy using fitting may be a good idea to apply, if the noise is present. The latter enables one to create a smooth, slowly varying function, thus effectivey reducing the effect of the noise in the derivative. Here I would like to propose a simple idea of how to treat the problem by the fitting. Since the author gave no list, to be able to show code, I will create such a list. Since the author is primarily interested in a power law-type of dependence, this will be the cases of the lists below.

The first example is a simple power law to which I add a bit of noise:

lst1 = Table[{x, x^(3/2) + RandomReal[{-30, 30}]}, {x, 1, 100}];
ListPlot[lst1]

yielding this:

This list can be transformed into a Log-Log list:

lst1Log = lst1 /. {x_, y_} -> {Log[x], Log[y]} // N;

The next step is to fit it to a suitable function. In principle, by looking at its plot one can easily guess that at large values of t=Log[x] it has a linear asymptote. This agrees with the author's expectation of a power law. It is not astonishing, since I created this list to be such. However, the authors expectation and the form of the plot gives me a hint to fit it to a linear function:

     model1 = a + b*t;
ff1 = FindFit[Select[lst1Log, #[[1]] > 3 &], model1, {a, b}, t]
Show[{
  ListPlot[lst1Log, PlotStyle -> Blue],
  Plot[model1 /. ff1, {t, 0, 5}, PlotStyle -> Red]
  }]
    (*   {a -> 0.167, b -> 1.46}   *)

where a and b are the fitting parameters, bbeing the exponent of the power law. The result differs from the power law b=1.5 since it is influenced by stattered points at small t. To improve it one can cut away the list points at smaller t: Select[lst1Log, #[[1]] > 4 &] and make a fit starting from the value b=1.46 obtained within the previous step:

model1 = a + b*t; ff1 = FindFit[Select[lst1Log, #[1] > 4 &], model1, {a, {b, 1.46}}, t] Show[{ ListPlot[lst1Log, PlotStyle -> Blue], Plot[model1 /. ff1, {t, 0, 5}, PlotStyle -> Red] }]

(* {a -> -0.0245, b -> 1.51} *)

Here how it looks like:

The difference from 1.5 is due to the noise.

The second example is the list containing two power laws and exhibiting a slightly more complex dependence:

    lst2 = Table[{x, 5 x^(1/2) + 0.001*x^2.7 + RandomReal[{-10, 10}]}, {x,
         4, 1000}];
    lst2Log = lst2 /. {x_, y_} -> {Log[x], Log[y]} // N; 
ListPlot[lst2]

It is the same story, but looking at the plot of the Log-Log list (that is, lst2Log) one can guess that it can be fitted by a slightly different function:

model2 = (a + b*t^(d + 1))/(1 + c*t^d)

Here a, b, c and d are parameters. In such functions, however, the exponents (such as d) are usually difficult to fit. For this reason instead of the fitting it is convenient to manipulate with the parameter dwhile fitting all the rest:

Manipulate[
 model2 = (a + b*t^(d + 1))/(1 + c*t^d);
 ff2 = FindFit[lst2Log, model2, {{a, 3}, b, c}, t];
 Show[{

   ListPlot[lst2Log, PlotStyle -> {Blue, PointSize[0.008]}],
   Plot[model2 /. ff2, {t, 0, 7}, PlotStyle -> Red]
   }], 
{{d, 3.58}, 2, 5, Appearance -> "Labeled"}]

By manipulating I found that the value d=3.58 gives a good approximation. It is this value that I placed as an initial one for convenience. This is what you should see on the screen:

Now one can get the approximation more precisely, by using the values of the parameters obtained in the previous step as their initial values:

model2 = (a + b*t^(4.6 + eps))/(1 + c*t^(3.6 + eps));
ff2 = FindFit[Select[lst2Log, #[[1]] > 3.5 &], 

  model2, {{a, 2.55}, {b, 0.007}, {c, 0.003}, {eps, 0.2}}, t]
Show[{
  ListPlot[Select[lst2Log, #[[1]] > 3.5 &], 
   PlotStyle -> {Blue, PointSize[0.008]}],
  Plot[model2 /. ff2, {t, 0, 7}, PlotStyle -> Red]
  }]

(*  {a -> 2.72, b -> 0.00462, c -> 0.00228, eps -> 0.258}  *)

The exponent in infinity will in this case be b/c. However you can here find the exponent value in each point by calculating the derivative of the model2in the desired point.

Generally, in this approach the fitting gives you a simple smooth approximation of the desired function, the derivative of which you can further find.

Have fun!

Edit 1: To illustrate how to find the regimes. Within the example 2 you can define the derivative of ydenoted below as der[t] and plot it:

Clear[der];
der[t_] := Evaluate[D[model2 /. ff2, t] // Simplify];

Plot[der[t], {t, 0, 7}, PlotStyle -> Red, 
 AxesLabel -> {Style["t=Log[x]", Italic, 12], 
   Style["der=\!\(\*FractionBox[\(dLog[y]\), \(dt\)]\)", Italic, 12]}]

Here you will find 3 regimes, from 1 to several decades each. In the middle a regime with a derivative described by an almost straight tilted line is visible. It is flanked by two power law rerimes at small and large t values.

Edit 2: To address your question with your list.

In this case I would do as follows. Here is your list, that should be somehow downloaded

lst = Drop[
Import[FullPathToYourList, "Table"], 1] ;

I have dropped the first element of the list, since there is something strange in it. But it makes a negligible influence on the result. and lst1=lst/. {x_, y_} -> {Log[x], Log[y]} Now, looking at the ListPlot[lst1] one sees that it consists of two straight regions corresponding to two different regimes (i) at Log[x]<4 and (ii) at Log[x]>4. In principle this means that you could express it in the form of a formula y=a e^{-\frac{x}{d}} x^{\alpha }+b \left(1-e^{-\frac{x}{d}}\right) x^{\beta } (I hope that the TeX expression here works). Here d is a cut-off close to 10^4. However, it is not easy to fit such a formula and to do it on the huge interval that takes place here.

For this reason it is more easy to fit the log-log list lst1, and to do it for t=Log[x]<4´separately from that for t>4. Below I do it:

   model1 = a + b*x;

model2 = c + k*x;

ff1 = FindFit[Select[lst, #[[1]] < 4 &], model1, {a, b}, x]
ff2 = FindFit[Select[lst, #[[1]] > 10 &], model2, {c, k}, x]
Show[{
  ListPlot[lst, PlotStyle -> Blue, 
   AxesLabel -> {Style["t=Log[x]", 12, Italic], 
     Style["Log[y]", 12, Italic]}],
  Plot[model1 /. ff1, {x, 1, 4}, PlotStyle -> Red],
  Plot[model2 /. ff2, {x, 4, 14}, PlotStyle -> Green]

  }]

(*  {a -> 1.54653449015407020, b -> 0.760890314643084274}

{c -> 3.002027525233422865, k -> 0.3272780474976788419}   *)

On the screen you will see the following plot:

Now, the fitting parameters b and k above are the exponents you are looking for. So I would say that the exponents are 0.76 and 0.33.

Pay attention that in the ff2 I fit the part of the list with t>10. It is because the interval 4 looks as a transition part of the curve. It is for this reason the solid green line deviates from the blue points in this interval. This is, however, up to you to choose how and where it is better to fit. It strongly depends on the nature of the problem.

Have fun!