# Fitting data with interdependence

I have an issue, of which I'm not sure it can be solved, with or without Mathematica. I feel as if it should be possible, but I'm quite clueless as to how.

Some introduction is required. I'm working in a research lab, and we are trying to characterize some samples that were fabricated. These are chips with capacitances on them, three different types. Lets call them $C_1$, $C_2$ and $C_J$. What we then have is 5 devices on which $C_1$ and $C_2$ are constant and $C_J$ is varied, and 4 devices on which $C_J$ is constant and $C_1$ and $C_2$ are varied. However, this is fabrication, we are not sure of the exact values of these parameters. This is what the experiment wants to characterize. They are both being varied by a parameter which we do know: the number of fingers we write on the chip, lets call them $F_{1,2}$ and $F_J$. So the idea is to find the functions $C_{1,2}(F_{1,2})$ and $C_J(F_J)$. Moreover, and this is a little messy, in these measurements it is the case that $C_1$ and $C_2$ always consist of the same amount of fingers, they simply have a slightly different length and thus lead to different capacitance values.

On these devices, measurements are performed. However, we can't measure the capacitance directly; we measure some energy related parameters, $\omega_{1,2}$ and $J$. This is where it gets a bit tricky, because they have a complicated definition, which might be easiest to just give in code:

Cmatrix = ( {
{C1 + CJ + CJ, -CJ, 0, -CJ},
{-CJ, C2 + CJ + CJ, -CJ, 0},
{0, -CJ, C1 + CJ + CJ, -CJ},
{-CJ, 0, -CJ, C2 + CJ + CJ}
} );

InvCmat = Inverse[Cmatrix];

ω1 = Sqrt[1/L InvCmat[[1, 1]]];
ω2 = Sqrt[1/L InvCmat[[2, 2]]];

J = 1/2 InvCmat[[1, 2]]/Sqrt[InvCmat[[1, 1]] InvCmat[[2, 2]]] Sqrt[ω1*ω2];


I should add that in the above L is a known parameter, with value 1.84.

So, what I have is essentially nine datapoints:

{{20, 20, 2, 7.70438, 7.50004, 0.156103}, {20, 20, 3, 7.60618,
7.30572, 0.212588}, {20, 20, 4, 7.46063, 7.18521, 0.273479}, {20,
20, 6, 7.25941, 6.79769, 0.354176}, {20, 20, 10, 6.93511, 6.59471,
0.494317}, {13, 13, 3, 8.8022, 8.20026, 0.297579}, {16, 16, 3,
8.22736, 7.78623, 0.252157}, {20, 20, 3, 7.58361, 7.359,
0.212081}, {24, 24, 3, 7.08458, 6.90351, 0.178109}}


Each each containing (in order of occurance) $F_1$, $F_2$, $F_J$, $\omega_1$, $\omega_2$, $J$. From these datapoints I'd like to find a model for $C_i(F_i)$ and $C_J(F_J)$. Maybe an exact model is not possible, but an approximation should be. Based on physical intuition the relationships should be linear, but it might get a little messy due to fabricational issues.

So, given all this, would anyone have an idea on how to start?

I should also add that the fact that $C_1 \neq C_2$ is an error in the fabrication. Ideally this would not be the case, which simplifies the problem for sure, but even then I'd not know how to do such a fitting routine. So if instead you want to set $C_1 = C_2$ in the initial matrix and describe a method from there it would also be valuable.

Edit: I forgot to add. All values of $C$ are positive and real and $probably$ on the order of $10^{-4}$.

• I upvoted, but perhaps you may get a good insight here – Dr. belisarius Sep 1 '15 at 16:33
• Thanks, I'll give that a try as well. – user129412 Sep 1 '15 at 17:39
• Might be useful to check responses and comments here. – Daniel Lichtblau Sep 1 '15 at 17:42
• @user129412 The way you have the code written for J only Sqrt[InvCmat[[1, 1]] is in the denominator. All the rest of the terms are in the numerator. Can you validate that this was your intention. If you want all of the terms to the right of the / sign to be in the denominator you need to put parenthesis around them. – Jack LaVigne Sep 1 '15 at 19:48
• @JackLaVigne Scratch what I wrote before. The way it is currently written is in fact correct. The two InvCMat terms are in the square root, the omega's do not belong to the denominator. – user129412 Sep 1 '15 at 19:53

You pose a two part question: Part 1 find the "C's" from omega and J:

Cmatrix = ({{C1 + CJ + CJ, -CJ, 0, -CJ}, {-CJ, C2 + CJ + CJ, -CJ,
0}, {0, -CJ, C1 + CJ + CJ, -CJ}, {-CJ, 0, -CJ, C2 + CJ + CJ}});
InvCmat = Inverse[Cmatrix];
L = 1.84
\[Omega]1 = Sqrt[1/L InvCmat[[1, 1]]];
\[Omega]2 = Sqrt[1/L InvCmat[[2, 2]]];
J = 1/2 InvCmat[[1, 2]]/
Sqrt[InvCmat[[1, 1]] InvCmat[[2, 2]]] Sqrt[\[Omega]1*\[Omega]2];
NMinimize[{({\[Omega]1, \[Omega]2, J}
- {7.70438, 7.50004, 0.156103})^2 // Total, C1 > 0, C2 > 0,
CJ > 0}, {C1, C2, CJ}]


{7.89631*10^-31, {C1 -> 0.00841431, C2 -> 0.00892173, CJ -> 0.000386302}}

Here are plots of Ci vs Fi , (i=1,2,J)

• Hm, I see. So you would say use this procedure to find the C's from omega and J, and then perhaps try and see what kind of fit can be done to these C's based on physical knowledge, with the C's and the number of fingers? – user129412 Sep 1 '15 at 19:35
• Its actually simple from there, find your 9 sets of C's and do a separate fit to C1 vs F1, C2 vs F2, etc ( and you are correct they are darn near perfectly linear ) – george2079 Sep 1 '15 at 19:41
• I'll go ahead and confirm this for myself, but that would be great, great news. Oh NMinimize, you're so powerful. – user129412 Sep 1 '15 at 19:45
• And youre right, good linear fits. C2 is showing two anomalies, unless I did something wrong, which is a bit strange considering how consistent the rest is, but that's something for me to think about. Did you also do this fit by any chance, or just for one of the three? Because you mentioned it was a good fit, but I'm assuming you just checked one of them. – user129412 Sep 1 '15 at 20:27
• I added my plots. Looks linear with experimental scatter to me, but of course I don't know anything about the precision you'd expect from your experiment. – george2079 Sep 1 '15 at 20:54

One can use FindMinimum and an objective function to find parameters that minimize the fit between the data and reconstructed data (i.e., forward model) provided one has equations for the forward model.

This will be a bit long winded but here goes.

Step 1. Define the model

We will assume that the capacitor is a linear function of the finger data (you are free to use almost any model).

C1 = c1*F1 + b1

where C1 is the unknown capacitance, F1 is the known finger, and c1 and b1 are the slope and intercept.

We will use a matrix that incorporates this model. We use ReplaceAll to convert your original matrix into the model matrix.

cmodel = ({
{C1 + CJ + CJ, -CJ, 0, -CJ},
{-CJ, C2 + CJ + CJ, -CJ, 0},
{0, -CJ, C1 + CJ + CJ, -CJ},
{-CJ, 0, -CJ, C2 + CJ + CJ}
}) /.
{C1 -> c1*F1 + b1,
C2 -> c2*F2 + b2,
CJ -> cj*Fj + bj}


Which results in

{
{b1 + c1 F1 + 2 (bj + cj Fj), -bj - cj Fj,
0, -bj - cj Fj}, {-bj - cj Fj,
b2 + c2 F2 + 2 (bj + cj Fj), -bj - cj Fj, 0}, {0, -bj - cj Fj,
b1 + c1 F1 + 2 (bj + cj Fj), -bj - cj Fj}, {-bj - cj Fj,
0, -bj - cj Fj, b2 + c2 F2 + 2 (bj + cj Fj)}
}


Step 2. Define the forward model functions

invCmodel = Inverse[cmodel];

ω1 = Sqrt[1/1.84 invCmodel[[1, 1]]];

ω1rec[c1_, b1_, c2_, b2_, cj_, bj_, F1_, F2_, Fj_] = ω1;

ω2 = Sqrt[1/1.84 invCmodel[[2, 2]]];

ω2rec[c1_, b1_, c2_, b2_, cj_, bj_, F1_, F2_, Fj_] = ω2;

J = 1/2 invCmodel[[1, 2]] Sqrt[ω1*ω2]/
Sqrt[invCmodel[[1, 1]] invCmodel[[2,2]]];

Jrec[c1_, b1_, c2_, b2_, cj_, bj_, F1_, F2_, Fj_] = J;


Note the use of Set (i.e, =) in defining the functions rather than the normal SetDelayed (i.e., :=).

Step 3. Define an objective function

We will take the measured data and compute reconstructed data. It is easy to make this more terse but a longer version may be easier to understand.

data is the set of measured data with F1 in column 1, ... ω1 in column 4 ...

The variable parameterVector represents the list {c1, b1, c2, b2, cj, bj}.

In the code funny lines that use the Sequence function enable us to use lists as arguments to the functions (e.g., ω1rec).

Reconstructed data is computed, then the residual is equated to the difference between the measured and reconstructed data. The objective function is the sum of the squares of the residuals.

objectiveFunction[data_, parameterVectors_] := Module[{
dataReconstructed,
residual
},

(* Apply the forward model equations to the
finger data to get reconstructed measurements *)

dataReconstructed = Map[Function[fingers,
{
ω1rec[Sequence @@ parameterVectors, Sequence @@ fingers],
ω2rec[Sequence @@ parameterVectors, Sequence @@ fingers],
Jrec[Sequence @@ parameterVectors, Sequence @@ fingers]
}],
data[[All, 1 ;; 3]]
];

(* Compute the residual *)

residual = dataReconstructed - data[[All, 4 ;; 6]];

(* The objective function is the sum of
the residuals squared *)

Total@Total@Map[#^2 &, residual, {2}]

]


Step 4. Apply FindMinimum

I did this twice (I'll only show you the second time). I discovered that the offset for Cj was less than zero.

FindMinimum can use constraints and works best if you supply it with reasonable initial values.

The solution, sol, represents the objective function and values for the parameters provided as rules.

sol = FindMinimum[
{
objectiveFunction[data, {c1, b1, c2, b2, cj, bj}],
c1 > 0 && c1 ∈ Reals && b1 > 0 && b1 ∈ Reals &&
c2 > 0 && c2 ∈ Reals && b2 > 0 && b2 ∈ Reals &&
cj > 0 && cj ∈ Reals && bj > -1*10^-4 && bj ∈ Reals
},
{{c1, 3.4*10^-4}, {b1, 1.6*10^-3}, {c2, 3*10^-4}, {b2,
3*10^-3}, {cj, 1.9*10^-4}, {bj, -4*10^-5}}
]


the result is

{0.0539887, {c1 -> 0.000340998, b1 -> 0.00157195, c2 -> 0.000301931,
b2 -> 0.00316187, cj -> 0.000192957, bj -> -0.0000434626}}


Step 5. Check the results

We will generate the reconstructed data to see how it compares to the measured data. I elected to use the With function which can be helpful in reducing the size of the code. You con't have to type in all the parameter values. You can use the output from FindMinimum (i.e., sol[[2]]) with ReplaceAll to minimize the typing.

With[
{
dataMeas = data[[All, 1 ;; 3]],
parameterVectors = {c1, b1, c2, b2, cj, bj} /. sol[[2]]
},
dataRec = Map[Function[fingers,
{
ω1rec[Sequence @@ parameterVectors, Sequence @@ fingers],
ω2rec[Sequence @@ parameterVectors, Sequence @@ fingers],
Jrec[Sequence @@ parameterVectors, Sequence @@ fingers]
}],
dataMeas
]
]


Here is a table of the reconstructed (left) and measured data (right).

Now let's plot it.

Column[{
Show[
ListPlot[
Transpose[{
data[[All, 4]],
dataRec[[All, 1]]
}],
PlotStyle -> {PointSize[Large], Red},
ImageSize -> 250,
PlotRange -> {{7, 9}, {7, 9}},
AxesLabel -> Style["Rec ω1", Bold, 12]
],
Plot[x, {x, 7, 9}, PlotStyle -> {Thin, Black}]
],
Spacer[{1, 10}],
Show[
ListPlot[
Transpose[{
data[[All, 5]],
dataRec[[All, 2]]
}],
PlotStyle -> {PointSize[Large], Red},
ImageSize -> 250,
PlotRange -> {{6.5, 8.5}, {6.5, 8.5}},
AxesLabel -> Style["Rec ω2", Bold, 12]
],
Plot[x, {x, 6.5, 8.5}, PlotStyle -> {Thin, Black}]
],
Spacer[{1, 10}],
Show[
ListPlot[
Transpose[{
data[[All, 6]],
dataRec[[All, 3]]
}],
PlotStyle -> {PointSize[Large], Red},
ImageSize -> 250,
PlotRange -> {{0.1, 0.6}, {0.1, 0.6}},
AxesLabel -> Style["Rec J", Bold, 12]
],
Plot[x, {x, 0.1, 0.6}, PlotStyle -> {Thin, Black}]
]
}]


The results look good:

Until new measurements are made we can model the capacitors as:

C1[F1_] := 0.000341*F1 + 0.001572
C2[F2_] := 0.0003019*F2 + 0.003162
Cj[Fj_] := 0.0001930*Fj - 0.00004236


I can unfortunately not just comment because I lack the reputation, which I would have found more appropriate. This case, to me, seems more like a question of actual fitting than the programing so I'll give you a few ideas that I had.

First of all, it could be useful to have analytic expressions for the invcmat elements. This has a simple reason. If those expressions are symmetric in C1, C2 your measured data is not giving you enough information to say which of your "identical" condensators has which value, you can just say that one of them has one value while the other one has the other - if the fits converge at all. This would mean that you need an additional measurement to tell them apart and then know which fitted value in each case belongs to which of those.

Second, my fitting routine would be to try a simple least square of differences fit simultaneously fitting \omega1, \omega2 and J with the three capacities as parameters (I know that professional software like origin can do this, I am unsure about mathematicas capacities in this point). This could be done for every measurement and then the capacities could be plotted over the number of fingers and you could then think about which function could be appropriate for the look of the data points (separately for each condensator). However, it should be noted that with 4/5 points a lot of functions will yield a graphically pleasing fit, so without any kind of theoretical model, this could end in a bit of guesswork.

I hope I did not misunderstand your question and thus oversimplify it. If I did so, I am sorry in advance for wasting your time.

• Thanks for the comment! I think you make some good suggestions. First of all, about the symmetry of the situation. I do in fact believe that it will end up being symmetric, but in fact I know which set of data belongs to which capacitor. One of the two will always have the lower value of $\omega$, so this is very clearly known in advance, even if the equations are symmetric. Is that what you meant? As for the fitting routine. So you'd suggest looking into simultaniously fitting the two omegas and J to the data, to find the C's, and then trying to link these C's to the number of fingers, right? – user129412 Sep 1 '15 at 19:18
• Yeah, this clears that point. I would then use the link provided by Daniel in the comments above to try the fitting of all datasets with mathematica, if this is what you are going to work with :) – Inari Sep 1 '15 at 19:21
• Yes, I did have a look at Daniels comment, it does seem a little complicated at this point so I'm trying to work it out. Thanks for the help! – user129412 Sep 1 '15 at 19:37