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
