Hi Ian.
I'm glad that your test was successful, albeit slow. The important optimization
is to replace my slow fourier transform by a fast fourier transform. Your
example has 100 items in the time series, and the real transform RT asks the
fourier transform FT to perform 98 vector multiplications by a 198*198 complex
matrix, (98*198*198)=3841992, which is embarassingly inefficient. The FFT will
speed it up by a factor ((%2&^.)198)=26.
Another optimization is to incorporate the computation of A and B into the
computation of C, rather than starting from scratch after computing C.
The choice of index origin for C is a matter of taste.
Why would a Hartley transform be faster than the Fast Fourier Transform?
It remains to estimate the inaccuracies of A and B. It is the root-mean-square
value of the noise vector X-A+B*C<#X. Note that the noise has only (#X)-2
degrees of freedom, because A and B have been estimated, so the
root-mean-square is
RMS=:(+/%#-2:)&.:*:
- Bo
>________________________________
> Fra: Ian Clark <earthspo...@gmail.com>
>Til: Programming forum <programming@jsoftware.com>
>Sendt: 2:41 tirsdag den 26. juni 2012
>Emne: Re: [Jprogramming] Kalman filter in J?
>
>...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
>
>
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
