3
$\begingroup$

I have vector-valued data with normally distributed noise that I want to analyse. Assumptions on data can be made; in this specific example, it's expected that components of fitted vector values always have a unit sum.

Here's a toy example of manually finding a fit with these constraints:

With[{source = {x^2 / 2, 1 - x^2 / 2 - x, x}},
 With[{
   coeff = {a, b, c, d, e, f},
   eqn =
    {a x^2 + b x + c, d x^2 + e x + f, 
     1 - (a x^2 + b x + c) - (d x^2 + e x + f)},
   data = 
    Table[{x, 
      source + RandomVariate[NormalDistribution[0, .1], 3]}, {x, 0, 1, .025}]},
  Sum[SquaredEuclideanDistance[v[[2]], eqn /. x -> v[[1]]], {v, data}] //
     eqn /. Last@Minimize[# /. Abs -> Identity, coeff, Reals] & // Echo //
   Show[
    ListPlot[
     Transpose[{{#1, #2[[1]]}, {#1, #2[[2]]}, {#1, #2[[3]]}} & @@@ data]],
    Plot[Join[source, #], {x, 0, 1}, Evaluated -> True]] &]]

{-0.0162628+0.0842304 x+0.43683 x^2,0.99686 -0.908144 x-0.608563 x^2,0.0194033 +0.823914 x+0.171733 x^2}

enter image description here

Now I have a curve fit for each vector component, but what I would really want is corresponding FittedModel objects with "MeanPredictionBands" and such properties. Is there any way to cross this gap, practically creating a model with the above fitting polynomial as "BestFit"?

EDIT: Looking at InputForm of FittedModel objects would suggest this is feasible. Most of the arguments seem straight-forward, although some on the end of the definition are not that obvious.

Improve this question
$\endgroup$
3
  • $\begingroup$ Possible duplicate: mathematica.stackexchange.com/questions/87927/…. $\endgroup$
    – JimB
    Mar 9, 2019 at 15:15
  • 2
    $\begingroup$ The potential downside to the possible duplicate is that one would need to be able to assume that the error variances are the same for all 3 equations. If that's not a feasible assumption for real-world data, then using LogLikelihood and FindMaximum functions would allow for potentially unequal variances. $\endgroup$
    – JimB
    Mar 9, 2019 at 15:21
  • 1
    $\begingroup$ The answer I gave below is wrong if the triplets of observations sum exactly to 1 as opposed to just the expected values summing to one. If the former, then by definition the 3 observations are correlated with each other and that doesn't yet seem to be addressed by either of the two current answers. (The estimates of the coefficients won't change much - if at all - but the estimates of the standard errors of the coefficients and the predictions will be off. $\endgroup$
    – JimB
    Mar 11, 2019 at 1:44

2 Answers 2

Reset to default
3
$\begingroup$

While you know how to put this altogether from scratch, I suspect that some will not be able to construct the confidence or prediction bands reliably in such a manner. If it is a reasonable assumption that the error variances are the same for all models, then just modifying the input data is about all that's necessary to get the desired output from NonlinearModelFit. (I make too many mistakes so I need to go that route.)

First, include in some dummy variables in the data:

(* Sample data using dummy variables:
1,0,0 => data from model 1
0,1,0 => data from model 2
0,0,1 => data from model 3 *)
data={{1, 0, 0, 0., -0.0988754}, {0, 1, 0, 0., 0.936639}, {0, 0, 1, 0., -0.168444}, {1, 0, 0, 0.025, -0.0201999}, {0, 1, 0, 0.025, 1.02862}, {0, 0, 1, 0.025, -0.0562009}, {1, 0, 0, 0.05, -0.0120367}, {0, 1, 0, 0.05, 0.870238}, {0, 0, 1, 0.05, -0.0936423}, {1, 0, 0, 0.075, -0.128804}, {0, 1, 0, 0.075, 0.985538}, {0, 0, 1, 0.075, 0.125768}, {1, 0, 0, 0.1, 0.078006}, {0, 1, 0, 0.1, 0.775311}, {0, 0, 1, 0.1, 0.0928273}, {1, 0, 0, 0.125, 0.121525}, {0, 1, 0, 0.125, 0.71475}, {0, 0, 1, 0.125, 0.0477607}, {1, 0, 0, 0.15, 0.115014}, {0, 1, 0, 0.15, 0.851672}, {0, 0, 1, 0.15, 0.0933949}, {1, 0, 0, 0.175, -0.0260633}, {0, 1, 0, 0.175, 0.721085}, {0, 0, 1, 0.175, 0.2533}, {1, 0, 0, 0.2, 0.0929319}, {0, 1, 0, 0.2, 0.741054}, {0, 0, 1, 0.2, 0.246308}, {1, 0, 0, 0.225, -0.121485}, {0, 1, 0, 0.225, 0.558944}, {0, 0, 1, 0.225, 0.180121}, {1, 0, 0, 0.25, 0.196925}, {0, 1, 0, 0.25, 0.739204}, {0, 0, 1, 0.25, 0.216802}, {1, 0, 0, 0.275, 0.00288085}, {0, 1, 0, 0.275, 0.531013}, {0, 0, 1, 0.275, 0.364702}, {1, 0, 0, 0.3, -0.0748257}, {0, 1, 0, 0.3, 0.452435}, {0, 0, 1, 0.3, 0.246644}, {1, 0, 0, 0.325, -0.102033}, {0, 1, 0, 0.325, 0.552273}, {0, 0, 1, 0.325, 0.307822}, {1, 0, 0, 0.35, 0.012419}, {0, 1, 0, 0.35, 0.493619}, {0, 0, 1, 0.35, 0.109975}, {1, 0, 0, 0.375, 0.174347}, {0, 1, 0, 0.375, 0.403849}, {0, 0, 1, 0.375, 0.640423}, {1, 0, 0, 0.4, 0.0267628}, {0, 1, 0, 0.4, 0.50906}, {0, 0, 1, 0.4, 0.249802}, {1, 0, 0, 0.425, 0.145036}, {0, 1, 0, 0.425, 0.401913}, {0, 0, 1, 0.425, 0.365026}, {1, 0, 0, 0.45, 0.0917381}, {0, 1, 0, 0.45, 0.452611}, {0, 0, 1, 0.45, 0.191856}, {1, 0, 0, 0.475, 0.05115}, {0, 1, 0, 0.475, 0.461663}, {0, 0, 1, 0.475, 0.370011}, {1, 0, 0, 0.5, 0.275475}, {0, 1, 0, 0.5, 0.20153}, {0, 0, 1, 0.5, 0.479379}, {1, 0, 0, 0.525, 0.184857}, {0, 1, 0, 0.525, 0.26697}, {0, 0, 1, 0.525, 0.586871}, {1, 0, 0, 0.55, 0.236133}, {0, 1, 0, 0.55, 0.379411}, {0, 0, 1, 0.55, 0.61064}, {1, 0, 0, 0.575, 0.241561}, {0, 1, 0, 0.575, 0.243509}, {0, 0, 1, 0.575, 0.508213}, {1, 0, 0, 0.6, 0.220454}, {0, 1, 0, 0.6, 0.118133}, {0, 0, 1, 0.6, 0.595092}, {1, 0, 0, 0.625, 0.12719}, {0, 1, 0, 0.625, 0.126479}, {0, 0, 1, 0.625, 0.644849}, {1, 0, 0, 0.65, 0.114451}, {0, 1, 0, 0.65, 0.144761}, {0, 0, 1, 0.65, 0.581061}, {1, 0, 0, 0.675, 0.197143}, {0, 1, 0, 0.675, 0.186537}, {0, 0, 1, 0.675, 0.705692}, {1, 0, 0, 0.7, 0.231923}, {0, 1, 0, 0.7, -0.0542245}, {0, 0, 1, 0.7, 0.939961}, {1, 0, 0, 0.725, 0.318287}, {0, 1, 0, 0.725, 0.0657639}, {0, 0, 1, 0.725, 0.609313}, {1, 0, 0, 0.75, 0.177536}, {0, 1, 0, 0.75, 0.0269016}, {0, 0, 1, 0.75, 0.741112}, {1, 0, 0, 0.775, 0.357059}, {0, 1, 0, 0.775, -0.191036}, {0, 0, 1, 0.775, 0.830267}, {1, 0, 0, 0.8, 0.188352}, {0, 1, 0, 0.8, -0.0579005}, {0, 0, 1, 0.8, 0.695232}, {1, 0, 0, 0.825, 0.271996}, {0, 1, 0, 0.825, -0.304167}, {0, 0, 1, 0.825, 0.832492}, {1, 0, 0, 0.85, 0.408437}, {0, 1, 0, 0.85, -0.0930686}, {0, 0, 1, 0.85, 0.717563}, {1, 0, 0, 0.875, 0.570663}, {0, 1, 0, 0.875, -0.286103}, {0, 0, 1, 0.875, 0.96456}, {1, 0, 0, 0.9, 0.340535}, {0, 1, 0, 0.9, -0.25579}, {0, 0, 1, 0.9, 0.672478}, {1, 0, 0, 0.925, 0.0811443}, {0, 1, 0, 0.925, -0.332707}, {0, 0, 1, 0.925, 0.841893}, {1, 0, 0, 0.95, 0.327199}, {0, 1, 0, 0.95, -0.43617}, {0, 0, 1, 0.95, 0.849182}, {1, 0, 0, 0.975, 0.359762}, {0, 1, 0, 0.975, -0.67714}, {0, 0, 1, 0.975, 0.950295}, {1, 0, 0, 1., 0.494625}, {0, 1, 0, 1., -0.533217}, {0, 0, 1, 1., 1.09384}};

(* Define equation using dummy variables *)
eqn[m1_, m2_, m3_, x_, a_, b_, c_, d_, e_, f_] :=
 {m1, m2, m3}.{a x^2 + b x + c, d x^2 + e x + f, 
   1 - (a x^2 + b x + c) - (d x^2 + e x + f)}

(* Run regression *)
coeff = {a, b, c, d, e, f};
nlm = NonlinearModelFit[data, eqn[m1, m2, m3, x, a, b, c, d, e, f], 
   coeff, {m1, m2, m3, x}];

(* Get mean prediction bands *)
mpb1 = nlm["MeanPredictionBands"] /. {m1 -> 1, m2 -> 0, m3 -> 0};
mpb2 = nlm["MeanPredictionBands"] /. {m1 -> 0, m2 -> 1, m3 -> 0};
mpb3 = nlm["MeanPredictionBands"] /. {m1 -> 0, m2 -> 0, m3 -> 1};

(* Show data, fit, and mean prediction bands *)
Show[ListPlot[data[[All, {4, 5}]]],
 Plot[{nlm[1, 0, 0, x], nlm[0, 1, 0, x], nlm[0, 0, 1, x]}, {x, 0, 1},
  PlotStyle -> {Blue, Red, Green}],
 Plot[{mpb1, mpb2, mpb3}, {x, 0, 1},
  PlotStyle -> {Blue, Blue, Red, Red, Green, Green}]]

Data, fit, and mean prediction bands

And this works even if the model is not a polynomial.

Improve this answer
$\endgroup$
1
$\begingroup$

I dug a bit into FittedModel and it would seem pretty straight-forward to reconstruct FittedModels for simple linear model fits of one-dimensional data. As a result, I created a PreFittedLinearModel function which takes data, fitted polynomial (plus options, such as Weights) as input, and returns a FittedModel object as a result:

ClearAll@PreFittedLinearModel;

Options[PreFittedLinearModel] = 
  FilterRules[Options[LinearModelFit], 
   ConfidenceLevel | Tolerance | VarianceEstimatorFunction | Weights |
     WorkingPrecision];

PreFittedLinearModel[data_, poly_, var_, OptionsPattern[]] :=
 With[{funs = Table[var^i, {i, 0, Exponent[poly, var]}]},
  FittedModel @@ {
    {"Linear", N@CoefficientList[poly, var], {{var}, funs}, {0, 0}},
    {N[OptionValue[Weights] /. Automatic -> Table[1, Length@data]]}, 
    data, Table[N@funs, {var, data[[All, 1]]}], 
    Function[Null, Internal`LocalizedBlock[{var}, #], {HoldAll}], 
    Sequence @@ (# -> OptionValue[#] & /@ {ConfidenceLevel, Tolerance,
         VarianceEstimatorFunction, WorkingPrecision})}]

Now I can check how (95%) confidence intervals compare with my original function:

With[{source = {x^2/2, 1 - x^2/2 - x, x}},
 With[{
    coeff = {a, b, c, d, e, f},
    eqn = {
      a x^2 + b x + c, d x^2 + e x + f, 
      1 - (a x^2 + b x + c) - (d x^2 + e x + f)},
    data = 
     Table[{x, 
       source + RandomVariate[NormalDistribution[0, .1], 3]}, {x, 0, 1, .025}]},
   Sum[SquaredEuclideanDistance[v[[2]], eqn /. x -> v[[1]]], {v, data}] //
     eqn /. Last@Minimize[# /. Abs -> Identity, coeff, Reals] & //
     PreFittedLinearModel[{#1, #2[[3]]} & @@@ data, #[[3]], x] &] //
   Plot[Join[{x}, #["MeanPredictionBands"]], {x, 0, 1}, 
    Evaluated -> True] &]

enter image description here

... but, as @JimB stated, this approach has issues in its core.

EDIT: Probably I should modify this to work with any linear combination of functions, and extract polynomial coefficient functions for this generalized version instead of assuming a polynomial. Polynomials just happened to be the tool I used originally...

Improve this answer
$\endgroup$

Your Answer

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

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