How can I calculate the confidence interval for "lsqnonneg" regression?

Hello All,

I did a "Isqnonneg" regression between vector(5,1) and matrix (5,3), I would like to find the confidence interval to my results, Please any help?

Thanks in advance,

Riyadh

[Riyadh Muttaleb]

Can I use "Lambda" that getting from "lsqnonneg" to find the confidence interval and how? and How can I evaluate the "lsqnonneg" result with lambda?

Thanks in advance

Answer 1

[John D'Errico]

The standard solutions for confidence intervals do not apply to bound constrained regression (lsqnonneg).

That normally leaves a bootstrap/jackknife solution, which can in theory compute confidence intervals on any such problem.

However, you have far too few data points (only 5, with 3 unknowns.) That suggests that any confidence intervals will be extremely wide, and VERY poorly estimable. For me to say this is a silly idea (i.e., a waste of time) might be taken wrongly. But IMHO, it is just that.

[Riyadh Muttaleb]

Even it is a silly idea, how can I do it?

Answer 2

[Walter Roberson]

 c = sym('c', [5 3], 'real');
  d = sym('d', [5 1], 'real');
  x = sym('x', [3 1], 'real');
  residue = sum((c*x - d).^2);
  solx1 = simplify( solve( diff(residue, x(1)), x(1)) );
  eq2 = simplify(subs(residue, x(1), solx1));
  solx2 = simplify( solve( diff(eq2, x(2)), x(2)) );
  eq3 = simplify(subs(eq2, x(2), solx2));
  solx3 = simplify( solve( diff(eq3, x(3)), x(3)) );
  X3 = solx3;
  X2 = (subs(solx2, x(3), solx3));
  X1 = ( subs( subs(solx1, x(2), X2), x(3), X3) );

Now X1, X2, and X3 are the general symbolic solution to the minimum residue for a 5 x 3 system, but with no range limitations. You can convert your actual values into rationals and substitute them in to get exact solutions. You can add a symbolic fudge factor to each of your actual values, and study the derivative of your X1, X2, X3 with respect to each of the fudge factors. For example,

C = sym(randi([-10 10], [5, 3]));        %actual values
D = sym(randi([0 20], [5, 1]));          %actual values
dC = sym('deltaC_', [5, 3]);             %per-variable fudge factor
dD = sym('deltaD_', [5, 1]);             %per-variable fudge factor

fudged_X3 = subs(subs(X3, c, C+dC), d, D+dD);

and now you can study fudged_X3 with respect to deltaC_4_2 for example.

In particular you would be interested in delta values that are +/- eps() of the actual value.

The expressions will be quite long. You might want to do statistical studies, substituting in a number of random values in the appropriate eps range and seeing how the X1, X2, X3 are affected.

Note: the above general solution could involve negative values. At the moment I am not sure of the best way to minimize the residue for that case.

[Riyadh Muttaleb]

Thank you for your answer,

I did lsqnonneg between these two matrices:

a=[1,0.9969,0.9947,0.9931,0.9918;10.9393,0.8999,0.8711,0.8485;1,0.5345,0.3484,0.2516,0.1935]
b=[3.4202,3.6494,3.8940,4.0478,4.14241]

and the result is

c=[3.8443,0,0]

I could not apply your code mathematically to these values to find the confidence interval,could you please help me,

I so appreciate your time and cooperation,

Thanks,

[Walter Roberson]

Your second row of a only has four columns. Should it be 1,0.9393 instead of 10.9393 ?

[Riyadh Muttaleb]

yes, exactly

[Walter Roberson]

Your data is 3 x 5 and 1 x 5, whereas your original question was about 5 x 3 and 5 x 1.

Using the notation I used above,

C = sym(a.');
D = sym(b.');
dC = sym('deltaC_', [5, 3]);             %per-variable fudge factor
dD = sym('deltaD_', [5, 1]);             %per-variable fudge factor

fudged_X1 = subs(subs(X1, c, C+dC), d, D+dD);
fudged_X2 = subs(subs(X2, c, C+dC), d, D+dD);
fudged_X3 = subs(subs(X3, c, C+dC), d, D+dD);

and now you study fudged_X3 with regards to the variables dC and dD .

fudge_vars = [dC(:); dD(:)] .';
fudge_vars_cell = num2cell(fudge_vars);
zfudge = 0 * fudge_vars;
exact_X1 = subs(fudged_X1, fudge_vars, zfudge);
exact_X2 = subs(fudged_X2, fudge_vars, zfudge);
exact_X3 = subs(fudged_X3, fudge_vars, zfudge);

eps_C = reshape( eps(a.'), 1, []);
eps_D = reshape( eps(b.'), 1, []);
eps_vars = [eps_C, eps_D];
nvars = length(eps_vars);

nprobe = 100;
R = cellfun(@sym, mat2cell( (rand(nprobe, nvars) * 2 - 1) .* repmat(eps_vars, nprobe, 1), nprobe, ones(1,nvars) ), 'uniform', 0);

X1_uncertainty = subs(fudged_X1 - exact_X1, fudge_vars_cell, R);
X2_uncertainty = subs(fudged_X2 - exact_X2, fudge_vars_cell, R);
X2_uncertainty = subs(fudged_X3 - exact_X3, fudge_vars_cell, R);

This may take some time (seriously; if you use like 1000 for nprobe then be prepared for a long wait.) It took about 3 1/2 minutes for 100 for me for X1_uncertainty alone.

You can now do things like

X1ud = double(X1_uncertainty);
mean(X1ud)
std(X1ud)

and otherwise calculate confidence intervals.

The sym(a.') and sym(b.') do not do exact binary conversion: by default sym() applied to a number approximates the number as a rational fraction (or rational multiple of Pi.) For example a(1,1) will be approximated as 9969/10000 but b(5) will be approximated as 2331969516551977/562949953421312 . If you know your inputs to higher precision than the floating point you posted, then use those higher precisions when constructing C and D.

The above uncertainty measurement still does not take into account the non-negative requirement. I have not worked out yet the effects that can have on the formula.

[Riyadh Muttaleb]

Thank you for your efforts, when I tried to run the code, there was an error "Undefined function or variable 'X1'"

[Walter Roberson]

You skipped the code I had posted earlier.

File attached.

[Riyadh Muttaleb]

Thank you so much, so after I get "mean" and "std" I have to apply these equations to get the range of confidence interval:

ts = tinv([0.025  0.975],length(x)-1);      
CI = mean + ts*std;

The result is mean(-7.5e-14) and std(9.8e-14) that is so far from my lsqnonneg result which is (3.8443,0,0)

Thank you for your help and time

[Walter Roberson]

Yes. The mean(X1ud) and std(X1ud) that I show are the mean and standard deviation of the uncertainty; so you do need to transform them to get confidence intervals. In particular you need to add back in the exact_X1 or exact_X2 or exact_X3 if you want them in terms of range.

You will find that two of the exact values are positive and one is negative; as I have been indicating, the equations here do not take the non-negative requirement into account, as I have not figured out the best way to use that information.

[Riyadh Muttaleb]

Thank you for your clarification, I have tried to get the confidence interval range by I couldn't, please could you help me to get confidence interval range for 95%. I did by using the above functions but the result of leqnonneg 3.8443 was out of the range.

Thank you for your great efforts,