MeanPredictionBands not working in ResourceFunction["MultiNonlinearModelFit"]

Question

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?

Answer 1

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"]

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]]

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