Re: [OctDev] regress and rcoplot (partially written, please help)

[Paul K. <pki...@us...>]

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