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 <bojac...@yahoo.dk> 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
