4
$\begingroup$

I have a problem creating MeanPredictionBands for a FittedModel from the ResourceFunction["MultiNonlinearModelFit"] if the FittedModel includes ParametricFunctions. A simple example:

data1 = {#, 0.05 E^(0.075 #)} & /@ Range[10];
data2 = {#, 0.075 E^(0.075 #)} & /@ Range[10];

model = ParametricNDSolveValue[{x'[t] == u x[t], x[0] == x0}, x, {t, 0, 100}, {u, x0}];

fit = ResourceFunction["MultiNonlinearModelFit"][{data1, 
   data2}, {model[u, x0][t], model[u, x02][t]}, {u, x0, x02}, {t}];

fit["ParameterTable"];
fit["MeanPredictionBands"]

Experimental`NumericalFunction::nnum: The function value Switch[Round[\[FormalN]],1,<<1>>,2,ParametricFunction[1,Internal`Bag[<1>],0,1,False,{{u$214682,x0$214683},System`Utilities`HashTable[<3>],{},{},{1,2},{Automatic,0,0},{0,1}},{NDSolve`base$214690,NDSolve`NDSolveParametricFunction[0,{ParametricNDSolveValue,Internal`Bag[<2>],None,ParametricNDSolveValue},<<7>>,{},<<1>>]}][u,x02][t]] is not a number at {u,x0,x02} = {0.0749998,0.0499999,0.0749999}.

The MultiNonlinearModelFit works and gives a good solution in the ParameterTable, but the fit[MeanPredictionBands] gives an error. It seems that the issue is with the ParametricFunction in the FittedModel. I have been able to get the MeanPredictionBands for a FittedModel from the MultiNonlinearModelFit if it does not contain ParametricFunctions. How can I get the mean prediction bands in this case?

$\endgroup$
8
  • $\begingroup$ Please add the text of the error you get by editing your question. $\endgroup$
    – MarcoB
    Sep 12, 2022 at 11:29
  • 2
    $\begingroup$ Note that MultiNonlinearModelFit only allows (at least the last time I looked) a common error variance. While estimates of parameters might not be affected so much, if the error variance of the different models differ, then the prediction bands can be greater affected. Your example above has no error term for either model which doesn't make for a very realistic proxy for your real data. $\endgroup$
    – JimB
    Sep 12, 2022 at 18:26
  • $\begingroup$ Thank you JimB. My real data have usually multiple measurements for each point and I weigh the data with the variance. Mostly my models are very similar and typically only vary in initial conditions, and I expect similar error variance in the models. $\endgroup$
    – JackySnoep
    Sep 13, 2022 at 7:10
  • $\begingroup$ @JimB MultiNonlinearModelFit allows for different weights to be applied to each individual data point across the different datasets. You need to specify them as Weights -> {{w11, w12, ...}, {w21, w22, ...}, ...}. $\endgroup$ Oct 10, 2022 at 7:22
  • 1
    $\begingroup$ @SjoerdSmit (Preface: As stated below MultiNonlinearModelFit is great.) Using the Weight option will allow for independent errors with different variances but only if one "knows" the ratios of the variances. In biological applications it's almost certainly wishful thinking that one knows the ratios and that the errors would be independent. So it's about the "ability to assume" rather than the "willingness to assume". I would hope that someday NonlinearModelFit and LinearModelFit would allow for mixed models (i.e., models with more than one error term and correlated errors). $\endgroup$
    – JimB
    Oct 10, 2022 at 13:31

1 Answer 1

4
$\begingroup$

It appears as though NonlinearModelFit/FittedModel are unable to compute the required gradients themselves, so you'll have to update the definition of MultiNonlinearModelFit:

Options[MultiNonlinearModelFit] = {AccuracyGoal->Automatic,ConfidenceLevel->19/20,EvaluationMonitor->None,Gradient->Automatic,MaxIterations->Automatic,
    Method->Automatic,PrecisionGoal->Automatic,StepMonitor->None,Tolerance->Automatic,VarianceEstimatorFunction->Automatic,
    Weights->Automatic, WorkingPrecision->Automatic, "DatasetIndexSymbol"->\[FormalN]
}; (* \[Equal] Join[Options[NonlinearModelFit], {"DatasetIndexSymbol" -> \[FormalN]}] *)

MultiNonlinearModelFit[datasets_, form_, fitParams_, independents : Except[_List], opts : OptionsPattern[]] := 
    MultiNonlinearModelFit[datasets, form, fitParams, {independents}, opts];
 
MultiNonlinearModelFit[datasets_, form : Except[_?AssociationQ], fitParams_, independents_, opts : OptionsPattern[]] := 
    MultiNonlinearModelFit[datasets, <|"Expressions" -> form, "Constraints" -> True|>, fitParams, independents, opts];
 
MultiNonlinearModelFit[
    datasets : {__?(MatrixQ[#1, NumericQ] &)}, 
    assoc : KeyValuePattern[{
        "Expressions" -> expressions_,
        "Constraints" -> constraints_
    }] /; AssociationQ[assoc],
    fitParams_List, 
    independents_List,
    opts : OptionsPattern[]
] := Module[{
    fitfun, weights,
    numSets = Length[datasets],
    precision = Precision @ datasets,
    augmentedData,
    indexSymbol = OptionValue["DatasetIndexSymbol"],
    grad
},
    augmentedData = Join @@ MapIndexed[
        Join[ConstantArray[N[#2, precision], Length[#1]], #1, 2]&,
        datasets
    ];
    fitfun = With[{
        conditions = Flatten @ Map[
            {#, Indexed[expressions, #]}&, 
            Range[numSets]
        ]
    }, 
        Switch @@ Prepend[conditions, Round[indexSymbol]]
    ];
    grad=D[fitfun, {Replace[fitParams, {v_, ___} :> v , 1]}];
    weights = Replace[
        OptionValue[Weights],
        {
            (list_List)?(VectorQ[#1, NumericQ]& ) /; Length[list] === numSets :> 
                Join @@ MapThread[ConstantArray, {list, Length /@ datasets}], 
            list : {__?(VectorQ[#1, NumericQ] & )} /; Length /@ list === Length /@ datasets :>
                Join @@ list, 
            "InverseLengthWeights" :> Join @@ Map[
                ConstantArray[N[1 / #1, precision], #1]&,
                Length /@ datasets
            ]
        }
    ]; 
    NonlinearModelFit[
        augmentedData,
        If[TrueQ[constraints], fitfun, {fitfun, constraints}], 
        fitParams,
        Flatten[{indexSymbol, independents}],
        Weights -> weights, 
        Sequence @@ FilterRules[{opts}, Options[NonlinearModelFit]],
        Gradient->grad
    ]
];

This is copied almost straight from the definition notebook of MultiNonlinearModelFit (available via ResourceFunction["MultiNonlinearModelFit", "DefinitionNotebook"], except that I manually supply the gradient of the function:

I compute it as (thanks @SjoerdSmit for the significant simplification!)

grad=D[fitfun, {Replace[fitParams, {v_, ___} :> v , 1]}];

This is effectively computing the gradients of each fit function separately and combining them using the same Switch[Round[...], ...] strategy that MultiNonlinearModelFit is using for the fitting itself.

The gradient is then supplied to NonlinearModelFit in the last expression, which means we can now compute the prediction bands:

fit["MeanPredictionBands"]

enter image description here

We can now insert a value for the formal index to select the bands for one of the functions:

data1 = {#, 0.05 E^(0.075 #) + 0.005 RandomReal[]} & /@ Range[10];
data2 = {#, 0.075 E^(0.075 #)} & /@ Range[10];

model = ParametricNDSolveValue[{x'[t] == u x[t], x[0] == x0}, 
   x, {t, 0, 100}, {u, x0}];

fit = MultiNonlinearModelFit[{data1, data2}, {model[u, x0][t], 
    model[u, x02][t]}, {u, x0, x02}, {t}];

Show[Plot[
  fit["MeanPredictionBands"] /. \[FormalN] -> 1 // Evaluate, {t, 0, 
   10}, PlotStyle -> Red, Filling -> 1 -> {2}], ListPlot[data1]]

enter image description here

(Note that I have added some noise to data1 to make the bands visible)

$\endgroup$
4
  • 1
    $\begingroup$ Thank you Lukas, I need to study the definition a bit more but your solution looks good. I find the MultiNonlinearModelFit function very useful, and with the explicit "grad" definition you used, I hope to be able to use the "MeanPredictionBands" for a wider set of FittedModels. $\endgroup$
    – JackySnoep
    Sep 12, 2022 at 14:46
  • 1
    $\begingroup$ +1 for the modification of MultiNonlinearModelFit (and despite my multiple complaints about MultiNonlinearModelFit, it is wonderfully constructed and useful function). However, that single complaint I have ("assuming a common error variance") results in this case with the resulting prediction interval being too small for data1 and too large for data2. $\endgroup$
    – JimB
    Sep 12, 2022 at 18:34
  • $\begingroup$ Thanks for figuring this out! Is there a reason for not just computing the gradient with grad = D[fitfun, {Replace[fitParams, {v_, ___} :> v , {1}]}]? It seems to give the exact same results for me that way. $\endgroup$ Oct 10, 2022 at 11:56
  • $\begingroup$ @SjoerdSmit Not really - I could have sworn that I first tried something like that, but evidently, I failed. I can't see a reason why it should ever fail. Thanks! $\endgroup$
    – Lukas Lang
    Oct 10, 2022 at 17:08

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.