🚀 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

In article <1456@mrsvr.UUCP> kohli@gemed.ge.com (Mr. Bad Judgment) writes:
[...previous discussion deleted]
>	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
>Prony's method?  I'm doing that, too).  Marple's book suggests
[]
>	X(t) = A(t)X(t-1) + n(t).
>
>	X(t) may be expected to grow without bound if |A(t)| >
>1.0.  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?  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
>Opinions?  Facts?
>Thanks again,
>Jim Kohli
>GE Medical Systems
>P.S.  Thank you Theo, Rob, and Stephen, the discussion has been
>interesting and enlightening (so far!).

I assume by damped sinusoids you mean sinusoids with exponentially decaying
amplitude? (Are you by any chance doing MRI data? We've done a lot of
experimentation on modeling MRI)
Let's assume that we have an exponentially decaying unit step. Then the time
series of the output is:
			        n
			x(n) = a  , a < 1.0, 0 <= n < inf.

The AR(1) model is:	       -1
			1 + a z

and if a < 1.0, the thing is stable. However, if we look at the thing 
backwards, we see an exponentially increasing series. If we model this
sequence, we get an AR coefficient of 1/a*, that is, outside the unit circle.
The coefficient is also conjugated from the forward model coefficient (I
think. I'm doing this from memory), so the overall angle of the pole in the
Z-plane stays the same, it's just been reflected outside the unit circle.
Obviously, this situation is unstable.
If your algorithm applies forward _and_ backward prediction to decaying data,
your results will not be correct, since the model cannot provide a stable
solution in both directions. The coefficients will be the result of a 
best-fit in both directions, and will not mean much of anything. This
is a problem I found in using the Burg algorithm; anything that's not
stationary doesn't get modeled worth &*@Q$.

About the pole radius scaling - if you don't care what the absolute amplitude
is (ie in spectral analysis, you're usually only interested in relative
magnitude), there is no problem. The pole scaling results in reduced
resolution in the frequency domain, similar to what windowing the data
does. If you're trying to model the data exactly enough to be able to
recreate it, a different approach will be necessary.

Hope this helps,

Theo Smit (smit@enel.ucalgary.ca)