@Bo -Here is your solution running in my test rig... A= 100 bias constant (popn mean of X before the step) B= 20 scaling factor for height of step C= 90 epoch of actual occurrence of the Heaviside step D= 10 scaling factor for Gaussian noise component (sigma) N= 100 number of epochs in time series X step detected at: (wait...) 90 ABC X 99.8246 20.0353 90
ABC takes approx 10 s to run on my iMac. CC (--your 'C') takes approx 5 s. I'm actually quite impressed with how accurately it estimates A, B, C, compared with other methods (which don't do A or B). I'm writing up my experiments and will post a wiki link. Whilst I like your symmetrisation idea, might instead using the discrete Hartley Transform (which works entirely in the real domain) speed things up? Ian On Mon, Jun 25, 2012 at 12:00 PM, Bo Jacoby <bojac...@yahoo.dk> wrote: > Hi Ian. Here is your solution. > > ABC=: 3 : 0 > c=.C y > y=.(0,>:c){RT"1&(,H&#)y > y=.y,(%/1{"1 y)*1{y > y=.y,(0{y)-2{y > y=.y,(#|:y){.{.{:y > y=.4 2{{:"1 RT"1 y > y,c > ) > ABC 101 103 100 210 230 200 > 101.15 115.101 2 > > > - Bo > > > > >>________________________________ >> Fra: Ian Clark <earthspo...@gmail.com> >>Til: Programming forum <programming@jsoftware.com> >>Sendt: 6:46 mandag den 25. juni 2012 >>Emne: Re: [Jprogramming] Kalman filter in J? >> >>Interesting approach, Bo. >> >>Turned my attention to CUSUM (an adaptive filter) which is rather easy >>to compute and gives good results, but needs parameters setting "by >>eye". I've yet to try a Bayesian approach, which I suspect will be the >>most sensitive of all. >> >>But I did play with your FT, using *:@(10&o.) to plot a power >>spectrum. I'm using rather noisier data than your example: my step is >>only roughly (sigma) high. Nevertheless with the onset of the step I >>can clearly see a high-frequency spike appear at the far end of the >>transformed time-series, due to the transformed Heaviside fn. Same >>problem as with the CUSUM statistic however: designing a detector >>which can be adjusted for false negatives and positives, and doesn't >>take too many subsequent samples to detect the step. >> >>Will try out your C function in my context, and try to tweak it. Your >>solution detects when the step happens, solving 1 in my stated >>objectives: 1-4. I fear however that I am more interested in 2-4. >> >> >> >>On Sun, Jun 24, 2012 at 6:06 PM, Bo Jacoby <bojac...@yahoo.dk> wrote: >>> 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 >>---------------------------------------------------------------------- >>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