Re: [Gretl-users] Loop recursion: Speed comparison with Matlab
On Tue, 24 Oct 2017, Sven Schreiber wrote:
[not sure if gretl-users is the right place for this...]
Am 24.10.2017 um 09:01 schrieb Artur Tarassow:
> Dear all,
>
> I made a little experiment this morning as I was a bit puzzled that the
> SB() function ...
Artur, first a thing that hasn't to do with speed. Are you sure the following
code from SB() is doing what it's intended to do:
<SB>
if randgen1(u,0,1)<1/b
u[t] = ceil(n*randgen1(u,0,1))
</SB>
Because each call to randgen1() should be giving you a new random number, if
you know what I mean. In the original code (now in SBslow()) this looks
different.
Just checking my understanding: the algorithm as given in current SB()
for generating the vector of indices, u, into the n-vector of data
using block size b is:
* Set u[1] to a random point in [1,n].
for t = 2 to n
* With probability 1/b, set u[t] to a new random point in [1,n], or
* with probability 1 - 1/b, set u[t] to u[t-1] + 1.
If that's right, here's faster code:
<hansl>
function matrix SB (const matrix x "Data to be resampled",
int b[0::] "Block size")
scalar n = rows(x)
if n == 0
funcerr "data input error, check data"
endif
if b == 0
b = round(1.75*(n^(1/3)))
endif
matrix u = zeros(n,1)
# start resampling
u[1] = randint(1,n)
matrix cond = muniform(n,1) .< 1/b
loop t=2..n -q
u[t] = cond[t] ? randint(1,n) : u[t-1] + 1
endloop
return (x|x)[u]
end function
</hansl>
On my system, that reduces the Example 1 time from about 80 seconds to
about 40. Octave 4.2.1, BTW, takes over 12 minutes to run the Matlab
version.
Allin