[Numpy-discussion] Recurrence relationships
josef.pktd@gmai...
josef.pktd@gmai...
Wed May 6 09:28:46 CDT 2009
On Wed, May 6, 2009 at 10:21 AM, <josef.pktd@gmail.com> wrote:
> On Wed, May 6, 2009 at 10:00 AM, Talbot, Gerry <Gerry.Talbot@amd.com> wrote:
>> Sorry, I guess I wasn't clear, I meant:
>>
>> for n in xrange(1,N):
>> y[n] = A*x[n] + B*y[n-1]
>>
>> So y[n-1] is the result from the previous loop iteration.
>>
>
> I was using scipy.signal for this but I have to look up what I did
> exactly. I think either signal.correlate or using signal.lti.
>
No, its signal.lfilter, below is a part of a script I used to simulate
and estimate an AR(1) process, which is similar to your example.
I haven't looked at it in a while but it might give you the general idea.
Josef
# Simulate AR(1)
#--------------
# ar * y = ma * eta
ar = [1, -0.8]
ma = [1.0]
# generate AR data
eta = 0.1 * np.random.randn(1000)
yar1 = signal.lfilter(ar, ma, eta)
etahat = signal.lfilter(ma, ar, y)
np.all(etahat == eta)
# find error for given filter on data
print 'AR(2)'
for rho in [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.79, 0.8, 0.81, 0.9]:
etahatr = signal.lfilter(ma, [1, --rho], yar1)
print rho,np.sum(etahatr*etahatr)
print 'AR(2)'
for rho2 in np.linspace(-0.4,0.4,9):
etahatr = signal.lfilter(ma, [1, -0.8, -rho2], yar1)
print rho2,np.sum(etahatr*etahatr)
def errfn(rho):
etahatr = signal.lfilter(ma, [1, -rho], yar1)
#print rho,np.sum(etahatr*etahatr)
return etahatr
def errssfn(rho):
etahatr = signal.lfilter(ma, [1, -rho], yar1)
return np.sum(etahatr*etahatr)
resultls = optimize.leastsq(errfn,[0.5])
print 'LS ARMA(1,0)', resultls
resultfmin = optimize.fmin(errssfn, 0.5)
print 'fminLS ARMA(1,0)', resultfmin
More information about the Numpy-discussion
mailing list