Does anyone know how to efficiently implement a recurrence relationship in numpy such as:
y[n] = A*x[n] + B*y[n1]
Thanks,
Gerry
On Wed, May 6, 2009 at 6:44 AM, Talbot, Gerry Gerry.Talbot@amd.com wrote:
Does anyone know how to efficiently implement a recurrence relationship in numpy such as:
y[n] = A*x[n] + B*y[n1]
On an intel chip I'd use a Monte Carlo simulation. On an amd chip I'd use:
x = np.array([1,2,3]) y = np.array([4,5,6]) y = x[1:] + y[:1] y
array([6, 8])
Sorry, I guess I wasn't clear, I meant:
for n in xrange(1,N): y[n] = A*x[n] + B*y[n1]
So y[n1] is the result from the previous loop iteration.
Gerry
Original Message From: numpydiscussionbounces@scipy.org [mailto:numpydiscussionbounces@scipy.org] On Behalf Of Keith Goodman Sent: Wednesday, May 06, 2009 9:53 AM To: Discussion of Numerical Python Subject: Re: [Numpydiscussion] Recurrence relationships
On Wed, May 6, 2009 at 6:44 AM, Talbot, Gerry Gerry.Talbot@amd.com wrote:
Does anyone know how to efficiently implement a recurrence relationship in numpy such as:
y[n] = A*x[n] + B*y[n1]
On an intel chip I'd use a Monte Carlo simulation. On an amd chip I'd use:
x = np.array([1,2,3]) y = np.array([4,5,6]) y = x[1:] + y[:1] y
array([6, 8]) _______________________________________________ Numpydiscussion mailing list Numpydiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
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[n1]
So y[n1] 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.
Josef
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[n1]
So y[n1] 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
Le mercredi 06 mai 2009 à 10:21 0400, josef.pktd@gmail.com a écrit :
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[n1]
So y[n1] 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.
Josef
Isn't it what scipy.signal.lfilter does ? y=scipy.signal.lfilter([A],[1,B],x)
You may be careful with initial conditions...
The application is essentially filtering 1D arrays, typically N is
20e6, the required result is y[1:N].
Gerry
Original Message From: numpydiscussionbounces@scipy.org [mailto:numpydiscussionbounces@scipy.org] On Behalf Of Alan G Isaac Sent: Wednesday, May 06, 2009 10:25 AM To: Discussion of Numerical Python Subject: Re: [Numpydiscussion] Recurrence relationships
On 5/6/2009 10:00 AM Talbot, Gerry apparently wrote:
for n in xrange(1,N): y[n] = A*x[n] + B*y[n1]
So, x is known before you start? How big is N? Also, is y.shape (N,)? Do you need all of y or only y[N]?
Alan Isaac
_______________________________________________ Numpydiscussion mailing list Numpydiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
On Wed, May 6, 2009 at 10:44 PM, Talbot, Gerry Gerry.Talbot@amd.com wrote:
Does anyone know how to efficiently implement a recurrence relationship in numpy such as:
y[n] = A*x[n] + B*y[n1]
That's the direct implement of a linear filter with an infinite impulse response. That's exactly what scipy.signal.lfilter is for,
cheers,
David
participants (6)

Alan G Isaac

David Cournapeau

Fabrice Silva

josef.pktd＠gmail.com

Keith Goodman

Talbot, Gerry