inprob.C

void tille_incprob(double a[], int *pn, int *plen){

     double a_sum=0;
     int i,n,len,l,l1;

     n=*pn;
     len=*plen;
     for(i=0;i<len; i++){
       a_sum+=a[i];
     }
     l=0;
     for(i=0; i<len; i++){
       a[i]=a[i]*n/a_sum;
       if (a[i]>=1) l++;
     }
     if(l>0){
       l1=0;
       while(l!=l1){
	 a_sum=0;
	 for(i=0;i<len; i++){
	   if (a[i]>0 && a[i]<1)
	     a_sum+=a[i];
	 }
	 for(i=0;i<len; i++){
	   if (a[i]>0){
	     if (a[i]<1) {
	       a[i]=(n-l)*a[i]/a_sum;
	     } else {
	       a[i]=1;
	     }
	   }
	 }
	 l1=l;
	 l=0;
	 for(i=0; i<len; i++){
	   if (a[i]>=1) l++;
	   }
       }
     }
     return;
}
     

inprob.R

tille_incprob<-function(p,n){
    n<-as.integer(round(n))
    len<-as.integer(length(p))
    storage.mode(p)<-"numeric"
    .C("tille_incprob", p, n, len)[[1]]
}



sample_tille <-function (pik, eps = 1e-06) 
{
    if (any(is.na(pik))) 
        stop("there are missing values in the pik vector")
    n = sum(pik)
    n = .as_int(n)
    list = pik > eps & pik < 1 - eps
    pikb = pik[list]
    N = length(pikb)
    s = pik
    if (N < 1) 
        stop("the pik vector has all elements outside of the range [eps,1-eps]")
    else {
        n = sum(pikb)
        sb = rep(1, N)
        b = rep(1, N)
        for (i in 1:(N - n)) {
            a = tille_incprob(pikb, N - i)
            v = 1 - a/b
            b = a
            p = v * sb
            p = cumsum(p)
            u = runif(1)
            for (j in 1:length(p)) if (u < p[j]) 
                break
            sb[j] = 0
        }
        s[list] = sb
    }
    s
}