Hi Ian The Fourier Transform FT transforms a complex n-vector X into a complex n-vector Y.If X is real, then Y is symmetric. And if X is symmetric then Y is real. So I think it is a good idea to symmetrize X before taking the FT. (This makes a slow program even slower. Optimize later!)
NB.sm=: symmetrize sm=:[:,}:"1&(,:|.) sm i.6 NB. see what I mean 0 1 2 3 4 5 4 3 2 1 NB.RT=: Real Transform. (9 o. removes every trace of complexity) RT=:9 o.#{.FT&sm NB.H=:Heaviside steps H=:i.&<:</i. NB.C=:Index of the last item before the change C=:[:(i.<./)[:+/"1[:*:[:}.[:(-"1{.)[:}."1[:(%{."1)[:}."1[:RT"1],[:H# NB.test C 101 204 202 200 203 0 C 101 104 202 200 203 1 C 101 104 102 200 203 2 C 101 104 102 100 203 3 This is not a Kálmán filter, but it solves the problem. - Bo >________________________________ > Fra: Ian Clark <earthspo...@gmail.com> >Til: Programming forum <programming@jsoftware.com> >Sendt: 3:20 fredag den 22. juni 2012 >Emne: Re: [Jprogramming] Kalman filter in J? > >Thanks, Bo. I'll look at it when I get back, probably over the weekend. > >Ian > >On Fri, Jun 22, 2012 at 12:44 AM, Bo Jacoby <bojac...@yahoo.dk> wrote: >> Hi Ian. >> The following is a tool, but not yet a solution. >> >> FT is a slow Fourier Transform. When the problem has been solved it is time >> to optimize. This FT has the nice property that FT^:_1 is identical to FT. >> >> NB.FT =: Fourier Transformation >> FT =: +/&(+*((%:%~_1&^&+:&(%~(*/])&i.))&#)) >> NB.rd =: round complex number. >> rd =: (**|)&.+. >> NB.H =: Heaviside step function >> H =: (]#0:),-#1: >> 10 H 3 >> 0 0 0 1 1 1 1 1 1 1 >> NB. FT = FT^:_1 >> rd FT FT 10 H 3 >> 0 0 0 1 1 1 1 1 1 1 >> >> >> NB. For fixed C the values of A and B to minimize +/*:X-A+B*(#X)H C are >> found by solving (2{.FT X) = 2{.FT A+B*(#X)H C >> >> - Bo >>>________________________________ >>> Fra: Ian Clark <earthspo...@gmail.com> >>>Til: Programming forum <programming@jsoftware.com> >>>Sendt: 18:15 torsdag den 21. juni 2012 >>>Emne: Re: [Jprogramming] Kalman filter in J? >>> >>>Yes, Bo, that looks like a fair statement of the problem, in terms of >>>least-squares minimization. >>> >>> >>>On Thu, Jun 21, 2012 at 4:46 PM, Bo Jacoby <bojac...@yahoo.dk> wrote: >>>> Consider the step functions >>>> >>>> H=:(]#0:),-#1: >>>> 10 H 0 >>>> 1 1 1 1 1 1 1 1 1 1 >>>> 10 H 3 >>>> 0 0 0 1 1 1 1 1 1 1 >>>> 10 H 10 >>>> 0 0 0 0 0 0 0 0 0 0 >>>> The problem is to find constants A,B,C to minimize the expression >>>> >>>> +/*:X-A+B*(#X)H C >>>> Right? >>>> -Bo >>>>>________________________________ >>>>> Fra: Ian Clark <earthspo...@gmail.com> >>>>>Til: Programming forum <programming@jsoftware.com> >>>>>Sendt: 14:34 torsdag den 21. juni 2012 >>>>>Emne: Re: [Jprogramming] Kalman filter in J? >>>>> >>>>>Thanks, Bo and Ric. Yes, I'd tried searching jwiki on "Kalman" and >>>>>Pieter's page was the only instance. >>>>> >>>>>I neglected to mention that in the most general case I can't really be >>>>>confident that X1 and X2 have the same distribution function, let >>>>>alone the same variance. But looking at it again, I see that under the >>>>>restrictions I've placed the problem simplifies immensely to fitting a >>>>>step-function H to: X=: K+G+H. If I can just do that, I'll be happy >>>>>for now. >>>>> >>>>>Repeated application of FFT should allow me to subtract the noise >>>>>spectrum F(G), or at least see a significant change in the overall >>>>>spectrum emerge after point T, and that might handle the more general >>>>>cases as well. >>>>> >>>>>Anyway it's simple-minded enough for me, and worth a try. FFT is, >>>>>after all, "fast" :-) >>>>> >>>>>But won't an even faster transform do the trick, such as (+/X)? On the >>>>>above model, X performs a drunkard's walk around a value M1 until some >>>>>point T, after which it walks around M2. Solution: simply estimate M1 >>>>>and M2 on an ongoing basis. >>>>> >>>>>I get the feeling I ought to be searching on terms like "edge >>>>>detection", "step detection" and CUSUM. >>>>> >>>>>Anyway, there's enough here to try. >>>>> >>>>>On Thu, Jun 21, 2012 at 10:28 AM, Ric Sherlock <tikk...@gmail.com> wrote: >>>>>> Hi Ian, >>>>>> >>>>>> A quick search of the J wiki finds this: >>>>>> http://www.jsoftware.com/jwiki/Stories/PietdeJong >>>>>> Sounds like he might have what you're after? >>>>>> >>>>>> Cheers, >>>>>> Ric >>>>>> >>>>>> On Thu, Jun 21, 2012 at 4:33 PM, Ian Clark <earthspo...@gmail.com> wrote: >>>>>>> Can anyone help? Has anyone written a Kalman filter in J? >>>>>>> >>>>>>> I'm not a specialist in either statistics or control theory, so I'm >>>>>>> only guessing a Kalman filter is what I need. Though I do have a >>>>>>> passing acquaintance with the terms: stochastic control and linear >>>>>>> quadratic Gaussian (LQG) control. I am aware that a "Kalman filter" >>>>>>> (like ANOVA) is more a topic than a black-box. >>>>>>> >>>>>>> So let me explain what I want it for. >>>>>>> >>>>>>> I have a time series X which I am assuming can be modelled like this: >>>>>>> >>>>>>> X=: K + G + (X1,X2) >>>>>>> >>>>>>> where >>>>>>> >>>>>>> K is constant >>>>>>> G is Gaussian noise >>>>>>> X1 is a random variable with mean: M1 and variance: V1 >>>>>>> X2 is a random variable with mean: M2 and variance: V2 >>>>>>> >>>>>>> Typically X is a sequence of sensor readings, but may also be >>>>>>> measurements from a series of user trials conducted on a working >>>>>>> prototype, which suffers a design-change at a given point T. >>>>>>> >>>>>>> Simplifying assumptions (which unfortunately I may need to relax in due >>>>>>> course): >>>>>>> >>>>>>> (a) X is not multivariate >>>>>>> (b) X1 and X2 are Gaussian >>>>>>> (c) V1=V2 (only the mean value changes, not the variance). >>>>>>> >>>>>>> The problem: >>>>>>> >>>>>>> 1. Estimate T=: 1+#X1 -- the point at which X1 gives way to X2. >>>>>>> >>>>>>> 2. Given T, estimate (M2-M1) -- the "underlying improvement", if any, >>>>>>> of the change to the prototype. >>>>>>> >>>>>>> 3. (subcase of 2.) Given T, test the null hypothesis: M1=M2, viz that >>>>>>> there has been no underlying improvement. >>>>>>> >>>>>>> 4. Estimate U=: #X2 -- the minimum number of samples needed after T in >>>>>>> order to achieve 1-3 above with 95% confidence. >>>>>>> >>>>>>> In other words, detect the signal-in-noise: M1-->M2, and do so in >>>>>>> real-time. >>>>>>> >>>>>>> Because of 4, the need to estimate T and (M2-M1) on an ongoing basis, >>>>>>> I can't do a randomised block design. I gather that a Kalman filter, >>>>>>> or some sort of adaptive filter, will handle this problem. >>>>>>> >>>>>>> But maybe something simpler will turn out good enough? >>>>>>> >>>>>>> Supposing I can get hold of a "black box" Kalman filter, I propose to >>>>>>> test it out on generated data and compare its performance to some >>>>>>> simple-minded approach, like estimating M1 / M2 from a simple moving >>>>>>> average of the last U samples, or applying the F-test to 2 sets of U >>>>>>> samples taken either side of T. >>>>>>> >>>>>>> But since the technique aims to be published, or at least critically >>>>>>> scrutinised (and maybe incorporated in a software product), I'd rather >>>>>>> depend on a state-of-art packaged solution than reinvent the wheel: a >>>>>>> large and very well-turned wheel it appears to me. >>>>>>> >>>>>>> Ian Clark >>>>>>> ---------------------------------------------------------------------- >>>>>>> For information about J forums see http://www.jsoftware.com/forums.htm >>>>>---------------------------------------------------------------------- >>>>>For information about J forums see http://www.jsoftware.com/forums.htm >>>>> >>>>> >>>>> >>>> ---------------------------------------------------------------------- >>>> For information about J forums see http://www.jsoftware.com/forums.htm >>>---------------------------------------------------------------------- >>>For information about J forums see http://www.jsoftware.com/forums.htm >>> >>> >>> >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forums.htm >---------------------------------------------------------------------- >For information about J forums see http://www.jsoftware.com/forums.htm > > > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm