octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Fitting with leasqr and uncertainties


From: CdeMills
Subject: Re: Fitting with leasqr and uncertainties
Date: Tue, 21 Jan 2014 15:53:23 -0800 (PST)

Mike Miller wrote
> Hi Pascal, I'm not familiar with this algorithm yet, but the name
> reminds me of some items that have come up before.
> 
> Can you comment on how this differs from the Matlab function prony in
> the signal processing toolbox [1]? There is no prony function in the
> signal package yet, but there should be eventually. Is there any
> overlap with what you have created here? Is this exactly what prony
> should be, or similar and can be adapted to the Matlab calling
> interface?
> 
> Also see bug #38063 [2] for some discussion on the expfit function in
> the optim package, where the Prony method was discussed also. Is there
> some overlap with that function as well or bits that can be borrowed?
> 
> [1] http://www.mathworks.com/help/signal/ref/prony.htmlThanks,
> [2] http://savannah.gnu.org/bugs/?38063

I'll try to be concise. The ML approach seems to be : find two vector, A and
B, such that
y(z) matches B(z) * u(z)/A(z); that is,  minimize the norm of (y(z) * A(z) -
B(z) u(z)). This can be solved as system of equation linear in the
coefficients of A and B:
[u  q u q^2 u q^m u  ... -q y -q^2 y ... -q^n y][b1 b2 ... bm a1 a2 ... an]
= y
where q is the unit shift operator. This is solved as an ordinary LS
problem, EXCEPT that the last n rows are arranged as a Toeplitz matrice AND
are noised too. Such approach is thus very sensitive about noise. R. Branham
showed that this noise is transformed into bias.

My approach use the same formulation, but the resolution is quite different
1) While constructing the LHS matrix, I add terms as [1 t t^2 ...] as
supplemental explanatory terms; this permits to model offset, time drift,
and more. 
2) The LHS matrix is solved in a mixed LS-TLS approach. I apply a QR
transform, and the last n rows and columns of R are solved using a TLS
algorithm; this is coherent with the noise structure. I stop at this stage,
as I'm only interested in the A vector.
3) The roots of this A vector are the system poles, that I directly
translate back into exponentials, sine and so on
4) The B part is solved by synthetising directly the exponentials and sine,
as well as adding terms as constant, time drift, and so on. This matrix is
returned by the pronysq as 'LB', and the B coefficients are computed by LS
(they do not contain noise in the ordinary sense)
5) I added a further refinement with a parameter called the 'Prediction
Horizon'; it is very effective at decorrelating noise.
6) Using error propagation technique, I compute estimates for the noise
covariance estimates of Z and B. The problem is non-linear in Z and linear
in B

To compare:
ML: only works on impulse response
pronysq: works on impulse as well as step-like signals. In this case,
transient and steady-state behaviour are returned separetelly

ML: returns B(z)/A(z)
pronysq: do the partial fraction decomposition and expansion of B(z)/A(z)

ML: solve the noised part of the LHS matrix by LS
pronysq: solve the noised part of the LHS using TLS

I also tried something similar to expfit; I gave up:
1) too sensitive to initial conditions
2) too sensitive to noise
3) tricky to achieve convergence; most of the times wandering around the
expected solution. Error decay stops in the vicinity of the nominal value

About overlap: 
1) pronysq is more generic than ML; their output are not compatibles
(explicit partial fraction expansion)
2) pronysq is not iterative and a lot of efforts have been put to reduce the
noise sensitivity. 

Regards

Pascal



--
View this message in context: 
http://octave.1599824.n4.nabble.com/Re-Fitting-with-leasqr-and-uncertainties-tp4661222p4661230.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]