Difference between calculated inclusion probability and what is returned by sampling function?

I have a (small) population from which I wish to sample. I assign probabilities proportional to $y$. I enumerate the possible samples and then determine the probability of each sample occurring based on the product of the probabilities for each $y_i$ in the sample. I add up the probabilities for the samples that contain the $y_1$ and I believe (incorrectly?) that under the assumption of independence (i.e. with replacement sampling) this gives me the inclusion probability for $y_1$. I look at the inclusion probabilities returned by the inclusionprobabilities function in the sampling package and I get a different answer. I do not understand why, is someone able to explain?

library(survey)
library(sampling)
library(gtools)

set.seed(123)
y <- c(1190,26751,68570,34536)
p <- y/sum(y)

df <- data.frame(permutations(n=length(y), r=2, v=1:length(y), repeats.allowed = T))
df$p <- p[df$X1] * p[df$X2]; df # X1 and X2 denote the index of the y value that is included # in the sample. X1 X2 p 1 1 1 0.00008245932 2 1 2 0.00185367169 3 1 3 0.00475145854 4 1 4 0.00239312195 5 2 1 0.00185367169 6 2 2 0.04167022794 7 2 3 0.10681198947 8 2 4 0.05379697926 9 3 1 0.00475145854 10 3 2 0.10681198947 11 3 3 0.27378782541 12 3 4 0.13789611111 13 4 1 0.00239312195 14 4 2 0.05379697926 15 4 3 0.13789611111 16 4 4 0.06945282329 samplesSet <- data.frame(df[1 == df$X1  | 1 == df$X2, ]) sum(samplesSet$p)

pik <- inclusionprobabilities(y, 2)
data.frame(pik=pik,name=1:length(y))

Update: Thanks both @whuber and @StasK. It is clear that the inclusion probabilities reflect sampling without replacement. However, I am uncertain what the inclusion probabilities returned by inclusionprobabilities are. They seem to be calculated as:

$$n \frac{y_i}{\sum_{i=1}^{N} y_i}$$

and have an adjustment to ensure that no probability is greater than 1 and also that the sum of the probabilities corresponds to the sample size.

If I assume that my population is $y=\{1,2,3\}$ such that the probabilities of selection are $\frac{1}{6}$, $\frac{2}{6}$ and $\frac{3}{6}$ and then I take a sample of 2, I calculate the inclusion probabilities to be $\frac{5}{12}$, $\frac{11}{15}$ and $\frac{17}{20}$ respectively. Clearly, these are not what is returned by inclusionprobabilities and so my question now is have I calculated the inclusion probabilities incorrectly or is the inclusionprobabilities function returning something that represents the inclusion probabilities but isn't actually the inclusion probabilities?

myn <- 2
a <- c(1,2,3)
p <- myn * a/sum(a); p
 0.3333333 0.6666667 1.0000000

inclusionprobabilities(a, myn)
 0.3333333 0.6666667 1.0000000

Thanks.

• The help is indeed abysmal. Protect yourself by testing this function on simple arguments with known answers. For instance, inclusionprobabilities(1:2,2) returns the vector 1 1. What does that tell you about the assumed form of sampling? Could this possibly reflect sampling with replacement? (Such prophylactic testing is essential when learning to use any package--the ultimate arbiter of questions like this is what the computer does, not what the help pages seem to say!)
– whuber
Jan 23 '15 at 16:53
• Sampling with unequal probabilities is really weird, and it does not always give you the answer you expect, although that mostly has to do for sampling with replacement. Selection probabilities, and especially the pairwise selection probabilities, depend on the particular sampling algorithm, see Hanif & Brewer (1983) -- nearly impossible to find -- and Tille 2006. If inclusionprobabilities() indeed refer to sampling WOR, you need to filter on X1!=X2. Jan 24 '15 at 4:31
• Thank you @whuber. The return of 1 1 tells me that it does not reflect sampling with replacement. I have done some additional tests and updated my original question. It is still not clear to me what the inclusionprobabilities() function is returning. Jan 27 '15 at 4:59

print(sum(df$p[df$X1==k|df\$X2==k]))