Hi Ian
Congratulation and thanks for test and analysis!
The result changes from C=90 (correct) in the old ABC to C=89 (incorrect) in
the new ABC. Also the old ABC result B=41 and the new ABC result B=37 should
both be B=20. That is bewildering.
Consider the test data.
X=:0.2+0.4*133<i.250
Without noise the new ABC reproduces correctly
ABC X
0.2 0.4 133
Random rounding to one bit per data point
round=:>[:?0"0
gives fair approximations
ABC round X
0.148556 0.451145 131
ABC round X
0.113223 0.51251 124
ABC round X
0.209905 0.344298 136
ABC round X
0.189825 0.393188 136
ABC round X
0.215108 0.418674 134
ABC round X
0.247703 0.401055 133
(This should be done using 'for.', but I don't know how!)
I expect that FFTW is much faster than SFT for bigger values of N.
The three fourier transform verbs have these characteristics.
1000 transform the H-vectors one at a time ("1)
2000 transform the whole H-matrix at once
0100 slow matrix multiplication (+/ . *)
0200 fast matrix multiplication (FFT tricks)
0010 do not normalize the U-matrix
0020 normalize the U-matrix
0001 no complex conjugation
0002 complex conjugation
0022 the transform equals the invers transform
1122 FT, old ABC
2122 SFT, new ABC
1211 FFTW
We need 2222: transform the whole H-matrix at once, use the fast matrix
multiplication, normalize, and conjugate.
require 'math/fftw'
|file name error: script
| 0!:0 y[4!:55<'y'
How is it done?
- Bo
>________________________________
> Fra: Ian Clark <[email protected]>
>Til: Programming forum <[email protected]>
>Sendt: 22:22 mandag den 2. juli 2012
>Emne: Re: [Jprogramming] Kalman filter in J?
>
>Tried out both your old and new ABC in my rig, with your original FT,
>with SFT and with fftw.
>
>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
>
> NB. Old ABC ...
>
> timer 'smoutput ABC X' [ FTS=:1 NB. old FT
>98.6649 41.0013 90
>10.1967
> timer 'smoutput ABC X' [ FTS=:2 NB. new SFT
>98.6649 41.0013 90
>2.20782
> timer 'smoutput ABC X' [ FTS=:3 NB. fftw (adjusted to be
>self-inverse)
>98.6649 41.0013 90
>0.0496461
> timer 'smoutput ABC X' [ FTS=:4 NB. fftw (raw)
>19535.7 8118.25 90
>0.0446591
>
>(fftw means, I understand, "the fastest Fourier transform in the West" :-)
>
> NB. New "shined up" ABC ...
>
> timer 'smoutput ABC X' [ FTS=:1
>98.6542 36.7969 89
>5.09876
> timer 'smoutput ABC X' [ FTS=:2
>98.6542 36.7969 89
>0.0406119
> timer 'smoutput ABC X' [ FTS=:3
>98.6542 36.7969 89
>0.02442
> timer 'smoutput ABC X' [ FTS=:4
>19533.5 7285.79 89
>0.022739
>
>The speed gain when using fftw (FTS=3 and 4) is considerable, eg 2.2
>--> 0.05 with the old ABC. There is very little extra gain from using
>fftw "raw", ie without the adjustment to make FT self-inverse. Without
>this adjustment (FTS=4), ABC returns A and B multiplied by unwanted
>factors, but still estimates C to be the same as before.
>
>Even without fftw, the new ABC is much faster, eg 2.2 --> 0.04, even
>when using SFT. In this instance the further improvement from using
>fftw is not all that spectacular, eg 0.04 --> 0.02. All variants seem
>to estimate A and C quite well, but with this instance of X, B is
>overestimated. However it usually gets B closer.
>
>The overhead of using smoutput within timer is small (approx 0.005).
>
>The flag: FTS determines the verb actually called by FT, which I've
>arranged to be called throughout both versions of ABC. In fact FT
>calls ft<n> for FTS=:<n> .
>
>require 'math/fftw'
>ft1=: +/&(+*((%:%~_1&^&+:&(%~(*/])&i.))&#))
>U=: %:%~_1^(%~+:&(*/~)&i.) NB. unitary mx
>ft2=: ++/ . *U&{:&$
>ft3=: [: + fftw % [: %: #
>ft4=: fftw
>FT=: 3 : '(ft=: 0:`(ft1"1)`ft2`(ft3"1)`(ft4"1) @. FTS) y'
>
>(fftw happens to return the complex conjugate of your SFT, hence '[:
>+' in the defn of ft3.)
>
>Note that ft1, ft3 and ft4 need to be suffixed by "1 to make them
>behave like ft2 when acting on a 2D matrix: y, ie to make them
>transform each row of y independently.
>
>My suspicion so far is that, with fftw, Bo's FT step-detector may well
>out-perform the other statistics-based methods I've tried to date. I
>guess this is because it uses all the available information from X,
>making no use of confidence limits or short sample means (which waste
>information), nor even thresholds for acceptable false positives /
>negatives.
>
>I shall write this up in my "experimental report" in the wiki, plus
>the code to reproduce my results.
>
>Ian
>
>
>
>On Sat, Jun 30, 2012 at 4:59 PM, Bo Jacoby <[email protected]> wrote:
>> Hi Ian.
>> I shined the program up a bit. Rather than transforming each row I tranform
>> the whole matrix at once. It is faster, but still slow. The test function
>> used is rounded down to zero or up to one. Of course it is hard to
>> reconstruct ABC from a bitvector, but I think it does a reasonable job.
>> Please set me know how fast is is on your test examples, and if you manage
>> to replace the Slow Fourier Transform (SFT) by Fast Fourier Transform (FFT).
>> - Bo
>>
>>
>> NB.H =: Heaviside steps
>> H =: i.&<:</i.
>> NB.U =: Unitary matrix generator
>> U =: %:%~_1^(%~+:&(*/~)&i.)
>> NB.SFT=: Slow Fourier Transformation
>> SFT=: ++/ . *U&{:&$
>> NB.sm =: symmetrize
>> sm =: [:,"1/}:"1&(,:|."1)
>> NB.RFT=: Real Fourier Transform.
>> RFT=: 9 o.${.SFT&sm
>> NB.M1 =: zero origin index of the last item before the step
>> M1 =: [:(i.<./)[:+/"1[:*:[:}.[:(-"1{.)[:}."1[:(%{."1)}."1
>> M2 =: [:{:"1[:RFT{.,:~([:{:$){.[:{.{:
>> M3 =: [:(}.,-/){.,:([:%/1&{"1)*1&{
>> M4 =: M2&M3
>> M5 =: 3 : 'x,~M4 y{~0,>:x=.M1 y'
>> M6 =: RFT&(,H&{:&$)
>> NB.ABC=: constant,stepsize,stepposition
>> ABC=: M5&M6
>> round=:>[:?0$~$
>> ABC X=.round 0.2+0.4*33<i.200
>> 0.264787 0.273477 31
>> 10 20$X
>> 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 0 1 0 0
>> 0 0 0 1 0 0 0 0 1 0 0 0 1 0 1 1 1 1 0 1
>> 1 0 1 1 0 0 1 1 1 0 1 0 1 1 0 1 0 0 1 0
>> 1 0 1 1 1 1 0 0 0 1 0 0 0 1 0 1 0 1 1 1
>> 0 1 1 1 0 0 1 1 0 0 0 1 1 0 1 0 0 1 0 0
>> 0 1 1 0 0 0 1 0 1 0 1 1 1 1 1 0 1 0 0 1
>> 0 1 0 0 0 0 1 1 1 1 1 1 0 1 1 0 1 0 1 1
>> 0 0 1 1 0 1 0 0 1 1 0 1 0 1 1 1 1 0 1 1
>> 1 0 1 0 0 0 1 1 0 1 1 0 1 0 0 0 1 0 1 0
>> 0 1 0 1 0 0 1 1 0 0 1 1 1 0 1 1 1 0 0 1
>>
>>
>...snipped
>----------------------------------------------------------------------
>For information about J forums see http://www.jsoftware.com/forums.htm
>
>
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
