🚀 go-pugleaf

In article <1456@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
>rob@kaa.eng.ohio-state.edu (Rob Carriere) writes:
><In article <2074@cs-spool.calgary.UUCP> smit@enel.ucalgary.ca (Theo Smit)
><writes: 
><[stability force deleted]
><>Mostly this is used to counter arithmetic errors; if your system is unstable
><>the poles will likely be well outside the unit circle and you have no hope
   in 
><>heck of doing much about it.
><
><If your system is unstable and you force the estimate to be stable you have a
><meaningless estimate.
><
>	That is true.  On the other hand, the system I am
>modeling is stable, consisting exclusively of damped sinusoids
>of various frequencies and additive noise (so why not use

Then a stability force can make sense.

>Prony's method?  I'm doing that, too).  Marple's book suggests

Sounds like a good idea.  Especially comparing results can be helpful.

>(and rightfully so) that if the results of doing AR analysis
>generates poles which have magnitude greater than one, attempts
>to model the system will be disappointing because of the
>instability in the model.  This may be easily confirmed upon

Not only that, but many algorithms that play with these models are numerically
unstable if the model is not stable.  So you get hosed from two sides :-)

>[...]  If the AR coefficients are the coefficients of
>polynomials whose roots are the poles are the system (we know
>they are), what are the consequences of scaling the
>coefficients?  

Scaling the coefficients (multiply the k-th coeff by alpha^k for some
0<alpha<1) has the effect of moving all the poles in by alpha (ie:

    i phi                   i phi
 r e        ---->  alpha r e

So geometric relations between the poles are preserved, but everything is
fitted into a smaller circle.

Another possibility is to calculate all the poles and move the unstable ones
just enough to stabilize them.  This sounds nice, but I have slightly horrible
experiences with it (and would be interested to hear from people who have good
ones).

Then there's calculating the partial covariances (or partial autocorrelations,
or whatever you want to call them today) and setting the magnitudes there.
Never tried this one.  Any takers?

>It seems to me that the frequency components
>should be preserved, and if |variance of n(t)| were also scaled
>down, the modeled X(t) sequence would be stable, but scaled
>incorrectly.  Is this wrong?  I would like to have stable
>estimates for this known-stable system, but most of the
>recommended AR techniques in Marple generate unstable
>parameters (Marple makes note of this possibility in a couple
>of places, but I didn't see any recommendations about how to
>fix the problem.  I *think* there's an implication that the
>"next better" algorithm should be used).

The general attitude seems to be that "as along as you futz only a little, it
doesn't matter".  Of course that doesn't help you much if you have to futz by
a lot...  I do not know of any good, general, rigorous treatment of this (very
real) problem (again, anybody out there know better?  Speak up please!)

I would say that all these ad-hoc methods are highly suspect if you have to
move any pole by more than a little.  At least, my experience has been that in
such cases you might as well throw away the estimate and try again.

There are methods that are guaranteed to give you stable estimates (they are
in Marple, I'm pretty sure), but the price you pay is bias in the estimates.
As your data-length increases, the bias goes down, but then, so does the
chance of getting an unstable estimate out of the other methods.  In short,
getting a good (or even reasonable) estimate out of short data records is
_hard_. 

If you are indeed using short records, you may want to try modifying your
algorithm to incorporate a SVD.  In both ARMA and Prony algorithms, you are
solving an equation of the form:

Ya = -z

Create the matrix [ z Y ] and do an SVD on it.  Say you get USV.  Now if you
are expecting the "correct" order to be n, keep n singular values.  Call this
partially zero-ed S matrix S'.  Compute a to be the LS solution of

Y'a = -z'

where 

[ z' Y' ] = US'V

The way you do this is start out with forming your Y and z for a much higher
order (even up to Number-of-data-points/2) and then use the SVD to go down to
order you want.  Of course your a polynomial will still be high order, ie have
a whole slew of spurious roots.  You need to get rid of those.  If the system
is known stable you can just throw away al unstable poles.  For backward
prediction that should take care of the spurious guys.

SR