...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
