...I should add that my CC returns 1+(the value of your 'C'). This is for compatibility with how I generate the step fn, using (the constant) C.
On Tue, Jun 26, 2012 at 1:37 AM, Ian Clark <earthspo...@gmail.com> wrote: > @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