1
$\begingroup$

I have some 2-dimensional data of dimension n x m that I suspect has the following structure:

syntheticData[i_, j_] := (a[[i]] + b[[i]] t[[j]] + c[[i]] t[[j]]^2 ) (1 + 
RandomReal[{-0.2, 0.2}])

where a, b and c are n-dimensional vectors and t is an m-dimensional vector.

I am trying to use NonLinearModelFit to obtain all the components of a, b and c (it would be nice to get t too, but I could find its value later if I knew a, b and c), but when I use t as the independent variable for the fit I do not know how to tell Mathematica that the value of t must be the same for each row.

For completeness, the code to simulate data similar to mine would be:

a = Table[RandomReal[{0, 5}], 4];
b = Table[RandomReal[{10, 30}], 4];
c = Table[RandomReal[{0, 10}], 4];
t = Table[RandomReal[{0.5, 1.5}], 50];

fullTable = Table[syntheticData[i, j], {i, 4}, {j, 50}]

and I am interested in finding a, b and c (with the associated uncertainty, if possible)

Edit:

My data has some missing data, further complicating the problem. It looks something like this (although each row has at least two non-zero entries)

missingTable = Table[If[RandomReal[] < 0.3, "", fullTable[[i, j]]], {i, 4}, {j, 50}]

My approach right now is to write dummy variables as

{a, b, c} =  Transpose @ Table[ToExpression[# <> ToString@index] & /@ {"a", "b", 
 "c"}, {index, 4}]; 
t = Table[ToExpression["t" <> ToString@index], {index, 50}];

and then minimize the following expression with respect to them:

NMinimize[
         Sum[
            (Sum[ If[NumberQ @ fulltable[[n, i]], 
                    (fulltable[[n, i]] - a[[i]] + b[[i]] t[[j]] + c[[i]] t[[j]]^2 )^2,
     0], {i, 4}])^2, 
                {j, 50}], 
         Join[a, b, c, t]]

The choice of what function to minimize is somewhat arbitrary, but the reasoning behind it is that as I assume the t to be the same for all the points in a given row, I square the sum of the residuals of these components to prioritize this fit.

Another possibility would be to define the residuals as

fulltable[[n, i]] / (a[[i]] + b[[i]] t[[j]] + c[[i]] t[[j]]^2) - 1

which could account for the errors being proportional to the magnitud and not linear.

If anyone has a better idea I would love to try it!

$\endgroup$
2
  • 1
    $\begingroup$ You don't need to generate variable names in this way. Consider this toy example Minimize[Sum[(s[i] - 1)^2, {i, 1, 3}], Array[s, 3]] $\endgroup$
    – mikado
    May 22, 2022 at 10:02
  • $\begingroup$ Thank you, that is a much better system than what I was using! I had tried using Part in the past and I would get errors, but using the argument for the indexes seems like the way to go $\endgroup$
    – edmateosg
    May 22, 2022 at 15:06

1 Answer 1

1
$\begingroup$

I've had a look at your problem, and I'm not going to try to offer a full solution.

The general approach I would take to this type of problem is "Alternating Least Squares". That is, I would attempt to solve alternately for {a,b,c} and t, at each step using my best guess for the the other variable set.

For example, given an estimate for t the following function will estimate {a,b,c}

estimateabc1 =.
estimateabc1[t_List, table_?VectorQ] := 
 LinearModelFit[Transpose[{t, table}], {t1, t1^2}, t1]

estimateabc2 =.
estimateabc2[t_List, table_?MatrixQ] := estimateabc1[t, #] & /@ table

estimateabc =.
estimateabc[t_List, table_?MatrixQ] := 
 Transpose[#["BestFitParameters"] & /@ estimateabc2[t, table]]

The {a,b,c} estimates can be updated with

{aest, best, cest} = estimateabc[test, fullTable]

However, looking at this step, it appears that your data doesn't allow you to get good estimates of {a,b,c}. Bearing this in mind, there are a number of steps I would take.

  1. Check whether your measurement errors are really as big as this model suggests.
  2. Perhaps try fitting a linear model rather than a quadratic model. If nothing else, a linear model might give you a decent starting point for fitting a quadratic.
  3. Note that nothing in your problem allows you to fix the time origin. I would approach this by fixing the mean of t as zero (following every estimation step).
  4. Of course, to apply alternate least squares you need an intial estimate for either t or {a,b,c}. In many problems, random initialisation will work. If different random errors give you different solutions, take the one with the smallest residual errors.

I hope this helps.

EDIT

In response to the updates:

  1. I would try alternating least squares with random (or manual initiation) before trying anything complicated for initiation.
  2. With an iterative approach, I would handle missing values by setting them to their currently estimated value. In this way, they should have no effect on the estimates.
  3. If your data is positive and your errors are multiplicative, I would consider taking the log of the data so you have additive errors. Of course, this changes the meaning of terms in your model too.
$\endgroup$
1
  • $\begingroup$ What I have been trying to do to fit them simultaneously is to define residuals and use NMinimize to fit them simultaneously, and I had thought of then using these results as the starting point for the alternating least squares, as you suggest. A second problem that I have is that some entries of my data table are missing, further complicating the problem. I will edit my question to better reflect this $\endgroup$
    – edmateosg
    May 21, 2022 at 16:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.