6:29 a.m.
I agree Artur. I was recently in need for a resample without replacement
function.
I have not check your code but I had a look on what I did and it could
be a refinement (following also Sven's answer).
Say a matrix A is nxm.
set k<=m.
rvec = ranking( mrandgen(u,0,1,m,1) )'
A=[,rvec[1:k]]
(keeps k randomly selectled columns of A)
I would like to see such a function in extra.gfn
Yiannis
--
**************************************************************************************************
Ioannis A. Venetis, Associate Professor
University of Patras, School of Economic Sciences and Business
Administration,
Department of Economics
University Campus, Rio, 26504, Greece Tel: 0030 2610 969964
http://www.econ.upatras.gr/en/people/teaching-staff/venetis-ioannis
[1]https://sites.google.com/site/ioannisavenetis/
[2]**************************************************************************************************
Links:
------
[1] http://application.econ.upatras.gr/en/node/130
[2] https://sites.google.com/site/ioannisavenetis/
7:47 a.m.
On Tue, 17 Sep 2019, Ioannis A. Venetis wrote:
I agree Artur. I was recently in need for a resample without
replacement
function.
I have not check your code but I had a look on what I did and it could
be a refinement (following also Sven's answer).
Say a matrix A is nxm.
set k<=m.
rvec = ranking( mrandgen(u,0,1,m,1) )'
A=[,rvec[1:k]]
(keeps k randomly selectled columns of A)
I would like to see such a function in extra.gfn
I like it. We could have something like this in extra.gfn:
<hansl>
function matrix sample_wo_rep(matrix X, scalar k)
scalar n = rows(X)
if k > n
funcerr "what?"
endif
matrix sel = ranking(muniform(n,1))
return X[sel[1:k],]
end function
### example
X = mnormal(10, 3)
print X
eval sample_wo_rep(X, 6)
</hansl>
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
9:19 a.m.
On Tue, 17 Sep 2019, Sven Schreiber wrote:
Am 17.09.2019 um 09:47 schrieb Riccardo (Jack) Lucchetti:
>
> matrix sel = ranking(muniform(n,1))
> return X[sel[1:k],]
>
the doc for ranking says "...plus one half the number of elements...",
so it can be non-integer, no?
In theory, yes. But the odds of having a tie among pseudo-random uniforms
are so small I wouldn't worry about it.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
9:32 a.m.
Am 17.09.2019 um 11:19 schrieb Riccardo (Jack) Lucchetti:
On Tue, 17 Sep 2019, Sven Schreiber wrote:
> Am 17.09.2019 um 09:47 schrieb Riccardo (Jack) Lucchetti:
>>
>> matrix sel = ranking(muniform(n,1))
>> return X[sel[1:k],]
>>
> the doc for ranking says "...plus one half the number of elements...",
> so it can be non-integer, no?
In theory, yes. But the odds of having a tie among pseudo-random
uniforms are so small I wouldn't worry about it.
Hm, I don't know. Small, of course, but if millions of future gretl
users start running huge simulations based on that, Murphy's law says
there's going to be an error. As you know, the probability isn't zero
for computer-represented numbers.
But anyway, doesn't ranking involve expensive sorting? My gut feeling is
that using mrandgen(i,...) to produce "more" than enough random indices,
then simply discarding the duplicates with uniq(), could be faster. (Yes
you need an additional check whether enough indices are left over, and
you might have to do extra work, but most of the time it wouldn't be
required.)
cheers
sven
11:03 a.m.
On Tue, 17 Sep 2019, Sven Schreiber wrote:
Am 17.09.2019 um 11:19 schrieb Riccardo (Jack) Lucchetti:
> On Tue, 17 Sep 2019, Sven Schreiber wrote:
>
>> Am 17.09.2019 um 09:47 schrieb Riccardo (Jack) Lucchetti:
>>>
>>> matrix sel = ranking(muniform(n,1))
>>> return X[sel[1:k],]
>>>
>> the doc for ranking says "...plus one half the number of elements...",
>> so it can be non-integer, no?
>
> In theory, yes. But the odds of having a tie among pseudo-random
> uniforms are so small I wouldn't worry about it.
Hm, I don't know. Small, of course, but if millions of future gretl
users start running huge simulations based on that, Murphy's law says
there's going to be an error. As you know, the probability isn't zero
for computer-represented numbers.
Yes. But doubles are 8 byte large. But the probability of extracting k
different numbers out of n choices without repetition is
$\frac{(n-1)!}{(n-k)! n^{k-1}}$ (if I recall correctly), so when n = 2^32
I guess we're pretty safe for any reasonable value of k.
But anyway, doesn't ranking involve expensive sorting? My gut
feeling is
that using mrandgen(i,...) to produce "more" than enough random indices,
then simply discarding the duplicates with uniq(), could be faster. (Yes
you need an additional check whether enough indices are left over, and
you might have to do extra work, but most of the time it wouldn't be
required.)
Sorting via qsort (which is what we use) is normally O(n log n), so I
guess it's about the same speed as uniq(). Maybe it could be interesting
to write a speed comparison script.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
11:07 a.m.
yOn Tue, 17 Sep 2019, Riccardo (Jack) Lucchetti wrote:
Yes. But doubles are 8 byte large. But the probability of extracting
k
different numbers out of n choices without repetition is
$\frac{(n-1)!}{(n-k)! n^{k-1}}$ (if I recall correctly), so when n = 2^32 I
guess we're pretty safe for any reasonable value of k.
Whoops, sorry, it's 2^64.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
7:13 p.m.
On Tue, 17 Sep 2019, Sven Schreiber wrote:
[D]oesn't ranking involve expensive sorting? My gut feeling is
that
using mrandgen(i,...) to produce "more" than enough random indices,
then simply discarding the duplicates with uniq(), could be faster.
Here's what we do for datasets with "smpl n --random". We first see
which is smaller, the number of observations to be selected or the
number to be dropped. Working with the smaller, we generate random
integers in the appropriate range to select observations to include or
skip, discarding duplicates and proceeding till we have enough.
I've put a matrix version of this into git. It's provisionally called
msample and it takes a matrix and the scalar n as arguments.
matrix X = seq(1,100)'
eval msample(X, 20)
But one further point occurs to me: When people run such a function
does it meet expectation if the observations/rows are a random subset,
but in their original order? Or would users expect the order to be
shuffled? (Of course it's easy enough to shuffle the rows of a matrix
yourself if you want.)
Allin
9:12 p.m.
Am 18.09.2019 um 21:13 schrieb Allin Cottrell:
But one further point occurs to me: When people run such a function
does it meet expectation if the observations/rows are a random subset,
but in their original order? Or would users expect the order to be
shuffled? (Of course it's easy enough to shuffle the rows of a matrix
yourself if you want.)
No, drawing without replacement should not be order-preserving, just
like doing it with replacement. The last index could be drawn first and
then should come first, and so on. (If that is actually your question.)
cheers
sven
9:48 p.m.
On Wed, 18 Sep 2019, Sven Schreiber wrote:
Am 18.09.2019 um 21:13 schrieb Allin Cottrell:
> But one further point occurs to me: When people run such a function
> does it meet expectation if the observations/rows are a random subset,
> but in their original order? Or would users expect the order to be
> shuffled? (Of course it's easy enough to shuffle the rows of a matrix
> yourself if you want.)
No, drawing without replacement should not be order-preserving, just
like doing it with replacement. The last index could be drawn first and
then should come first, and so on.
OK -- but if the order of the observations doesn't matter, as in a
cross-sectional dataset, it's computationally cheaper to build a
random array of selected rows then extract them in order.
Anyway, the (still provisional) built-in msample() function now (in
git) shuffles the order by default but takes a third, optional
boolean argument to preserve order: if it has a non-zero value the
selected rows appear in their original order.
Any other thoughts on this?
Allin
6:45 a.m.
On Wed, 18 Sep 2019, Allin Cottrell wrote:
On Wed, 18 Sep 2019, Sven Schreiber wrote:
> Am 18.09.2019 um 21:13 schrieb Allin Cottrell:
>
>> But one further point occurs to me: When people run such a function
>> does it meet expectation if the observations/rows are a random subset,
>> but in their original order? Or would users expect the order to be
>> shuffled? (Of course it's easy enough to shuffle the rows of a matrix
>> yourself if you want.)
>
> No, drawing without replacement should not be order-preserving, just
> like doing it with replacement. The last index could be drawn first and
> then should come first, and so on.
OK -- but if the order of the observations doesn't matter, as in a
cross-sectional dataset, it's computationally cheaper to build a random array
of selected rows then extract them in order.
Anyway, the (still provisional) built-in msample() function now (in git)
shuffles the order by default but takes a third, optional boolean argument to
preserve order: if it has a non-zero value the selected rows appear in their
original order.
Any other thoughts on this?
I personally don't like very much having both msample() and resample(),
but I can't see an easy way out of this.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
9:02 a.m.
Am 19.09.2019 um 08:45 schrieb Riccardo (Jack) Lucchetti:
On Wed, 18 Sep 2019, Allin Cottrell wrote:
> Any other thoughts on this?
I personally don't like very much having both msample() and resample(),
but I can't see an easy way out of this.
I agree about the function proliferation point.
The second argument to resample must currently be >= 2, though. Perhaps
a value of 0 could be interpreted as "don't replace". If this
order-preserving thing is really needed, this then could become a third
choice argument?
cheers
sven
9:18 a.m.
On Thu, 19 Sep 2019, Sven Schreiber wrote:
If this order-preserving thing is really needed, this then could
become
a third choice argument?
Actually, I'm struggling to find a case when the order-preserved output
would be indispensable. Allin, do you have anything in mind?
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
2:37 p.m.
On Thu, 19 Sep 2019, Riccardo (Jack) Lucchetti wrote:
On Thu, 19 Sep 2019, Sven Schreiber wrote:
> If this order-preserving thing is really needed, this then could become a
> third choice argument?
Actually, I'm struggling to find a case when the order-preserved output would
be indispensable. Allin, do you have anything in mind?
I wasn't envisaging preservation of order as indispensable, just
thinking of uses where order didn't matter -- in which case I thought
it was potentially a good deal quicker to leave the order as is
(particularly if the sample size is a large fraction of the available
rows). But I now think I have a better way of combining selection and
shuffling, so I'm ready to drop the idea of an order-preserving
option.
As for the possibility of folding the without-replacement variant into
resample(), is the idea that without-replacement would also support
the blocksize option? (It would be easier if it didn't!)
Allin
2:55 p.m.
On Thu, 19 Sep 2019, Allin Cottrell wrote:
As for the possibility of folding the without-replacement variant
into
resample(), is the idea that without-replacement would also support
the blocksize option? (It would be easier if it didn't!)
Unfortunately, it may do (I think: but it'd be nice if Stefano Fachin, our
bootstrap guru, would give his opinion on the matter).
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
3:03 p.m.
Am 19.09.2019 um 16:55 schrieb Riccardo (Jack) Lucchetti:
On Thu, 19 Sep 2019, Allin Cottrell wrote:
> As for the possibility of folding the without-replacement variant into
> resample(), is the idea that without-replacement would also support
> the blocksize option? (It would be easier if it didn't!)
Unfortunately, it may do (I think: but it'd be nice if Stefano Fachin,
our bootstrap guru, would give his opinion on the matter).
His opinion and examples of use cases would certainly be welcome. But
OTOH I think you're right, it's probably good to generalize this
properly now if it is touched at all.
Actually, maybe the simple solution to introduce a new third argument
isn't so bad even if the middle argument isn't always used. One could
write something like
Xnew = resample(X, , 0)
to switch off the replacement (which would remain the default) and
simply leave the middle argument at its default (no blocks).
cheers
sven
3:42 p.m.
On Thu, 19 Sep 2019, Sven Schreiber wrote:
Am 19.09.2019 um 16:55 schrieb Riccardo (Jack) Lucchetti:
> On Thu, 19 Sep 2019, Allin Cottrell wrote:
>
> > As for the possibility of folding the without-replacement variant into
> > resample(), is the idea that without-replacement would also support
> > the blocksize option? (It would be easier if it didn't!)
>
> Unfortunately, it may do (I think: but it'd be nice if Stefano Fachin,
> our bootstrap guru, would give his opinion on the matter).
His opinion and examples of use cases would certainly be welcome. But
OTOH I think you're right, it's probably good to generalize this
properly now if it is touched at all.
On reflection it seems to me that sampling by blocks is conceptually
quite tricky in the without-replacement case. In a toy example,
suppose you have 10 observations and a blocksize of 3. The first
random selection takes, say, a block starting at t = 3. Now we can't
start a subsequent block at t = 1 or 2 because there won't be enough
contiguous observations. Each draw takes blocksize + p observations
out of the pool, where p is a random variable, so the max number of
draws we can make also becomes a random variable.
Or in this case would we take "without replacement" as applying to
blocks and not individual observations? (That is, it's OK to reuse
observations so long as all the blocks are unique.)
Allin
3:17 p.m.
On Thu, 19 Sep 2019, Riccardo (Jack) Lucchetti wrote:
On Thu, 19 Sep 2019, Allin Cottrell wrote:
> As for the possibility of folding the without-replacement variant into
> resample(), is the idea that without-replacement would also support
> the blocksize option? (It would be easier if it didn't!)
Unfortunately, it may do (I think: but it'd be nice if Stefano Fachin, our
bootstrap guru, would give his opinion on the matter).
OK, we'll wait for that. In the meantime there's a revised msample()
in git, without the order-preserving option and without a blocksize
option -- for testing if anyone wants to try it.
Allin
7:59 a.m.
Am 17.09.19 um 09:47 schrieb Riccardo (Jack) Lucchetti:
On Tue, 17 Sep 2019, Ioannis A. Venetis wrote:
> I agree Artur. I was recently in need for a resample without replacement
> function.
>
> I have not check your code but I had a look on what I did and it could
> be a refinement (following also Sven's answer).
>
> Say a matrix A is nxm.
>
> set k<=m.
>
> rvec = ranking( mrandgen(u,0,1,m,1) )'
>
> A=[,rvec[1:k]]
>
> (keeps k randomly selectled columns of A)
>
> I would like to see such a function in extra.gfn
>
I like it. We could have something like this in extra.gfn:
<hansl>
function matrix sample_wo_rep(matrix X, scalar k)
scalar n = rows(X)
if k > n
funcerr "what?"
endif
matrix sel = ranking(muniform(n,1))
return X[sel[1:k],]
end function
### example
X = mnormal(10, 3)
print X
eval sample_wo_rep(X, 6)
</hansl>
Interesting to see that this hansl-based version almost as fast as
compared to the native msample() function:
<hansl>
clear
set verbose off
function matrix sample_wo_rep (const matrix X,
const scalar k)
scalar n = rows(X)
if k > n
funcerr "what?"
endif
matrix sel = ranking(muniform(n,1))
return X[sel[1:k],]
end function
# Parameters
scalar R = 100
scalar C = 100
scalar SIZE = 5
if SIZE>R
stop
endif
scalar RUNS = 10000
matrix X = mnormal(R, C)
matrix runtime = zeros(RUNS, 2)
loop algo = 1..2 -q
loop i=1..RUNS -q
set stopwatch
if algo
matrix res = sample_wo_rep(X, SIZE)
else
matrix res = msample(X, SIZE)
endif
runtime[i, algo] = $stopwatch
endloop
endloop
printf "Mean runtime\n%12.9f\n", meanc(runtime)
printf "Ratio of Jack-to-msample() = %.3f\n",
meanc(runtime[,1])/meanc(runtime[,2])
</hansl>
Artur
9:20 a.m.
On Fri, 20 Sep 2019, Artur Tarassow wrote:
Interesting to see that this hansl-based version almost as fast as
compared
to the native msample() function:
Hm, not on my machine here. I get more or less the same time for various
combinations of R, C and SIZE.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
5:14 p.m.
On Fri, 20 Sep 2019, Artur Tarassow wrote:
[commenting on a function written by Jack]
Interesting to see that this hansl-based version almost as fast as
compared to the native msample() function:
<hansl>
clear
set verbose off
function matrix sample_wo_rep (const matrix X,
const scalar k)
scalar n = rows(X)
if k > n
funcerr "what?"
endif
matrix sel = ranking(muniform(n,1))
return X[sel[1:k],]
end function
# Parameters
scalar R = 100
scalar C = 100
scalar SIZE = 5
if SIZE>R
stop
endif
scalar RUNS = 10000
matrix X = mnormal(R, C)
matrix runtime = zeros(RUNS, 2)
loop algo = 1..2 -q
loop i=1..RUNS -q
set stopwatch
if algo
matrix res = sample_wo_rep(X, SIZE)
else
matrix res = msample(X, SIZE)
endif
runtime[i, algo] = $stopwatch
endloop
endloop
printf "Mean runtime\n%12.9f\n", meanc(runtime)
printf "Ratio of Jack-to-msample() = %.3f\n",
meanc(runtime[,1])/meanc(runtime[,2])
Unfortunately this test is running sample_wo_rep against itself: algo
runs from 1 to 2 so "if algo" is always true. If we substitute the
condition "if algo == 1" we get quite different results. For example,
with Artur's parameter values:
Mean runtime
0.000027866 0.000000776
Ratio of Jack-to-msample() = 35.932
Allin
9:27 a.m.
Am 20.09.19 um 19:14 schrieb Allin Cottrell:
On Fri, 20 Sep 2019, Artur Tarassow wrote:
[commenting on a function written by Jack]
> Interesting to see that this hansl-based version almost as fast as
> compared to the native msample() function:
>
> <hansl>
> clear
> set verbose off
>
> function matrix sample_wo_rep (const matrix X,
> const scalar k)
> scalar n = rows(X)
> if k > n
> funcerr "what?"
> endif
> matrix sel = ranking(muniform(n,1))
> return X[sel[1:k],]
> end function
>
> # Parameters
> scalar R = 100
> scalar C = 100
> scalar SIZE = 5
> if SIZE>R
> stop
> endif
> scalar RUNS = 10000
>
> matrix X = mnormal(R, C)
>
> matrix runtime = zeros(RUNS, 2)
>
> loop algo = 1..2 -q
> loop i=1..RUNS -q
> set stopwatch
> if algo
> matrix res = sample_wo_rep(X, SIZE)
> else
> matrix res = msample(X, SIZE)
> endif
> runtime[i, algo] = $stopwatch
> endloop
> endloop
>
> printf "Mean runtime\n%12.9f\n", meanc(runtime)
> printf "Ratio of Jack-to-msample() = %.3f\n",
> meanc(runtime[,1])/meanc(runtime[,2])
Unfortunately this test is running sample_wo_rep against itself: algo
runs from 1 to 2 so "if algo" is always true. If we substitute the
condition "if algo == 1" we get quite different results. For example,
with Artur's parameter values:
Mean runtime
0.000027866 0.000000776
Ratio of Jack-to-msample() = 35.932
Oh, good point. I've thought that <if algo> is equivalent to <if
algo==1>. I also obtain a factor of about 26. As expected the C based
version dominates.
Thanks,
Artur
4:13 p.m.
On Sat, 21 Sep 2019, Artur Tarassow wrote:
Am 20.09.19 um 19:14 schrieb Allin Cottrell:
[revised timing
comparison]
>
> Ratio of Jack-to-msample() = 35.932
Oh, good point. I've thought that <if algo> is equivalent to <if
algo==1>. I also obtain a factor of about 26. As expected the C
based version dominates.
Since we've collectively devoted a fair amount of time to this issue
(the sampling without replacement, not just the timing) I think we
should try to resolve it as expeditiously as possible.
Sorry, this is going to be a bit long but I'll try to be concise.
1) Choice of algorithm. Jack's suggestion is very clever (no
surprise!) but it does appear to be a good deal slower than the
algorithm used by msample(). It's not just hansl-versus-C. I tried
coding Jack's ranking-based algorithm in C as an alternative
back-end for msample(), and in that case my timings show the hansl
version taking only about 20 percent longer than C.
The msample algorithm is quite different:
* Create a "pool" array of integers, {0, 1, ..., r-1}, where r is
the number of rows in the source matrix, and set poolsize = r.
Set k = 0.
* For i = 0 to n-1, where n is the required number of draws,
- Generate a random unsigned integer, u, on [0,poolsize)
- Write row pool[u] of the source matrix into row k of the
matrix to be returned
- poolsize <- poolsize - 1
- k <- k + 1
- If u did not select the last member of the pool, shift
the "tail" of the pool (to the right of the selection) one
place to the left to close up the gap left by the selected
row. Use the C-library function memmove().
This emulates the actual drawing of balls from an urn without
replacement.
I was not sure how expensive the calls to memmove() would be, but in
fact this method seems to be quite fast. Comparing C implementation
with C implementation, it's at least an order of magnitude faster
than using ranking, and very much faster if you're drawing
relatively few rows from a matrix with many rows.
In addition, in running lots of cases for timing purposes, the
problem that concerned Sven with the ranking method -- namely the
possibility of getting ties -- showed itself a couple of times:
Invalid selection vector
*** error in function sample_wo_rep
X[sel[1:n],]
(So there must have been a *.5 somewhere.)
I therefore propose that we use the "balls from urn" algorithm.
2) Interface. It would be nice if we could generalize an existing
function, resample, rather than add another. But I really don't see
how that could work.
resample takes a matrix argument plus an optional block-length
argument. Its return value always has as many rows as the input.
msample takes a matrix argument plus a number of cases to be drawn,
which seems the natural interface for such a function. (And I now
take it that we don't want to try adding a block-length parameter
when sampling without replacement.)
Trying to meld these two seems to me artificial. The documentation
would have to be kinda convoluted, such that a user might well
think, "Look, you've got two different functions in here, why not
just give them different names?"
The name "msample" is not set in stone; if anyone has a better
suggestion, go for it. But as of now I think we need a new function.
Allin
6:04 p.m.
Am 21.09.2019 um 18:13 schrieb Allin Cottrell:
On Sat, 21 Sep 2019, Artur Tarassow wrote:
Since we've collectively devoted a fair amount of time to this
issue
(the sampling without replacement, not just the timing) I think we
should try to resolve it as expeditiously as possible.
Sorry, this is going to be a bit long but I'll try to be concise.
Thanks, Allin. What I still don't really fully understand is the
features that are needed (or are likely to be needed). Initially I said
we don't want ordered draws, but I'm not so sure that that couldn't also
be a sensible option. Something similar goes for the block length.
Maybe Artur's reference to R's sample function is a hint that we should
study more closely what that function does?
cheers
sven
6:48 p.m.
On Sat, 21 Sep 2019, Sven Schreiber wrote:
Am 21.09.2019 um 18:13 schrieb Allin Cottrell:
> On Sat, 21 Sep 2019, Artur Tarassow wrote:
> Since we've collectively devoted a fair amount of time to this issue
> (the sampling without replacement, not just the timing) I think we
> should try to resolve it as expeditiously as possible.
>
> Sorry, this is going to be a bit long but I'll try to be concise.
Thanks, Allin. What I still don't really fully understand is the
features that are needed (or are likely to be needed). Initially I said
we don't want ordered draws, but I'm not so sure that that couldn't also
be a sensible option.
I suppose it could be -- though my original preference for allowing
it was really just based on (what I took to be) its relative
computational simplicity, something I no longer think is very
important now that I've figured out an efficient algorithm for doing
the sampling + scrambling in one pass.
Something similar goes for the block length.
I think sampling by blocks without replacement is really problematic
and probably not to be messed with. Do you have a response to my
diagnosis in
https://www.mail-archive.com/gretl-users@gretlml.univpm.it/msg14200.html
The second possible approach I mention there would be
straightforward to implement but I seriously doubt it would have any
desirable statistical properties.
Maybe Artur's reference to R's sample function is a hint that
we should
study more closely what that function does?
The signature is
sample(x, size, replace = FALSE, prob = NULL)
where "size" is the number of draws and "prob" is an optional vector
of weights. As Artur said, sampling by blocks is not supported.
There's one more optional boolean argument, useHash, which is
available only if replace = FALSE, prob = NULL, and size <= n/2.
This indicates that "the hash-version of the algorithm should be
used" and is recommended for large n. Oddly, I don't see any
indication of whether useHash is the default when its conditions of
application are met.
Allin
7:01 p.m.
On Sat, 21 Sep 2019, Allin Cottrell wrote:
On Sat, 21 Sep 2019, Sven Schreiber wrote:
> What I still don't really fully understand is the features that
> are needed (or are likely to be needed).
Sorry for the bombardment, but one more observation.
Artur has told us that msample meets his needs (sampling by block
not required); Ioannis has given us a quite detailed account of his
needs and my understanding is that msample would meet those too; and
Stefan has made it plain that bootstrapping (where you might well
want sampling by block) has no use for sampling without replacement.
I realize these instances don't necessarily cover all possible uses
-- and I agree that it's worth thinking ahead. But at this point I
see a fairly strong presumption that msample would be sufficient.
Allin
7:42 p.m.
Am 21.09.2019 um 21:01 schrieb Allin Cottrell:
On Sat, 21 Sep 2019, Allin Cottrell wrote:
> On Sat, 21 Sep 2019, Sven Schreiber wrote:
>
>> What I still don't really fully understand is the features that are
>> needed (or are likely to be needed).
Sorry for the bombardment, but one more observation.
Artur has told us that msample meets his needs (sampling by block not
required); Ioannis has given us a quite detailed account of his needs
and my understanding is that msample would meet those too; and Stefan
has made it plain that bootstrapping (where you might well want sampling
by block) has no use for sampling without replacement.
I realize these instances don't necessarily cover all possible uses --
and I agree that it's worth thinking ahead. But at this point I see a
fairly strong presumption that msample would be sufficient.
I really would very much prefer if somebody more knowledgeable in this
area than myself had to take the counter-position. I'm also one of those
guys who stick to the bootstrap 99% of the times, and have much less
exposure to other resampling techniques.
Having said that, for example in section 4.1 of this article there seems
to be subsampling with a block size b: Linton/Maasoumi/Whang (2003),
https://pdfs.semanticscholar.org/9378/bb9e0ae9c7fbf669132ec07f77931a00b1b...
Or Politis & Romano (1994, Annals) see section 3 (stationary time
series...), where they explicitly seem to consider subsampling for
"blocks of size b of the consecutive data...".
(https://projecteuclid.org/euclid.aos/1176325770, open access pdf)
The general idea of the "subsampling for dependent data" method seems to
be: Since we have dependent data, we need to draw blocks. This idea is
similar to a block bootstrap. As we do subsampling, in principle we
would look at *all* possible subsamples (non-overlapping, no
replacement). However, the combinatorics may become too heavy --
implicitly here a lot of data are assumed. Thus instead of looking at
everything, in practice a bunch of random subsamples are drawn (without
replacement). This method works as an approximation given the
established theoretical convergence results.
So if somebody wants to do that, she needs random draws of contiguous
blocks without replacement. Since the aim is not to come close to
exhausting the available data (lots of them by assumption), your
argument (that the number of possible drawn blocks would be random) is
correct but irrelevant.
OK, sorry for any mistakes in this proclamation, I'm just an amateur
browsing the literature there. Given that I'm not using these methods
myself, I can also happily live with gretl functions that do not cover
everything in this area, such as the current msample. But I think in
principle those features would not be absurd.
FWIW,
sven
11:52 p.m.
On Sat, 21 Sep 2019, Sven Schreiber wrote:
[F]or example in section 4.1 of this article there seems to be
subsampling with a block size b: Linton/Maasoumi/Whang (2003),
https://pdfs.semanticscholar.org/9378/bb9e0ae9c7fbf669132ec07f77931a00b1b...
Or Politis & Romano (1994, Annals) see section 3 (stationary time
series...), where they explicitly seem to consider subsampling for
"blocks of size b of the consecutive data...".
(https://projecteuclid.org/euclid.aos/1176325770, open access pdf)
Thanks for the references.
I found Linton/Maasoumi/Whang heavy going but I think I understand
Politis & Romano OK (and since LMW are drawing on the latter maybe I
understand them well enough for current purposes too).
If so, it appears they're not doing random sampling, they're
considering all contiguous sub-samples of a certain size. This is
clear in section 1 of Politis & Romano where they're talking about
"subsets of the data that retain the dependence structure of the
data." They say:
<quote>
For example, if X_i,...,X_n represent n observations from a
stationary time series, the statistic is recomputed over the n - b +
1 subsets of size b of the form {X_i, X_{i+1},...X_{i+b-1}}.
</quote>
Let n = 6 and b = 2, so n - b + 1 = 5. It's clear what 5 subsets
they must be talking about:
(1,2) (2,3) (3,4) (4,5) (5,6)
There's nothing random here and it would be easy to generate such
subsets in hansl without the help of any special function.
Randomness might come in at a "meta" level if one wanted to select a
random subset of the full set of such samples -- presumably, without
replacement. That could be done easily using msample() on a vector
of contiguous index values (just {1,2,3,4,5} for the toy case
above). Given the selected indices, getting the samples would then
be a simple deterministic calculation.
Allin
8:55 a.m.
Am 22.09.2019 um 01:52 schrieb Allin Cottrell:
If so, it appears they're not doing random sampling, they're
considering
all contiguous sub-samples of a certain size.
...
There's nothing random here and it would be easy to generate such
subsets in hansl without the help of any special function.
I agree on the non-randomness with respect to all possible subsets.
Randomness might come in at a "meta" level if one wanted to select a
random subset of the full set of such samples -- presumably, without
replacement.
Yes, that's how I understood the randomness coming in.
That could be done easily using msample() on a vector of
contiguous index values (just {1,2,3,4,5} for the toy case above).
Given
the selected indices, getting the samples would then be a simple
deterministic calculation.
Not sure I understand that. Previously you said that taking blocks
without replacement would be tricky. Now you say it's simple?
(But I haven't had time right now to give it much thought, sorry.)
cheers
sven
6:05 p.m.
On Sun, 22 Sep 2019, Sven Schreiber wrote:
Am 22.09.2019 um 01:52 schrieb Allin Cottrell:
> If so, it appears they're not doing random sampling, they're considering
> all contiguous sub-samples of a certain size.
...
>
> There's nothing random here and it would be easy to generate such
> subsets in hansl without the help of any special function.
I agree on the non-randomness with respect to all possible subsets.
> Randomness might come in at a "meta" level if one wanted to select a
> random subset of the full set of such samples -- presumably, without
> replacement.
Yes, that's how I understood the randomness coming in.
That could be done easily using msample() on a vector of
> contiguous index values (just {1,2,3,4,5} for the toy case above). Given
> the selected indices, getting the samples would then be a simple
> deterministic calculation.
Not sure I understand that. Previously you said that taking blocks
without replacement would be tricky. Now you say it's simple?
In the Politis and Romano case under consideration, they're doing an
analysis that involves calculating a certain statistic m times, for
each of m contiguous subsamples of length b from the full data.
They're fine with these subsamples overlapping, and they're not
sticking them together into a combined sample of size m*b. As they
say, there will always be n - b + 1 such subsamples available. My
observation was just that if m < n subsamples were judged "enough" for
the job, one could easily select m of them using msample() on their
indices (the index for each subsample just being the observation at
which it starts). I mean (trivial example, for a matrix X with n
rows):
sample_starts = seq(1,n-b+1)
samples_used = msample(sample_starts, m)
matrices samples = array(m)
loop j=1..m
istart = samples_used[i]
istop = istart + b - 1
samples[i] = X[istart:istop,]
endloop
This is quite different from the case I couldn't get my head around,
namely drawing a single sample of size m > b with block length b and
no replacement of observations (and so no overlap of blocks).
Allin
12:17 p.m.
Am 22.09.2019 um 20:05 schrieb Allin Cottrell:
On Sun, 22 Sep 2019, Sven Schreiber wrote:
> Am 22.09.2019 um 01:52 schrieb Allin Cottrell:
>> Randomness might come in at a "meta" level if one wanted to select a
>> random subset of the full set of such samples -- presumably, without
>> replacement.
> Yes, that's how I understood the randomness coming in.
>
> That could be done easily using msample() on a vector of
>> contiguous index values (just {1,2,3,4,5} for the toy case above). Given
>> the selected indices, getting the samples would then be a simple
>> deterministic calculation.
In the Politis and Romano case under consideration, they're doing an
analysis that involves calculating a certain statistic m times, for
each of m contiguous subsamples of length b from the full data.
They're fine with these subsamples overlapping,
I think you're right, and
the overlap is a crucial point.
...
This is quite different from the case I couldn't get my head
around,
namely drawing a single sample of size m > b with block length b and
no replacement of observations (and so no overlap of blocks).
Yes, I give up ;-)
cheers
sven
>
12:28 p.m.
Am 23.09.2019 um 14:17 schrieb Sven Schreiber:
Yes, I give up ;-)
P.S.: I still think the name 'msample' is suboptimal. The meaning of the
initial m is not clear to me. If it refers to "matrix", I'd say it
should also work on a series. So perhaps "subsample", or "pick", or
whatever.
thanks
sven
12:29 p.m.
On Mon, 23 Sep 2019, Sven Schreiber wrote:
Am 23.09.2019 um 14:17 schrieb Sven Schreiber:
>
> Yes, I give up ;-)
>
P.S.: I still think the name 'msample' is suboptimal. The meaning of the
initial m is not clear to me. If it refers to "matrix", I'd say it
should also work on a series. So perhaps "subsample", or "pick", or
whatever.
"pick" is too generic, I think
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
7:18 p.m.
On Mon, 23 Sep 2019, Sven Schreiber wrote:
Am 23.09.2019 um 14:17 schrieb Sven Schreiber:
>
> Yes, I give up ;-)
>
P.S.: I still think the name 'msample' is suboptimal. The meaning of the
initial m is not clear to me. If it refers to "matrix", I'd say it
should also work on a series. So perhaps "subsample", or "pick", or
whatever.
Yes, the 'm' is supposed to suggest "matrix" -- but actually I think
that's right.
We already have "smpl n --random" for datasets, and allowing a
single series argument to <function whose name is under discussion>
would be a little strange. You couldn't get a series back, it would
have to be a vector -- unless you set n equal to the number of
observations, thereby just scrambling the series rather than
subsampling, which seems an unlikely project.
If for some reason you want to draw a subsample from a particular
series in isolation from the dataset as a whole, it seems OK to me
that you should "cast" it to a vector -- as in, for a series x,
matrix subx = msample({x}, n)
And in the unlikely event that you want to put into the dataset a
scrambled version of an entire series, you could do
series scrambled_x = msample({x}, $nobs)
("casting" the vector result back to a series).
Admittedly there's no rule saying that a function taking a series
argument can't return a vector but this would be unusual. Most such
functions return either a transformed series or a scalar statistic
(with the exception of values()).
Allin
7:36 a.m.
Am 23.09.2019 um 21:18 schrieb Allin Cottrell:
On Mon, 23 Sep 2019, Sven Schreiber wrote:
> Am 23.09.2019 um 14:17 schrieb Sven Schreiber:
>>
>> Yes, I give up ;-)
>>
> P.S.: I still think the name 'msample' is suboptimal. The meaning of the
> initial m is not clear to me. If it refers to "matrix", I'd say it
> should also work on a series. So perhaps "subsample", or "pick",
or
> whatever.
Yes, the 'm' is supposed to suggest "matrix" -- but actually I think
that's right.
We already have "smpl n --random" for datasets, and allowing a single
series argument to <function whose name is under discussion> would be a
little strange. You couldn't get a series back, it would have to be a
vector
OK, but now that we have established we don't need a blocksize for the
no-replacement case, why not go back to 'resample'?
1st arg: as before
2nd arg: 0 to indicate no replacement
3rd arg: (new and optional) sample size
Having the 2nd arg != 0 and the 3rd arg present would be an error. Of
course, using the 3rd arg might be considered sufficient and the 2nd arg
could be left empty (null) then.
cheers
sven
1915
days inactive
1922
days old
34 comments
5 participants
participants (5)
-
Allin Cottrell -
Artur Tarassow -
Ioannis A. Venetis -
Riccardo (Jack) Lucchetti -
Sven Schreiber