##################################################################
#                                                                #
# 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]