MATLAB: Functional derivative of anonymous function of parameters with known errors

functionsstatistical toleranciingstatistical tolerancing

Hi all, I have a fairly complicated function of time and several parameters:
f(t;p_1,…,p_k),
each of which has a known error {dp_1,…,dp_k} (which give 95% confidence intervals). I would like to plot this function along with confidence bands derived from these known bands.
Currently, I have my function written like this:
p_1 = 1;
p_2 = 2;
...
p_k = k;
dp_1 = .1;
...
f = @(t) somefunction(t,p_1,...,p_k);
Ideally, I would be able to generate two functions that I could then plot as bands around f, of the following form (D denotes the partial derivative)
f+ = Df/Dp_1 * dp_1 + ... + Df/Dp_k * dp_k;
f- = Df/Dp_1 * (-dp_1) + ... + Df/Dp_k * (-dp_k);
to approximate the functional derivative. Does anybody have suggestions on the best way to calculate these functions?
*******************************************
Alternatively, I know that it's easy to calculate the total derivative of f w/ respect to t (i.e., the non-parameter variable) — does anybody know of an easy way to establish relatively accurate upper and lower confidence bands on f(t) directly from f'(t) given confidence intervals for its parameters? It seems to me that you can't get around needing some kind of knowledge about partial derivatives, but I'm adding this addendum to the question in case I'm wrong about that.

Best Answer

This is the domain of a branch of statistics known as statistical tolerancing. Thus you have a function of some parameters that have some assumed distributions, then you want to compute the distribution of that function. That the "function" is an ODE is irrelevant. At any point in time, you can theoretically compute the distribution of the function at that point.
However, simple use of the first derivative is often woefully inadequate.
For example, given that x is normally distributed, with mean 0 and variance 1, and a very simple function f(x)=x^2, would the first derivative of f be sufficient to put useful bands around the result? Somehow, I doubt it.
There are many common schemes to solve the general tolerancing problem.
1. Monte Carlo simulation - thus generate random samples for the parameters, then compute the "function" for each. Confidence intervals will be based on percentiles at each time step.
2. First or second order Taylor series approximations. What your question was aimed at originally. This may take more work than you need to do, compared to the alternatives.
3. Modified Taguchi methods. These reduce the problem to an implicit numerical integration. But they make the problem trivial to solve. See this reference for a paper I wrote with Nick Zaino.
https://www.jstor.org/stable/1269802?seq=1#page_scan_tab_contents
Do I need to provide some examples? sigh. let me see...
Trivial is Monte Carlo:
%%Assume that u is normally distributed, with mean 1, standard deviation 0.1.
% Solve the problem 
%


%  y' = u*y
%
% over the domain [0,1], subject to the initial condition y(0), which
% is also assumed to be normally distributed, with mean of 0.5 &
% standard deviation of 0.12.
%
% So u and y(0) both have SOME value, but we don't know what is the value.
% Only that they are normally distributed with the given mean and variance.
% in fact, other distributions than normal can be allowed, but the solution
% must deal properly with the distribution.
%%Monte Carlo...
nsim = 1000;
y0 = randn(1,nsim)*0.12 + 0.5;
u = randn(1,nsim)*0.1 + 1;
nt = 20;
ysim = zeros(nt,nsim);
yp = @(t,y,U) U*y;
tspan = linspace(0,1,nt)';
for i = 1:nsim
  [tpred,ysim(:,i)] = ode45(@(t,y) yp(t,y,u(i)),tspan,y0(i));
end
% compute 5% and 95% percentiles of ysim
plo = 0.05;
phi = 0.95;
ysim = sort(ysim,2);
ylo = ysim(:,floor(plo*nsim));
yhi = ysim(:,floor(phi*nsim));
plot(tspan,[ylo,yhi],'-')
grid on
So the green line is the 95'th percentile, blue is the 5% line.
Remember that the point at t=0 is NOT an independent sample along the curve from the point at t=0.1, etc.
Other methods will produce similar predictions. The difference is that the modified Taguchi methods require solving the differential equation system only a relatively few times, even if a higher order method is employed.
Related Question