2
$\begingroup$

I have a system that transit between two different states. Each state output varies linearly with time, given by $m t+c$, where both lines intersect the $x$ axis in the same point. The output of this system is something like this:

Ramp

The idea is two obtain the two lines that describe this system. What I tried first is guess two lines and assign the points in the data to be fitted to this lines depending on how far they fall.

func1[x_] := 0.06*x - 300
func2[x_] := 0.002*x - 10

squared1 = Table[{x, (func1[x] - data[[x]])^2}, {x, 1, 
Length[data]}];
squared2 = Table[{x, (func2[x] - data[[x]])^2}, {x, 1, 
Length[data]}];

The squared function is giving me the squared error of each data point to each line. Then I create two list to which I assign the values to be fitted latter.

list1 = {};
list2 = {};
Table[If[squared1[[x, 2]] > squared2[[x, 2]], AppendTo[list2, {x,data[[x]]]}, AppendTo[list1, {x,data[[x]]}]], {x, Length[data]}];

In this way I created two list with the values that are closer to each of the guesses.

Then I fit each list with a line:

line1 = LinearModelFit[list1, x, x];
line2 = LinearModelFit[list2, x, x];

The idea will be to repeat this process until the results converge, but they don`t.

This is the result from the first fit:

first fit

The second is worst:

second fit

And they don't converge. I don`t believe is a problem with choosing the initial values but of convergence due to the assignation of the lists.

This question is related to ion channel biophysics. Single channels vary between open and closed states. The current that goes trough the system depends on the state of the channel and the voltage. If we observe a single channel while we change the voltage we'll see something like the image.

I would like to generalize this later in three ways:

1) One of the lines will be a parabola

2) I could get several time series like the one in the graph, with the idea of fitting all of them to the same set of parameters.

3) Having more than two lines.

Help in this secondary items is not necessary but appreciated.

The dataset file from which I constructed the graph is here: data

Or can me imported

data = Import["https://pastebin.com/raw/z321TFZ6"]

I did some more data with Gaussian noise: Gauss_noise

$\endgroup$
  • 1
    $\begingroup$ @rhermans Data is at the end of the post. $\endgroup$ – BPinto Aug 9 '18 at 21:48
  • 1
    $\begingroup$ There are some other things you might want to consider: (1) How to use the information in the random process (or deterministic process) for shifting from one curve to another, (2) For real-life data would you expect that the variances about the lines to be equal? and (3) Again I have to assume you've presented simulated data with a uniform (rather than Gaussian) error structure so that might affect what goodness-of-fit metric you might want to use to estimate the parameters. In other words, this is not just a data-only exercise but you really need to define the model that generates the data. $\endgroup$ – JimB Aug 10 '18 at 16:07
  • $\begingroup$ @JimB The data is expected to have the same Gaussian noise for both lines. The underlying process is a Markov chain that transit between two states. In order to to include this in the adjustment we will need to make assumptions about the underling model or make a fit, which is not my intention at this point. $\endgroup$ – BPinto Aug 10 '18 at 19:17
  • $\begingroup$ The noise is clearly not Gaussian. It appears to be uniformly distributed. Is that intentional? $\endgroup$ – mef Aug 16 '18 at 9:45
  • $\begingroup$ @mef you're right, in this case the noise is given by RandomReal. In the data I expect to analyze is gaussian. $\endgroup$ – BPinto Aug 17 '18 at 11:35
3
$\begingroup$

Concept

NMinimize for the Min Mean distance of the data to the $n$ fitting functions, all equally weighted.

Solution

data = Import["https://pastebin.com/raw/z321TFZ6"]

func[parms_List][x_] := 
 parms[[1]] + Sum[parms[[k]] x^(k - 1), {k, 2, Length[parms]}]
dist[f_][{x_, y_}] := Abs[y - f[x]]  (* A distance function *)
minPos[l_List] := {#, Position[l, #]} &[Min[l]] (* Min and it's pos*)

sol = Last@With[{n = 2, m = 2},
   NMinimize[
    Total[
     Mean /@ GatherBy[
        Map[
         Function[{L}, 
          minPos@Table[dist[func[Array[c[#, k] &, n]]][L], {k, m}]
         ]
         , data
        ]
        , Last][[All, All, 1]]
     ]
    , Flatten[Table[c[i, j], {i, n}, {j, m}]]]]

(* {c[1, 1] -> -57.221, c[1, 2] -> -330.999, c[2, 1] -> 0.0110719, 
 c[2, 2] -> 0.0663247} *)

Show[ListPlot[
  GatherBy[data, 
   Position[#, Min[#]] &[
     Table[dist[func[{c[1, k], c[2, k]} /. sol]][#], {k, 2}]] &]
  , PlotTheme -> "Scientific"
  , PlotStyle -> {Pink, Gray}
  ],
 Plot[
  Evaluate[
   Table[func[{c[1, k], c[2, k]}][x], {k, 2}] /. sol
   ]
  , {x, 0, 10000}
  , PlotStyle -> {Red, Blue}
  ]
 ]

Mathematica graphics

Explanations

I use some auxiliary functions

func[parms_List][x_] Function to fit. In this case a polynomial of arbitrary order given by the number of coefficients.

func[{a, b, c}][x]
(* a + b x + c x^2 *)

dist[f_][{x_, y_}] := Abs[y - f[x]] A distance function. You can play with other distances.

minPos[l_List] := {#, Position[l, #]} &[Min[l]] Gives the Min and it's Position

Other details:

Table[dist[func[Array[c[#, k] &, n]]][L], {k, m}] Creates the list of functions to fit.

With[{n = 2, m = 2}… Here n=2 is the number of coefficients for the polynomial (linear fit) , m=2 the number of polynomials to fit.

Mean /@ GatherBy[ … minPos@ Calculates the minimum distance to all the available functions and GatherBy the position of the Min, to then takes the Mean and then the Total so all curves weight the same, regardless the number of points in each.

Analysis

Conceptually there are many similarities with the answer by @bills and with my answer to this question.

The main difference with @bills's is that his requires an initial guess and separates the data first. I NMinimize without previous separation. I also give more flexibility to modify the fitting function and the number of curves to fit.

The main difference with my other answer is that here I solve the problem generated when the two models have a significantly different number of points. That is done by using Total[Mean/@ GatherBy[… instead of just Total therefore giving equal weight to each of the curves been fit.

A fundamental problem persist in the intersection where the points get skewed, probably a better approach would be to allow points to belong to many curves, but for that we would need to define a threshold. This cause the fit with the fewer points to be tilted. @bills's fit is better, but it depends on a good initial guess.

For convergence, there is room to play with different definitions for the distance function dist and different NMinimize Method->.

In terms of flexibility, it's easy to define other fitting functions func and scale up the numbers of functions to fit.

$\endgroup$
  • $\begingroup$ This really looks like what I'm looking for. Specially I think I could add the constrain of both curves intersecting x in the same point. I also don't see the skew when I run it. $\endgroup$ – BPinto Aug 10 '18 at 2:34
3
$\begingroup$

Here is an approach that provides guesses for the initial values of the parameters given that the functions are lines and that both intersect the y axis at the same point:

(* Define distance from observed point to the closest line *)
f[y_, x_, x0_, b1_, b2] := 
 Min[(y - b1 (x - x0))^2, (y - b2 (x - x0))^2]

(* Get initial values *)
x0Init = Sort[data, Abs[#1[[2]]] < Abs[#2[[2]]] &][[1, 1]];
sortedData = Sort[data, Abs[#1[[1]]] < Abs[#2[[1]]] &];
nLower = Range[Floor[0.05 Length[sortedData]]];
lowerData = sortedData[[nLower]];
nUpper = Range[Floor[0.95 Length[sortedData]], Length[sortedData]];
upperData = sortedData[[nUpper]];
lowerMinMax = MinMax[lowerData[[All, 2]]];
upperMinMax = MinMax[upperData[[All, 2]]];
b1Init = (upperMinMax[[2]] - lowerMinMax[[1]])/(Mean[nUpper] - Mean[nLower]);
b2Init = (upperMinMax[[1]] - lowerMinMax[[2]])/(Mean[nUpper] - Mean[nLower]);

(* Find solution *)
sol = FindMinimum[Mean[f[#[[2]], #[[1]], x0, b1, b2] & /@ data],
  {{x0, x0Init}, {b1, b1Init}, {b2, b2Init}}]
(* {499.183, {x0 -> 4998.07, b1 -> 0.0662429, b2 -> 0.00642343}} *)

(* Show results *)
Show[ListPlot[data],
 Plot[{b1 (x - x0), b2 (x - x0)} /. sol[[2]],
  {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}, PlotStyle -> Red]]

Data and fit to two lines

I'm sure that better estimates of the initial values can be obtained but this approach works for the supplied data.

Now while the fits "look" OK, one really needs to be able to provide standard errors to the estimates of the parameters. Probably a parametric bootstrap would be appropriate.

$\endgroup$
1
$\begingroup$

Your list1 and list2 are not segmenting the data properly. Here is one way to fix this:

func1[x_] := 0.06*x - 300
func2[x_] := 0.002*x - 10
data1 = Select[
   Table[If[(func1[i] - data[[i, 2]])^2 < (func2[i] - data[[i, 2]])^2,
      data[[i]]], {i, Length[data]}], # =!= Null &];
data2 = Select[
   Table[If[(func1[i] - data[[i, 2]])^2 >= (func2[i] - 
      data[[i, 2]])^2, data[[i]]], {i, Length[data]}], # =!= Null &];
line1 = LinearModelFit[data1, x, x];
line2 = LinearModelFit[data2, x, x];
Show[ListPlot[{data1, data2}], 
 Plot[{line1[x], line2[x]}, {x, 1, Length[data]}]]

enter image description here

$\endgroup$
  • $\begingroup$ What is this 55? $\endgroup$ – rhermans Aug 9 '18 at 22:38
  • $\begingroup$ Instead of fitting to two lines, I select the data as all the horizontal ones (which match the first line) and all the angled ones (which do not match). I selected the value by looking at the plot. This seemed simpler than selecting all the points closer to func1 and to func2. $\endgroup$ – bill s Aug 9 '18 at 22:41
  • $\begingroup$ @bills This is better than what I did, but you can see that the cutoff is a little bit rough. I can observe some yellow near the blue line, and the slope of the blue line doesn`t seems right. Maybe using as a selection a third line in between instead of an absolute value may do better. $\endgroup$ – BPinto Aug 9 '18 at 22:50
  • $\begingroup$ I've revised the answer, using the original formulation of closeness to the two functions func1 and func2. $\endgroup$ – bill s Aug 9 '18 at 22:55
  • $\begingroup$ I don't understand, you start knowing the functions in order to separate two groups and then fit to get the same functions. Seems circular. You can do without that, look my answer here: mathematica.stackexchange.com/a/146105/10397 $\endgroup$ – rhermans Aug 9 '18 at 23:01

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.