Finding fit parameters for a 2-dimensional matrix of data if the value of some parameters depend on the row and others on the column

Question

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!

Answer 1

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.