################################################################## # # # R Program for Calculating Second Order Inclusion Probabilities # # # # 1. For repeated simulations, it is preferred to calculate all # # pi_ij for i, j=1, ..., N and save the N x N pi_{ij} matrix # # # # (1) Input: N is the population size; n is the sample size # # p is the N x 1 vector of the size variable # # (2) Output: piij is the N x N matrix for pi_{ij} # # # # 2. The following program works well for N<=800 # # The program takes very long for N>=1500 # # # # Written by Changbao Wu, July 2004 # # # ################################################################## ## ## p is the normalized size vector (sum_{i=1}^N p_i=1) ## ################################################################## piij<-matrix(0,N,N) Lm<-rep(0,n) q<-rep(0,N) for(i in 1:N) q[i]<-p[i]/(1-n*p[i]) Lm[1]<-1 for(i in 2:n){ for(r in 1:(i-1)) Lm[i]<-Lm[i]+((-1)^(r-1))*sum(q^r)*Lm[i-r] Lm[i]<-Lm[i]/(i-1) } t1<-(n+1)-(1:n) Kn<-1/sum(t1*Lm/n^t1) Lm2<-rep(0,n-1) t2<-(1:(n-1)) t3<-n-t2 for(i in 2:N){ for(j in 1:(i-1)){ Lm2[1]<-1 Lm2[2]<-Lm[2]-(q[i]+q[j]) for(m in 3:(n-1)){ Lm2[m]<-Lm[m]-(q[i]+q[j])*Lm2[m-1]-q[i]*q[j]*Lm2[m-2] } piij[i,j]<-Kn*q[i]*q[j]*sum((t2+1-n*(p[i]+p[j]))*Lm2[t3]/n^(t2-1)) piij[j,i]<-piij[i,j] } piij[i,i]<-n*p[i] } piij[1,1]<-n*p[1]