Here's what I wrote in wsolve:
[nr,nc] = size(A);
if nc > nr, error("underdetermined system"); end
## system solution: A x = y => x = inv(A) y
## QR decomposition has good numerical properties:
## AP = QR, with P'P = Q'Q = I, and R upper triangular
## so
## inv(A) y = P inv(R) inv(Q) y = P inv(R) Q' y = P (R \ (Q' y))
## Note that b is usually a vector and Q is matrix, so it will
## be faster to compute (y' Q)' than (Q' y).
[Q,R,p] = qr(A,0);
x = R\(y'*Q)';
x(p) = x;
s.R = R;
s.R(:,p) = R;
s.df = nr-nc;
s.normr = norm(y - A*x);
Here's what I wrote in polyconf:
## Confidence intervals for linear system are given by:
## x' p +/- sqrt( Finv(1-a,1,df) var(x' p) )
## where for confidence intervals,
## var(x' p) = sigma^2 (x' inv(A'A) x)
## and for prediction intervals,
## var(x' p) = sigma^2 (1 + x' inv(A'A) x)
##
## Rather than A'A we have R from the QR decomposition of A, but
## R'R equals A'A. Note that R is not upper triangular since we
## have already multiplied it by the permutation matrix, but it
## is invertible. Rather than forming the product R'R which is
## ill-conditioned, we can rewrite x' inv(A'A) x as the equivalent
## x' inv(R) inv(R') x = t t', for t = x' inv(R)
## Since x is a vector, t t' is the inner product sumsq(t).
## Note that LAPACK allows us to do this simultaneously for many
## different x using sqrt(sumsq(X/R,2)), with each x on a different row.
##
## Note: sqrt(F(1-a;1,df)) = T(1-a/2;df)
function [y,dy] = confidence(A,p,S,alpha,typestr)
if nargin < 4, alpha = []; end
if nargin < 5, typestr = 'ci'; end
y = A*p(:);
switch typestr,
case 'ci', pred = 0; default_alpha=erfc(1/sqrt(2));
case 'pi', pred = 1; default_alpha=0.05;
otherwise, error("use 'ci' or 'pi' for interval type");
end
if isempty(alpha), alpha = default_alpha; end
s = t_inv(1-alpha/2,S.df)*S.normr/sqrt(S.df);
dy = s*sqrt(pred+sumsq(A/S.R,2));
Can you use these two functions to do what you need?
- Paul
On Dec 31, 2005, at 3:06 AM, William Poetra Yoga Hadisoeseno wrote:
> In the Matlab documentation from http://www.mathworks.com (for regress),
> it mentions QR factorization (using QR factorization instead of
> computing an inverse matrix). Can anyone explain to me about this?
|