@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
