#计算合法近邻数 nneighb <- function(x,a){ nn <- 1 sum0 <- 0 sumb <- 0 for(i in 1:length(x)){ #要判断的条件的当前值 sum0 <- sum0 + i*x[i] } #交换 x[i]和 x[j]的位置看是否符合条件 for(i in 1:(length(x)-1)){ for(j in (i+1):length(x)){ sumb = sum0 - i * x[i] - j * x[j] + i * x[j] + j * x[i] if(sumb >= a){nn = nn + 1} } } return(nn) } demo_permsub <- function(x,nout=10,n=100,a=1000){ sum1 <- 0 for(i in n){ #要判断的条件的当前值 sum1 <- sum1 + i*x[i] } ## 输出为 nout 行 n 列矩阵 resmat <- matrix(0,nout,n) nsucc <- 0 nfail <- 0 ntry <- 0 nn1 <- nneighb(x,a) while(nsucc < nout){ #已经成功找到的元素个数 #上一次成功后的累计失败次数 #总试投次数 #当前合法近邻数 ntry <- ntry +1 nfail <- nfail +1 if(nfail>1000){ cat('尝试失败超过 1000 次： nsucc',nsucc,'ntry=',ntry,'nfail=',nfail,'\n') print(x) break } #进行一次尝试，随机选取两个元素下标，交换对应元素的位置 #看是否合法元素。 

I1 <- ceiling(runif(1,0,1)*n) I2 <- ceiling(runif(1,0,1)*(n-1)) if(I2 == I1){ #第一个随机下标 #第二个随机下标 I2 <- I2 + 1 } #按照新的位置，计算目标的新值，比较是否超过界限 sum2 <- sum1 - I1*x[I1] -I2*x[I2] + I1*x[I2] + I2*x[I1] if(sum2 >= a){ #交换得到新的合法元素，看是否要转移 #先跳过去，如果不满足概率要求再跳回来 tmp1 <- x[I1] x[I1] <- x[I2] x[I2] <- tmp1 tran <- FALSE nn2 <- nneighb(x,a) if(nn2 <= nn1){ tran <- TRUE } else{ R <- runif(1,0,1) if(R < nn1 / nn2){ tran <- TRUE } } if(tran == TRUE){ nsucc <- nsucc +1 nfail <- 0 resmat[nsucc,] = x nn1 <- nn2 sum1 <- sum2 } else{ tmp1 = x[I1] x[I1] = x[I2] x[I2] = tmp1 nfail <- nfail +1 } } else nfail <-nfail + 1 } return(resmat) } 

#测试 #x：初始排列； n：排列的元素个数； a：界限； nout：需要输出的抽样次数 n=100 a=10000 nout=100 set.seed(200) x <- sample(1:n,n,replace = FALSE) demo_permsub(x,nout,n,a)