[Octave-cvsupdate] SF.net SVN: octave:[10445] trunk/octave-forge/extra/control-devel/inst/arx. m

[<par...@us...>]

Revision: 10445
          http://octave.svn.sourceforge.net/octave/?rev=10445&view=rev
Author:   paramaniac
Date:     2012-05-15 15:22:00 +0000 (Tue, 15 May 2012)
Log Message:
-----------
control-devel: arx: use Ljung's QR method

Modified Paths:
--------------
    trunk/octave-forge/extra/control-devel/inst/arx.m

Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m	2012-05-15 14:40:12 UTC (rev 10444)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m	2012-05-15 15:22:00 UTC (rev 10445)
@@ -95,14 +95,17 @@
 
 
 function theta = __theta__ (phi, y, i, n)
-                                                
-  ## Theta = Phi \ Y(n+1:end, :);                   # naive formula
     
-  if (numel (phi) == 1)
-    ## single-experiment dataset
-    theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i));
-  else
-    ## multi-experiment dataset
+  if (numel (phi) == 1)                             # single-experiment dataset
+    A = horzcat (phi{1}, y{1}(n(i)+1:end, i));      # [Phi, Y]
+    R0 = triu (qr (A, 0));                          # 0 for economy-size R (without zero rows)
+    R1 = R0(1:end-1, 1:end-1);                      # R1 is triangular - can we exploit this in R1\R2
+    R2 = R0(1:end-1, end);
+    theta = __ls_svd__ (R1, R2);                    # R1 \ R2
+    
+    ## Theta = Phi \ Y(n+1:end, :);                 # naive formula
+    ## theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i));
+  else                                              # multi-experiment dataset
     ## TODO: find more sophisticated formula than
     ## Theta = (Phi1' Phi + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...)
     

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.