# HG changeset patch
# User doublevvinged<doublevvinged@abv.bg>
# Date 1435403457 -10800
#      Sat Jun 27 14:10:57 2015 +0300
# Node ID 9112fdc69be6243a726efe22c53fa794c9600414
# Parent  b15b513a726577a806045e540659a23d3dab5af6
Function randsample can now take as input a weight vector when sampling with replacement.

diff -r b15b513a7265 -r 9112fdc69be6 inst/randsample.m
--- a/inst/randsample.m	Sat Jun 06 14:45:47 2015 +0100
+++ b/inst/randsample.m	Sat Jun 27 14:10:57 2015 +0300
@@ -1,4 +1,4 @@
-## Copyright (C) 2014 - Nir Krakauer
+## Copyright (C) 2014 - doublevvinged
 ##
 ## This progrm is free software; you can redistribute it and/or modify
 ## it under the terms of the GNU General Public License as published by
@@ -13,83 +13,221 @@
 ## You should have received a copy of the GNU General Public License
 ## along with this program. If not, see <http://www.gnu.org/licenses/>.
 
-## Author: Nir Krakauer <nkrakauer@ccny.cuny.com>
-
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{y} =} randsample (@var{v}, @var{k}, @var{replacement}=false [, @var{w}])
-## Elements sampled from a vector.
+## @deftypefn {Function File} {@var{y} =} randsample (@var{v}, @var{n}, @var{replacement}, @var{weights}  )
+## This function returns a vector @var{y}[n x 1] of size @var{n} by 1 elements
+## of random elements sampled from a vector @var{v}.
+## If @var{replacement}='true', replacement is enabled and
+## multiple occurances of the same element are allowed 
+## If in addition to that the non negative weights @var{weights}[n x 1] >0 are given,
+## then each sampel is taken with probability proportional to its weight
+## The weights are not necessary normalised.
+## The default @var{replacement} method is without replacement.
 ##
-## Returns @var{k} random elements from a vector @var{v} with @var{n} elements, sampled without or with @var{replacement}.
+## Example:
 ##
-## If @var{v} is a scalar, samples from 1:@var{v}.
+##p=[0.1 0.2 0.6];
 ##
-## If sampling with replacement, can specify a weight vector @var{w} of the same size as @var{v} such that the probablility of each element being sampled is proportional to @var{w}.
+##for i=1:1000
 ##
-## Randomization is performed using rand().
+## res(1:2,i)=randsample([2,4,7],2,'true', p );
+##
+##end
+##
+##res=res(:);
+##
+##hist(res)
 ##
 ## @seealso{randperm}
 ## @end deftypefn
+## Author: doublevvinged
 
-function y = randsample(v,k,replacement=false,w=[])
 
-  if isscalar (v)
-    n = v;
-    vector_v = false;
-  elseif isvector (v)
-    n = numel (v);
-    vector_v = true;
-  else
-    error ('Octave:invalid-input-arg', 'The input v must be a vector or positive integer.');
-  endif
+
+
+function samp=randsample(x,n,repl,p)
+%# by doublevvinged
+%# x is the sample vector
+%# n is the number of requested samples
+%# repl if sampling is with replacement
+%
+%# build in function closest_match
+%# this is a binary search that returns the closest lower match of a in sorted increasing x
+function yes=closest_match(a,x)
+yes=0;
+lo=1;
+up=length(x);
+i=ceil( (up+lo)/2 );
+k=up-lo;
+while k>1
+    if x(i)==a
+        yes=i;
+        break;
+    else 
+        if x(i)>=a
+            up=i;
+        else            
+            lo=i;           
+        end
+    end
+i=ceil( (up+lo)/2 );
+k=up-lo;
+
+end
+
+if up==length(x)
+   if(x(up)<a) 
+   yes=up;
+   else
+   yes=lo;
+   end
+else
+  yes=lo;
+end
    
-  if k < 0
-    error ('Octave:invalid-input-arg', 'The input k must be a positive integer.');
-  endif
+ 
+endfunction
 
 
-  if replacement #sample with replacement
-    if isempty (w) #all elements are equally likely to be sampled
-      y = round (n * rand(1, k) + 0.5);
-    else
-      w = w / sum(w);
-      w = [0 cumsum(w(:))'];
-      y = arrayfun(@(x) find(w <= x, 1, "last"), rand (1, k)); #distribute k uniform random deviates based on the given weighting
-    endif
-  else #sample without replacement
-    y = randperm (n, k);
-  endif
-  
-  if vector_v
-    y = v(y);
-  endif
-  
-endfunction 
- 
 
-%!test
-%! n = 20;
-%! k = 5;
-%! x = randsample(n, k);
-%! assert (size(x), [1 k]);
-%! x = randsample(n, k, true);
-%! assert (size(x), [1 k]);
-%! x = randsample(n, k, false);
-%! assert (size(x), [1 k]);
-%! x = randsample(n, k, true, ones(n, 1));
-%! assert (size(x), [1 k]);
-%! x = randsample(1:n, k);
-%! assert (size(x), [1 k]);
-%! x = randsample(1:n, k, true);
-%! assert (size(x), [1 k]);
-%! x = randsample(1:n, k, false);
-%! assert (size(x), [1 k]);
-%! x = randsample(1:n, k, true, ones(n, 1));
-%! assert (size(x), [1 k]);
-%! x = randsample((1:n)', k);
-%! assert (size(x), [k 1]);
-%! x = randsample((1:n)', k, true);
-%! assert (size(x), [k 1]);
-%! x = randsample((1:n)', k, false);
-%! assert (size(x), [k 1]);
-%! x = randsample((1:n)', k, true, ones(n, 1));
-%! assert (size(x), [k 1]);
+if exist('repl')
+  if strcmp(repl,'true')
+      if exist('p')
+%          # generate the look up table for the specified distribution p
+          %look_up(k,:)=[ lower interval of k-th element of x, upper interval of k-th element of x]
+          look_up=zeros(length(p),2);
+          look_up(1,:)=[0,p(1)];
+          for i=2:length(p)
+          look_up(i,1)=look_up(i-1,2);
+          look_up(i,2)=look_up(i,1)+p(i);
+          end
+
+          for j=1:n
+%          # a random index sample of elements of x from a distribution p
+          %ind= closest_match( unifrnd(0,look_up(size(look_up,1),2 )  ), look_up(:,1)  ) 
+          % Note that the probability distribution does not need to be normalised
+          samp(j)= x(  closest_match( unifrnd(0,look_up(size(look_up,1),2 )  ), look_up(:,1)  )  );
+          end
+      else %# exist('p') no sample distribution specified sample from uniform
+          look_up=zeros(length(x),2);
+          look_up(1,:)=[0,1/length(x)];
+
+            for i=2:length(x)
+            look_up(i,1)=look_up(i-1,2);
+            look_up(i,2)=look_up(i,1)+1/length(x);
+            end
+
+            for j=1:n
+            %ind= closest_match( unifrnd(0,look_up(size(look_up,1),2 )  ), look_up(:,1)  ) 
+            % Note that the probability distribution does not need to be normalised
+            samp(j)= x(  closest_match( unifrnd(0,look_up(size(look_up,1),2 )  ), look_up(:,1)  )  );
+            end            
+      end  # exist('p')
+
+   else % strcmp(repl,'true')
+        % NOTE THAT IF SAMPLING IS repl IS NOT 'true' SAMPLING IS WITHOUT REPLACEMENT
+        % sample WITHOUT replacement
+          if n>=length(x)
+          samp=0;
+          display('Requested sample size without replacement is greater than the size of the input')
+          break
+          end  
+%          # this is reduced to a random permutation
+          sampl=randperm(length(x),n);
+%          # map the indexes to the input vector
+          for j=1:n
+          samp(j)= x(  sampl(j)  );
+          end
+%          samp=transpose(samp);
+  end % strcmp(repl,'true')
+    
+else % exist('repl') 
+        % Default is WITHOUT replacement
+        if n>=length(x)
+        samp=0;
+        display('Requested sample size without replacement is greater than the size of the input')
+        break
+        end  
+%        # this is reduced to a random permutation
+        sampl=randperm(length(x),n);
+%        # map the indexes to the input vector
+        for j=1:n
+        samp(j)= x(  sampl(j)  );
+        end
+%        samp=transpose(samp);
+end # exist('repl')
+
+
+
+
+end # function samp=randsamp_doublevvinged(x,n,repl,p)
+
+
+%% Testing
+
+%p=[0.1 0.1 0.2 0.6]
+%randsample(1:4,2, 'true',p )
+% randsample('abcd',2, 'true',p )
+
+%p=[0.1 0.1 0.2 0.6]
+%k=1000
+%res=zeros(2,k);
+%for i=1:k
+%res(1:2,i)=randsample(1:4,2, 'true',p );
+%end
+%res=res(:);
+%hist(res)
+%
+%
+%p=[0.1 0.2 0.6]
+%k=1000
+%res=zeros(2,k);
+%for i=1:k
+%res(1:2,i)=randsample([2,4,7],2,'true', p );
+%end
+%res=res(:);
+%hist(res)
+%
+%
+%k=1000
+%res=zeros(2,k);
+%for i=1:k
+%res(1:2,i)=randsample([2,4,7],2,'true');
+%end
+%res=res(:);
+%hist(res)
+%
+%
+%
+%k=1000
+%res=zeros(2,k);
+%for i=1:k
+%res(1:2,i)=randsample([2,4,7],2,'false');
+%end
+%res=res(:);
+%hist(res)
+%
+%
+%
+%
+%p=[0.1 0.2 0.6 0.3]
+%k=1000
+%res=zeros(3,k);
+%for i=1:k
+%res(:,i)=randsample([2,4,6,7],3,'false',p);
+%end
+%res=res(:);
+%hist(res)
+%
+%
+%
+%p=[0.1 0.2 0.6 0.3]
+%k=1000
+%res=zeros(3,k);
+%for i=1:k
+%res(:,i)=randsample([2,4,6,7],3,'false',p);
+%end
+%res=res(:);
+%hist(res)
+
+
