有没有一种有效的方法可以在具有给定总和或平均值的范围内生成 N 个随机整数?

Is there an efficient way to generate N random integers in a range that have a given sum or average?

提问人:Peter O. 提问时间:4/24/2020 最后编辑:Peter O. 更新时间:7/11/2022 访问量:3676

问:

有没有一种有效的方法可以生成 N 个整数的随机组合,以便——

  • 每个整数都在区间 [, ],minmax
  • 整数的总和为 ,sum
  • 整数可以按任何顺序出现(例如,随机顺序),并且
  • 组合是从满足其他要求的所有组合中随机统一选择的?

对于随机组合,是否有类似的算法,其中整数必须按其值(而不是按任何顺序)排序?

(选择具有均值的适当组合是一种特例,如果 .这个问题等效于生成一个均匀的随机分区,该部分分为 N 个部分,每个部分都在区间 [, ] 中,并按其值以任何顺序或排序顺序出现,视情况而定。meansum = N * meansumminmax

我知道这个问题可以通过以下方式解决以随机顺序出现的组合(编辑 [4 月 27 日]:算法修改。

  1. 如果 或 ,则没有解决方案。N * max < sumN * min > sum

  2. 如果 ,只有一个解,其中所有数字都等于 。如果 ,只有一个解,其中所有数字都等于 。N * max == sumNmaxN * min == sumNmin

  3. 使用 Smith 和 Tromble 中给出的算法(“从单位单纯形采样”,2004 年)生成 N 个随机非负整数,其总和为 。sum - N * min

  4. 添加到以这种方式生成的每个数字中。min

  5. 如果任何数字大于 ,请转到步骤 3。max

但是,如果远小于 ,则此算法会很慢。例如,根据我的测试(上面涉及的特殊情况的实现),算法平均拒绝——maxsummean

  • 大约 1.6 个样本,如果 ,但N = 7, min = 3, max = 10, sum = 42
  • 大约 30.6 个样本,如果 .N = 20, min = 3, max = 10, sum = 120

有没有办法修改此算法以在满足上述要求的同时对大 N 有效?

编辑:

作为评论中建议的替代方案,生成有效随机组合(满足除最后一个要求外的所有要求)的有效方法是:

  1. 计算 、 和 给定的可能有效组合数。Xsumminmax
  2. 选择 ,一个统一的随机整数。Y[0, X)
  3. 将 (“unrank”) 转换为有效的组合。Y

但是,是否有计算有效组合(或排列)数量的公式,是否有办法将整数转换为有效组合?[编辑(4 月 28 日):排列而不是组合相同]。

编辑(4 月 27 日):

在阅读了Devroye的Non-Uniform Random Variate Generation(1986)之后,我可以确认这是生成随机分区的问题。此外,第 661 页的练习 2(尤其是 E 部分)与此问题相关。

编辑(4 月 28 日):

事实证明,我给出的算法是统一的,其中涉及的整数是以随机顺序给出的,而不是按其值排序的顺序。由于这两个问题都引起了人们的普遍兴趣,因此我修改了这个问题,以寻求这两个问题的规范答案。

以下 Ruby 代码可用于验证均匀性的潜在解决方案(其中是候选算法):algorithm(...)

combos={}
permus={}
mn=0
mx=6
sum=12
for x in mn..mx
  for y in mn..mx
    for z in mn..mx
      if x+y+z==sum
        permus[[x,y,z]]=0
      end
      if x+y+z==sum and x<=y and y<=z
        combos[[x,y,z]]=0
      end
    end
  end
end

3000.times {|x|
 f=algorithm(3,sum,mn,mx)
 combos[f.sort]+=1
 permus[f]+=1
}
p combos
p permus

编辑(4 月 29 日):重新添加了当前实现的 Ruby 代码。

下面的代码示例是用 Ruby 给出的,但我的问题与编程语言无关:

def posintwithsum(n, total)
    raise if n <= 0 or total <=0
    ls = [0]
    ret = []
    while ls.length < n
      c = 1+rand(total-1)
      found = false
      for j in 1...ls.length
        if ls[j] == c
          found = true
          break
        end
      end
      if found == false;ls.push(c);end
    end
    ls.sort!
    ls.push(total)
    for i in 1...ls.length
       ret.push(ls[i] - ls[i - 1])
    end
    return ret
end

def integersWithSum(n, total)
 raise if n <= 0 or total <=0
 ret = posintwithsum(n, total + n)
 for i in 0...ret.length
    ret[i] = ret[i] - 1
 end
 return ret
end

# Generate 100 valid samples
mn=3
mx=10
sum=42
n=7
100.times {
 while true
    pp=integersWithSum(n,sum-n*mn).map{|x| x+mn }
    if !pp.find{|x| x>mx }
      p pp; break # Output the sample and break
    end
 end
}

算法 随机 语言无关

评论

0赞 user58697 4/24/2020
你能澄清一下你的第三个要求吗?您是否需要所有可能的组合(包括均值错误的组合)或所有有效组合(即均值正确的组合)之间的一致性?
0赞 Peter O. 4/24/2020
所有有效组合,即满足其他要求的所有组合。
0赞 גלעד ברקן 4/26/2020
如果我们有一种方法可以对限制为 [min, max] 中 N 个整数的总和的分区进行计数和取消排序,那么随机选择其中一个分区并取消排序是否代表均匀分布,这会比您当前的方法更有效吗?总和和 N 可以有多大?
0赞 Peter O. 4/27/2020
我不知道你说的“取消总和的划分”是什么意思,我也不知道有证据证明这样做会导致这个问题意义上的均匀分布。对于这个问题,两者实际上都是无限的(在合理范围内)。我正在寻找一个规范的答案,因为 Stack Overflow 上提出的许多问题都出现了潜在的问题,包括这个问题和这个问题@גלעדברקןsumN
0赞 גלעד ברקן 4/27/2020
如果我们给每个可能的组合一个“等级”(或索引),以所有组合的有序排列,那么“取消排名”将意味着在给定其等级(当然还有 N、min 和 max)的情况下生成组合。为什么从所有可能的组合中选择一个不符合均匀分布?

答:

3赞 rossum 4/25/2020 #1

我没有测试过这个,所以它不是一个真正的答案,只是要尝试的东西,太长了,不适合评论。从一个满足前两个条件的数组开始,然后使用它,这样它仍然满足前两个条件,但随机性要高得多。

如果均值为整数,则初始数组可以是 [4, 4, 4, ...4] 或者也许是 [3, 4, 5, 3, 4, 5, ...5、8、0] 或类似的东西。对于平均值 4.5,请尝试 [4, 5, 4, 5, ...4, 5].

接下来,在数组中选择一对数字和 。可能第一个数字应该按顺序取,就像费舍尔-耶茨洗牌一样,第二个数字应该随机选择。按顺序选择第一个数字可确保每个数字至少被选中一次。num1num2

现在计算和 .这些是两个数字到 和 边界的距离。设置为两个距离中较小的一个。这是允许的最大更改,不会将一个或另一个数字置于允许的限制之外。如果为零,则跳过此对。max-num1num2-minmaxminlimitlimit

在 [1, ] 范围内选择一个随机整数:调用它。我从可选择范围内省略了 0,因为它没有效果。测试可能表明,通过包含它,您可以获得更好的随机性;我不确定。limitchange

现在设置 和 .这不会影响平均值,数组的所有元素仍在所需的边界内。num1 <- num1 + changenum2 <- num2 - change

您需要至少运行一次整个阵列。测试应该显示您是否需要多次运行它才能获得足够随机的东西。

预计到达时间:包含伪代码

// Set up the array.
resultAry <- new array size N
for (i <- 0 to N-1)
  // More complex initial setup schemes are possible here.
  resultAry[i] <- mean
rof

// Munge the array entries.
for (ix1 <- 0 to N-1)  // ix1 steps through the array in order.

  // Pick second entry different from first.
  repeat
    ix2 <- random(0, N-1)
  until (ix2 != ix1)

  // Calculate size of allowed change.
  hiLimit <- max - resultAry[ix1]
  loLimit <- resultAry[ix2] - min
  limit <- minimum(hiLimit, loLimit)
  if (limit == 0)
    // No change possible so skip.
    continue loop with next ix1
  fi

  // Change the two entries keeping same mean.
  change <- random(1, limit)  // Or (0, limit) possibly.
  resultAry[ix1] <- resultAry[ix1] + change
  resultAry[ix2] <- resultAry[ix2] - change

rof

// Check array has been sufficiently munged.
if (resultAry not random enough)
  munge the array again
fi

评论

0赞 Peter O. 4/29/2020
我已经测试过了,不幸的是,无论我进行多少次迭代,您的算法都不会形成所有解决方案的均匀分布。
0赞 rossum 4/29/2020
那好吧。无论如何,这都值得一试。:(
9赞 John McClane 4/30/2020 #2

这是我在 Java 中的解决方案。它功能齐全,包含两个生成器:用于未排序分区和已排序分区。您的生成器也在类中实现以进行比较。该类按顺序枚举所有可能的分区(未排序或已排序,具体取决于参数)。我已经为所有这些生成器添加了全面的测试(包括您的测试用例)。 在大多数情况下,实现是可以自我解释的。如果您有任何问题,我会在几天内回答。PermutationPartitionGeneratorCombinationPartitionGeneratorSmithTromblePartitionGeneratorSequentialEnumerator

import java.util.Random;
import java.util.function.Supplier;

public abstract class PartitionGenerator implements Supplier<int[]>{
    public static final Random rand = new Random();
    protected final int numberCount;
    protected final int min;
    protected final int range;
    protected final int sum; // shifted sum
    protected final boolean sorted;

    protected PartitionGenerator(int numberCount, int min, int max, int sum, boolean sorted) {
        if (numberCount <= 0)
            throw new IllegalArgumentException("Number count should be positive");
        this.numberCount = numberCount;
        this.min = min;
        range = max - min;
        if (range < 0)
            throw new IllegalArgumentException("min > max");
        sum -= numberCount * min;
        if (sum < 0)
            throw new IllegalArgumentException("Sum is too small");
        if (numberCount * range < sum)
            throw new IllegalArgumentException("Sum is too large");
        this.sum = sum;
        this.sorted = sorted;
    }

    // Whether this generator returns sorted arrays (i.e. combinations)
    public final boolean isSorted() {
        return sorted;
    }

    public interface GeneratorFactory {
        PartitionGenerator create(int numberCount, int min, int max, int sum);
    }
}

import java.math.BigInteger;

// Permutations with repetition (i.e. unsorted vectors) with given sum
public class PermutationPartitionGenerator extends PartitionGenerator {
    private final double[][] distributionTable;

    public PermutationPartitionGenerator(int numberCount, int min, int max, int sum) {
        super(numberCount, min, max, sum, false);
        distributionTable = calculateSolutionCountTable();
    }

    private double[][] calculateSolutionCountTable() {
        double[][] table = new double[numberCount + 1][sum + 1];
        BigInteger[] a = new BigInteger[sum + 1];
        BigInteger[] b = new BigInteger[sum + 1];
        for (int i = 1; i <= sum; i++)
            a[i] = BigInteger.ZERO;
        a[0] = BigInteger.ONE;
        table[0][0] = 1.0;
        for (int n = 1; n <= numberCount; n++) {
            double[] t = table[n];
            for (int s = 0; s <= sum; s++) {
                BigInteger z = BigInteger.ZERO;
                for (int i = Math.max(0, s - range); i <= s; i++)
                    z = z.add(a[i]);
                b[s] = z;
                t[s] = z.doubleValue();
            }
            // swap a and b
            BigInteger[] c = b;
            b = a;
            a = c;
        }
        return table;
    }

    @Override
    public int[] get() {
        int[] p = new int[numberCount];
        int s = sum; // current sum
        for (int i = numberCount - 1; i >= 0; i--) {
            double t = rand.nextDouble() * distributionTable[i + 1][s];
            double[] tableRow = distributionTable[i];
            int oldSum = s;
            // lowerBound is introduced only for safety, it shouldn't be crossed 
            int lowerBound = s - range;
            if (lowerBound < 0)
                lowerBound = 0;
            s++;
            do
                t -= tableRow[--s];
            // s can be equal to lowerBound here with t > 0 only due to imprecise subtraction
            while (t > 0 && s > lowerBound);
            p[i] = min + (oldSum - s);
        }
        assert s == 0;
        return p;
    }

    public static final GeneratorFactory factory = (numberCount, min, max,sum) ->
        new PermutationPartitionGenerator(numberCount, min, max, sum);
}

import java.math.BigInteger;

// Combinations with repetition (i.e. sorted vectors) with given sum 
public class CombinationPartitionGenerator extends PartitionGenerator {
    private final double[][][] distributionTable;

    public CombinationPartitionGenerator(int numberCount, int min, int max, int sum) {
        super(numberCount, min, max, sum, true);
        distributionTable = calculateSolutionCountTable();
    }

    private double[][][] calculateSolutionCountTable() {
        double[][][] table = new double[numberCount + 1][range + 1][sum + 1];
        BigInteger[][] a = new BigInteger[range + 1][sum + 1];
        BigInteger[][] b = new BigInteger[range + 1][sum + 1];
        double[][] t = table[0];
        for (int m = 0; m <= range; m++) {
            a[m][0] = BigInteger.ONE;
            t[m][0] = 1.0;
            for (int s = 1; s <= sum; s++) {
                a[m][s] = BigInteger.ZERO;
                t[m][s] = 0.0;
            }
        }
        for (int n = 1; n <= numberCount; n++) {
            t = table[n];
            for (int m = 0; m <= range; m++)
                for (int s = 0; s <= sum; s++) {
                    BigInteger z;
                    if (m == 0)
                        z = a[0][s];
                    else {
                        z = b[m - 1][s];
                        if (m <= s)
                            z = z.add(a[m][s - m]);
                    }
                    b[m][s] = z;
                    t[m][s] = z.doubleValue();
                }
            // swap a and b
            BigInteger[][] c = b;
            b = a;
            a = c;
        }
        return table;
    }

    @Override
    public int[] get() {
        int[] p = new int[numberCount];
        int m = range; // current max
        int s = sum; // current sum
        for (int i = numberCount - 1; i >= 0; i--) {
            double t = rand.nextDouble() * distributionTable[i + 1][m][s];
            double[][] tableCut = distributionTable[i];
            if (s < m)
                m = s;
            s -= m;
            while (true) {
                t -= tableCut[m][s];
                // m can be 0 here with t > 0 only due to imprecise subtraction
                if (t <= 0 || m == 0)
                    break;
                m--;
                s++;
            }
            p[i] = min + m;
        }
        assert s == 0;
        return p;
    }

    public static final GeneratorFactory factory = (numberCount, min, max, sum) ->
        new CombinationPartitionGenerator(numberCount, min, max, sum);
}

import java.util.*;

public class SmithTromblePartitionGenerator extends PartitionGenerator {
    public SmithTromblePartitionGenerator(int numberCount, int min, int max, int sum) {
        super(numberCount, min, max, sum, false);
    }

    @Override
    public int[] get() {
        List<Integer> ls = new ArrayList<>(numberCount + 1);
        int[] ret = new int[numberCount];
        int increasedSum = sum + numberCount;
        while (true) {
            ls.add(0);
            while (ls.size() < numberCount) {
                int c = 1 + rand.nextInt(increasedSum - 1);
                if (!ls.contains(c))
                    ls.add(c);
            }
            Collections.sort(ls);
            ls.add(increasedSum);
            boolean good = true;
            for (int i = 0; i < numberCount; i++) {
                int x = ls.get(i + 1) - ls.get(i) - 1;
                if (x > range) {
                    good = false;
                    break;
                }
                ret[i] = x;
            }
            if (good) {
                for (int i = 0; i < numberCount; i++)
                    ret[i] += min;
                return ret;
            }
            ls.clear();
        }
    }

    public static final GeneratorFactory factory = (numberCount, min, max, sum) ->
        new SmithTromblePartitionGenerator(numberCount, min, max, sum);
}

import java.util.Arrays;

// Enumerates all partitions with given parameters
public class SequentialEnumerator extends PartitionGenerator {
    private final int max;
    private final int[] p;
    private boolean finished;

    public SequentialEnumerator(int numberCount, int min, int max, int sum, boolean sorted) {
        super(numberCount, min, max, sum, sorted);
        this.max = max;
        p = new int[numberCount];
        startOver();
    }

    private void startOver() {
        finished = false;
        int unshiftedSum = sum + numberCount * min;
        fillMinimal(0, Math.max(min, unshiftedSum - (numberCount - 1) * max), unshiftedSum);
    }

    private void fillMinimal(int beginIndex, int minValue, int fillSum) {
        int fillRange = max - minValue;
        if (fillRange == 0)
            Arrays.fill(p, beginIndex, numberCount, max);
        else {
            int fillCount = numberCount - beginIndex;
            fillSum -= fillCount * minValue;
            int maxCount = fillSum / fillRange;
            int maxStartIndex = numberCount - maxCount;
            Arrays.fill(p, maxStartIndex, numberCount, max);
            fillSum -= maxCount * fillRange;
            Arrays.fill(p, beginIndex, maxStartIndex, minValue);
            if (fillSum != 0)
                p[maxStartIndex - 1] = minValue + fillSum;
        }
    }

    @Override
    public int[] get() { // returns null when there is no more partition, then starts over
        if (finished) {
            startOver();
            return null;
        }
        int[] pCopy = p.clone();
        if (numberCount > 1) {
            int i = numberCount;
            int s = p[--i];
            while (i > 0) {
                int x = p[--i];
                if (x == max) {
                    s += x;
                    continue;
                }
                x++;
                s--;
                int minRest = sorted ? x : min;
                if (s < minRest * (numberCount - i - 1)) {
                    s += x;
                    continue;
                }
                p[i++]++;
                fillMinimal(i, minRest, s);
                return pCopy;
            }
        }
        finished = true;
        return pCopy;
    }

    public static final GeneratorFactory permutationFactory = (numberCount, min, max, sum) ->
        new SequentialEnumerator(numberCount, min, max, sum, false);
    public static final GeneratorFactory combinationFactory = (numberCount, min, max, sum) ->
        new SequentialEnumerator(numberCount, min, max, sum, true);
}

import java.util.*;
import java.util.function.BiConsumer;
import PartitionGenerator.GeneratorFactory;

public class Test {
    private final int numberCount;
    private final int min;
    private final int max;
    private final int sum;
    private final int repeatCount;
    private final BiConsumer<PartitionGenerator, Test> procedure;

    public Test(int numberCount, int min, int max, int sum, int repeatCount,
            BiConsumer<PartitionGenerator, Test> procedure) {
        this.numberCount = numberCount;
        this.min = min;
        this.max = max;
        this.sum = sum;
        this.repeatCount = repeatCount;
        this.procedure = procedure;
    }

    @Override
    public String toString() {
        return String.format("=== %d numbers from [%d, %d] with sum %d, %d iterations ===",
                numberCount, min, max, sum, repeatCount);
    }

    private static class GeneratedVector {
        final int[] v;

        GeneratedVector(int[] vect) {
            v = vect;
        }

        @Override
        public int hashCode() {
            return Arrays.hashCode(v);
        }

        @Override
        public boolean equals(Object obj) {
            if (this == obj)
                return true;
            return Arrays.equals(v, ((GeneratedVector)obj).v);
        }

        @Override
        public String toString() {
            return Arrays.toString(v);
        }
    }

    private static final Comparator<Map.Entry<GeneratedVector, Integer>> lexicographical = (e1, e2) -> {
        int[] v1 = e1.getKey().v;
        int[] v2 = e2.getKey().v;
        int len = v1.length;
        int d = len - v2.length;
        if (d != 0)
            return d;
        for (int i = 0; i < len; i++) {
            d = v1[i] - v2[i];
            if (d != 0)
                return d;
        }
        return 0;
    };

    private static final Comparator<Map.Entry<GeneratedVector, Integer>> byCount =
            Comparator.<Map.Entry<GeneratedVector, Integer>>comparingInt(Map.Entry::getValue)
            .thenComparing(lexicographical);

    public static int SHOW_MISSING_LIMIT = 10;

    private static void checkMissingPartitions(Map<GeneratedVector, Integer> map, PartitionGenerator reference) {
        int missingCount = 0;
        while (true) {
            int[] v = reference.get();
            if (v == null)
                break;
            GeneratedVector gv = new GeneratedVector(v);
            if (!map.containsKey(gv)) {
                if (missingCount == 0)
                    System.out.println(" Missing:");
                if (++missingCount > SHOW_MISSING_LIMIT) {
                    System.out.println("  . . .");
                    break;
                }
                System.out.println(gv);
            }
        }
    }

    public static final BiConsumer<PartitionGenerator, Test> distributionTest(boolean sortByCount) {
        return (PartitionGenerator gen, Test test) -> {
            System.out.print("\n" + getName(gen) + "\n\n");
            Map<GeneratedVector, Integer> combos = new HashMap<>();
            // There's no point of checking permus for sorted generators
            // because they are the same as combos for them
            Map<GeneratedVector, Integer> permus = gen.isSorted() ? null : new HashMap<>();
            for (int i = 0; i < test.repeatCount; i++) {
                int[] v = gen.get();
                if (v == null && gen instanceof SequentialEnumerator)
                    break;
                if (permus != null) {
                    permus.merge(new GeneratedVector(v), 1, Integer::sum);
                    v = v.clone();
                    Arrays.sort(v);
                }
                combos.merge(new GeneratedVector(v), 1, Integer::sum);
            }
            Set<Map.Entry<GeneratedVector, Integer>> sortedEntries = new TreeSet<>(
                    sortByCount ? byCount : lexicographical);
            System.out.println("Combos" + (gen.isSorted() ? ":" : " (don't have to be uniform):"));
            sortedEntries.addAll(combos.entrySet());
            for (Map.Entry<GeneratedVector, Integer> e : sortedEntries)
                System.out.println(e);
            checkMissingPartitions(combos, test.getGenerator(SequentialEnumerator.combinationFactory));
            if (permus != null) {
                System.out.println("\nPermus:");
                sortedEntries.clear();
                sortedEntries.addAll(permus.entrySet());
                for (Map.Entry<GeneratedVector, Integer> e : sortedEntries)
                    System.out.println(e);
                checkMissingPartitions(permus, test.getGenerator(SequentialEnumerator.permutationFactory));
            }
        };
    }

    public static final BiConsumer<PartitionGenerator, Test> correctnessTest =
        (PartitionGenerator gen, Test test) -> {
        String genName = getName(gen);
        for (int i = 0; i < test.repeatCount; i++) {
            int[] v = gen.get();
            if (v == null && gen instanceof SequentialEnumerator)
                v = gen.get();
            if (v.length != test.numberCount)
                throw new RuntimeException(genName + ": array of wrong length");
            int s = 0;
            if (gen.isSorted()) {
                if (v[0] < test.min || v[v.length - 1] > test.max)
                    throw new RuntimeException(genName + ": generated number is out of range");
                int prev = test.min;
                for (int x : v) {
                    if (x < prev)
                        throw new RuntimeException(genName + ": unsorted array");
                    s += x;
                    prev = x;
                }
            } else
                for (int x : v) {
                    if (x < test.min || x > test.max)
                        throw new RuntimeException(genName + ": generated number is out of range");
                    s += x;
                }
            if (s != test.sum)
                throw new RuntimeException(genName + ": wrong sum");
        }
        System.out.format("%30s :   correctness test passed%n", genName);
    };

    public static final BiConsumer<PartitionGenerator, Test> performanceTest =
        (PartitionGenerator gen, Test test) -> {
        long time = System.nanoTime();
        for (int i = 0; i < test.repeatCount; i++)
            gen.get();
        time = System.nanoTime() - time;
        System.out.format("%30s : %8.3f s %10.0f ns/test%n", getName(gen), time * 1e-9, time * 1.0 / test.repeatCount);
    };

    public PartitionGenerator getGenerator(GeneratorFactory factory) {
        return factory.create(numberCount, min, max, sum);
    }

    public static String getName(PartitionGenerator gen) {
        String name = gen.getClass().getSimpleName();
        if (gen instanceof SequentialEnumerator)
            return (gen.isSorted() ? "Sorted " : "Unsorted ") + name;
        else
            return name;
    }

    public static GeneratorFactory[] factories = { SmithTromblePartitionGenerator.factory,
            PermutationPartitionGenerator.factory, CombinationPartitionGenerator.factory,
            SequentialEnumerator.permutationFactory, SequentialEnumerator.combinationFactory };

    public static void main(String[] args) {
        Test[] tests = {
                            new Test(3, 0, 3, 5, 3_000, distributionTest(false)),
                            new Test(3, 0, 6, 12, 3_000, distributionTest(true)),
                            new Test(50, -10, 20, 70, 2_000, correctnessTest),
                            new Test(7, 3, 10, 42, 1_000_000, performanceTest),
                            new Test(20, 3, 10, 120, 100_000, performanceTest)
                       };
        for (Test t : tests) {
            System.out.println(t);
            for (GeneratorFactory factory : factories) {
                PartitionGenerator candidate = t.getGenerator(factory);
                t.procedure.accept(candidate, t);
            }
            System.out.println();
        }
    }
}

您可以在 Ideone 上尝试此操作

评论

0赞 Peter O. 4/30/2020
感谢您的回答;它工作得很好。我在这里的另一个答案中描述了排列生成器;在你的帮助下回答了另一个问题;并将很快将您的算法包含在我关于随机生成方法的文章的 Python 示例代码中。
0赞 Joseph Wood 5/2/2020
只是为了清楚。该算法是否依赖于生成所有可能的分区/组合来进行采样?
0赞 John McClane 5/2/2020
@JosephWood 不,它依赖于计算所有这些。这在生成器初始化时只执行一次,并且非常有效,因为它利用了动态编程方法。
0赞 Peter O. 5/2/2020
动态规划如何解决将“sum”选择均匀随机划分为随机选择的 N 个整数的相关问题,从列表中替换(示例)或不替换示例),或者如何解决该问题
1赞 John McClane 12/19/2020
@Will 你是说?它是在构造函数中预先计算的表,然后在方法中用于生成随机分区。 计算从 0 到 range = max - min(含)的 n 个数字序列的数个序列,总和为 s。为了在我们已经找到具有更高指数的项之后生成第 i 项,我们将(这是某个区间内 for s 的总和)乘以随机数 unif。分布在 [0,1) 中,然后查找最高值 s(新的项和),使得乘积 t 小于 的累积和。distributionTableget()d.t.[n][s]d.t.[i + 1][s]d.t.[i][s]d.t.[i][s]
6赞 Peter O. 4/30/2020 #3

这是 John McClane 的 PermutationPartitionGenerator 中的算法,在本页的另一个答案中。它有两个阶段,即设置阶段和采样阶段,并在 [, ] 中生成随机变量,其总和为 ,其中数字以随机顺序列出。nminmaxsum

设置阶段:首先,使用以下公式生成解决方案表(其中位于 [0, ] 中,位于 [0, ] 中):t(y, x)ynxsum - n * min

  • t(0, j) = 1 如果 j == 0;否则为 0
  • t(i, j) = t(i-1, j) + t(i-1, j-1) + ... + t(i-1, j-(最大值-最小值))

这里,t(y, x) 存储数字之和(在适当范围内)等于 的相对概率。该概率相对于所有具有相同 的 t(y, x) 而言。yxy

抽样阶段:在这里,我们生成一个数字样本。设置为 ,然后对于每个位置,从 0 开始并向后工作到 0:nssum - n * minin - 1

  • 设置为 [0, t(i+1, s)) 中的统一随机整数。v
  • 设置为 。rmin
  • 从 中减去 t(i, s)。v
  • 当保持 0 或大于 0 时,从 中减去 t(i, s-1),将 1 加到 ,然后从 中减去 1。vvrs
  • 样本中位置处的数字设置为 。ir

编辑:

似乎通过对上述算法进行微不足道的更改,可以让每个随机变量使用单独的范围,而不是对所有变量使用相同的范围:

位置 ∈ [0, ) 处的每个随机变量都有一个最小值 min(i) 和一个最大值 max(i)。in

设 = - ∑min(i)。adjsumsum

设置阶段:首先,使用以下公式生成解决方案表(其中位于 [0, ] 中,位于 [0, ] 中):t(y, x)ynxadjsum

  • t(0, j) = 1 如果 j == 0;否则为 0
  • t(i, j) = t(i-1, j) + t(i-1, j-1) + ... + t(i-1, j-(最大(i-1)-最小(i-1)))

然后,采样阶段与之前完全相同,只是我们设置为 (而不是 ) 并设置为 min(i) (而不是 )。sadjsumsum - n * minrmin


编辑:

对于 John McClane 的 CombinationPartitionGenerator,设置和采样阶段如下。

设置阶段:首先,使用以下公式生成解决方案表(其中 is in [0, ]、 in [0, ]和 is in [0, ]):t(z, y, x)znymax - minxsum - n * min

  • t(0, j, k) = 1 如果 k == 0;否则为 0
  • t(i, 0, k) = t(i - 1, 0, k)
  • t(i, j, k) = t(i, j-1, k) + t(i - 1, j, k - j)

抽样阶段:在这里,我们生成一个数字样本。设置为 和 ,然后对于每个位置,从 0 开始并向后工作到 0:nssum - n * minmrangemax - minin - 1

  • 设置为 [0, t(i+1, mrange, s)) 中的统一随机整数。v
  • 设置为 min(,mrangemranges)
  • 从 中减去 。mranges
  • 设置为 。rmin + mrange
  • 从 中减去 t(, , )。imrangesv
  • 当保持 0 或大于时,将 1 加到 ,从中减去 1 和从 1 中减去 1,然后从 中减去 t(, , )。vsrmrangeimrangesv
  • 样本中位置处的数字设置为 。ir
3赞 Joseph Wood 5/4/2020 #4

正如 OP 所指出的,有效取消排名的能力非常强大。如果我们能够做到这一点,可以通过三个步骤生成均匀分布的分区(重申 OP 在问题中列出的内容):

  1. 计算长度为 N 的分区的总数 M,使各部分在 [, ] 范围内。summinmax
  2. 生成整数的均匀分布。[1, M]
  3. 将步骤 2 中的每个整数取消排列到其各自的分区中。

下面,我们只关注生成第 n 个分区,因为有大量关于在给定范围内生成均匀分布的整数的信息。这是一个简单的取消排名算法,它应该很容易翻译成其他语言(注意。我还没有弄清楚如何取消组合大小写的排名(即顺序问题))。C++

std::vector<int> unRank(int n, int m, int myMax, int nth) {

    std::vector<int> z(m, 0);
    int count = 0;
    int j = 0;

    for (int i = 0; i < z.size(); ++i) {
        int temp = pCount(n - 1, m - 1, myMax);

        for (int r = n - m, k = myMax - 1;
             (count + temp) < nth && r > 0 && k; r -= m, --k) {

            count += temp;
            n = r;
            myMax = k;
            ++j;
            temp = pCount(n - 1, m - 1, myMax);
        }

        --m;
        --n;
        z[i] = j;
    }

    return z;
}

主力函数由下式给出:pCount

int pCount(int n, int m, int myMax) {

    if (myMax * m < n) return 0;
    if (myMax * m == n) return 1;

    if (m < 2) return m;
    if (n < m) return 0;
    if (n <= m + 1) return 1;

    int niter = n / m;
    int count = 0;

    for (; niter--; n -= m, --myMax) {
        count += pCount(n - 1, m - 1, myMax);
    }

    return count;
}

此函数基于用户@m69_snarky_and_unwelcoming对有没有有效的算法进行有限数量的整数分区?上面给出的是对简单算法(没有记忆的算法)的轻微修改。这可以很容易地修改以合并记忆以提高效率。我们暂且不谈,重点介绍未排名的部分。

解释unRank

我们首先注意到,从长度为数字的 N 的分区(使得部分在 [, ] 范围内)到长度为 N 且部分在 [, ] 中的受限分区之间存在一对一映射。summinmaxsum - N * (min - 1)1max - (min - 1)

举个小例子,考虑 of length 的分区,使得 和 .这将具有与长度限制分区相同的结构,最大部分等于 。504min = 10max = 1550 - 4 * (10 - 1) = 14415 - (10 - 1) = 6

10   10   15   15   --->>    1    1    6    6
10   11   14   15   --->>    1    2    5    6
10   12   13   15   --->>    1    3    4    6
10   12   14   14   --->>    1    3    5    5
10   13   13   14   --->>    1    4    4    5
11   11   13   15   --->>    2    2    4    6
11   11   14   14   --->>    2    2    5    5
11   12   12   15   --->>    2    3    3    6
11   12   13   14   --->>    2    3    4    5
11   13   13   13   --->>    2    4    4    4
12   12   12   14   --->>    3    3    3    5
12   12   13   13   --->>    3    3    4    4

考虑到这一点,为了便于计数,如果您愿意,我们可以添加步骤 1a 将问题转换为“单位”情况。

现在,我们只有一个计数问题。正如@m69出色地显示的那样,通过将问题分解为更小的问题,可以很容易地实现对分区的计数。@m69提供的功能让我们完成了 90%,我们只需要弄清楚如何处理有上限的额外限制。这就是我们得到的地方:

int pCount(int n, int m, int myMax) {

    if (myMax * m < n) return 0;
    if (myMax * m == n) return 1;

我们还必须记住,随着我们的前进,这种情况会减少。如果我们看一下上面的第 6分区,这是有道理的:myMax

2   2   4   6

为了从现在开始计算分区的数量,我们必须继续将转换应用于“单位”情况。这看起来像:

1   1   3   5

在前面的步骤中,我们有一个最大值,现在我们只考虑一个最大值。65

考虑到这一点,取消分区排名与取消标准排列或组合的排名没有什么不同。我们必须能够计算给定部分中的分区数。例如,要计算以上述开头的分区数,我们所做的就是删除第一列中的:1010

10   10   15   15
10   11   14   15
10   12   13   15
10   12   14   14
10   13   13   14

10   15   15
11   14   15
12   13   15
12   14   14
13   13   14

转换为单位案例:

1   6   6
2   5   6
3   4   6
3   5   5
4   4   5

并致电:pCount

pCount(13, 3, 6) = 5

给定一个随机整数来取消排名,我们继续计算越来越小的部分中的分区数量(就像我们上面所做的那样),直到我们填满了索引向量。

例子

给定 、 、 和 ,下面是一个生成 20 个随机分区的 ideone 演示。输出如下:min = 3max = 10n = 7sum = 42

42: 3 3 6 7 7 8 8 
123: 4 4 6 6 6 7 9 
2: 3 3 3 4 9 10 10 
125: 4 4 6 6 7 7 8 
104: 4 4 4 6 6 8 10 
74: 3 4 6 7 7 7 8 
47: 3 4 4 5 6 10 10 
146: 5 5 5 5 6 7 9 
70: 3 4 6 6 6 7 10 
134: 4 5 5 6 6 7 9 
136: 4 5 5 6 7 7 8 
81: 3 5 5 5 8 8 8 
122: 4 4 6 6 6 6 10 
112: 4 4 5 5 6 8 10 
147: 5 5 5 5 6 8 8 
142: 4 6 6 6 6 7 7 
37: 3 3 6 6 6 9 9 
67: 3 4 5 6 8 8 8 
45: 3 4 4 4 8 9 10 
44: 3 4 4 4 7 10 10

词典索引在左边,未排名的分区在右边。

评论

1赞 Peter O. 5/4/2020
事实证明,这是一个非常好的选择,并且确实可以通过记忆变得有效。
1赞 גלעד ברקן 5/5/2020
对一对一映射的观察很好。
0赞 Lior Kogan 5/4/2020 #5

如果均匀生成 [l, x-1] 范围内的 0≤a≤1 个随机值,并均匀生成 [x, h] 范围内的随机值的 1-a,则预期均值为:

m = ((l+x-1)/2)*a + ((x+h)/2)*(1-a)

因此,如果你想要一个特定的 m,你可以用 a 和 x 来玩。

例如,如果设置 x = m:a = (h-m)/(h-l+1)。

为确保不同组合的概率更接近均匀,请从上述等式的一组有效解中随机选择 a 或 x。(x 必须在 [l, h] 范围内,并且应(接近)整数;N*a 也应该(接近)整数。

0赞 lostinafro 1/19/2021 #6

我为 Python-numpy 实现了(未排序的)算法,每个随机数都有单独的范围 [min, max]。也许它对使用 Python 作为主要编程语言的人很有用。

import numpy as np


def randint_sum_equal_to(sum_value: int, 
                         n: int, 
                         lower: (int, list) = 0, 
                         upper: (int,list) = None):

# Control on input
if isinstance(lower, (list, np.ndarray)):
    assert len(lower) == n
else:
    lower = lower * np.ones(n)
if isinstance(upper, (list, np.ndarray)):
    assert len(upper) == n
elif upper is None:
    upper = sum_value * np.ones(n)
else:
    upper = upper * np.ones(n)

# Trivial solutions
if np.sum(upper) < sum_value:
    raise ValueError('No solution can be found: sum(upper_bound) < sum_value')
elif np.sum(lower) > sum_value:
    raise ValueError('No solution can be found: sum(lower_bound) > sum_value')
elif np.sum(upper) == sum_value:
    return upper
elif np.sum(lower) == sum_value:
    return lower

# Setup phase
# I generate the table t(y,x) storing the relative probability that the sum of y numbers
# (in the appropriate range) is equal x.
t = np.zeros((n + 1, sum_value))
t[0, 0] = 1
for i in np.arange(1, n + 1):
    # Build the k indexes which are taken for each j following k from 0 to min(u(i-1)-l(i-1), j).
    # This can be obtained creating a repetition matrix of from t[i] multiplied by the triangular matrix
    # tri_mask and then sum each row
    tri_mask = np.tri(sum_value, k=0) - np.tri(sum_value, k=-(upper[i-1] - lower[i-1]))
    t[i] = np.sum(np.repeat(t[i-1][np.newaxis], sum_value, 0)*tri_mask, axis=1)

# Sampling phase
values = np.zeros(n)
s = (sum_value - np.sum(lower)).astype(int)
for i in np.arange(n)[::-1]:
    # The basic algorithm is the one commented:
    # v = np.round(np.random.rand() * t[i+1, s])
    # r = lower[i]
    # v -= t[i, s]
    # while (v >= 0) and (s > 0):
    #     s -= 1
    #     v -= t[i, s]
    #     r += 1
    # values[i] = r
    # ---------------------------------------------------- #
    # To speed up the convergence I use some numpy tricks.
    # The idea is the same of the Setup phase:
    # - I build a repeat matrix of t[i, s:1];
    # - I take only the lower triangular part, multiplying by a np.tri(s)
    # - I sum over rows, so each element of sum_t contains the cumulative sum of t[i, s - k]
    # - I subtract v - sum_t and count the element greater of equal zero,
    #   which are used to set the output and update s
    v = np.round(np.random.rand() * t[i+1, s])
    values[i] = lower[i]
    sum_t = np.sum(np.repeat(t[i, np.arange(1, s + 1)[::-1]][np.newaxis], s, 0) * np.tri(s), axis=1)
    vt_difference_nonzero = np.sum(np.repeat(v, s) - sum_t >= 0)
    values[i] += vt_difference_nonzero
    s -= vt_difference_nonzero
return values.astype(int)