提问人:Peter O. 提问时间:4/24/2020 最后编辑:Peter O. 更新时间:7/11/2022 访问量:3676
有没有一种有效的方法可以在具有给定总和或平均值的范围内生成 N 个随机整数?
Is there an efficient way to generate N random integers in a range that have a given sum or average?
问:
有没有一种有效的方法可以生成 N 个整数的随机组合,以便——
- 每个整数都在区间 [, ],
min
max
- 整数的总和为 ,
sum
- 整数可以按任何顺序出现(例如,随机顺序),并且
- 组合是从满足其他要求的所有组合中随机统一选择的?
对于随机组合,是否有类似的算法,其中整数必须按其值(而不是按任何顺序)排序?
(选择具有均值的适当组合是一种特例,如果 .这个问题等效于生成一个均匀的随机分区,该部分分为 N 个部分,每个部分都在区间 [, ] 中,并按其值以任何顺序或排序顺序出现,视情况而定。mean
sum = N * mean
sum
min
max
我知道这个问题可以通过以下方式解决以随机顺序出现的组合(编辑 [4 月 27 日]:算法修改。
如果 或 ,则没有解决方案。
N * max < sum
N * min > sum
如果 ,只有一个解,其中所有数字都等于 。如果 ,只有一个解,其中所有数字都等于 。
N * max == sum
N
max
N * min == sum
N
min
使用 Smith 和 Tromble 中给出的算法(“从单位单纯形采样”,2004 年)生成 N 个随机非负整数,其总和为 。
sum - N * min
添加到以这种方式生成的每个数字中。
min
如果任何数字大于 ,请转到步骤 3。
max
但是,如果远小于 ,则此算法会很慢。例如,根据我的测试(上面涉及的特殊情况的实现),算法平均拒绝——max
sum
mean
- 大约 1.6 个样本,如果 ,但
N = 7, min = 3, max = 10, sum = 42
- 大约 30.6 个样本,如果 .
N = 20, min = 3, max = 10, sum = 120
有没有办法修改此算法以在满足上述要求的同时对大 N 有效?
编辑:
作为评论中建议的替代方案,生成有效随机组合(满足除最后一个要求外的所有要求)的有效方法是:
- 计算 、 和 给定的可能有效组合数。
X
sum
min
max
- 选择 ,一个统一的随机整数。
Y
[0, X)
- 将 (“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
}
答:
我没有测试过这个,所以它不是一个真正的答案,只是要尝试的东西,太长了,不适合评论。从一个满足前两个条件的数组开始,然后使用它,这样它仍然满足前两个条件,但随机性要高得多。
如果均值为整数,则初始数组可以是 [4, 4, 4, ...4] 或者也许是 [3, 4, 5, 3, 4, 5, ...5、8、0] 或类似的东西。对于平均值 4.5,请尝试 [4, 5, 4, 5, ...4, 5].
接下来,在数组中选择一对数字和 。可能第一个数字应该按顺序取,就像费舍尔-耶茨洗牌一样,第二个数字应该随机选择。按顺序选择第一个数字可确保每个数字至少被选中一次。num1
num2
现在计算和 .这些是两个数字到 和 边界的距离。设置为两个距离中较小的一个。这是允许的最大更改,不会将一个或另一个数字置于允许的限制之外。如果为零,则跳过此对。max-num1
num2-min
max
min
limit
limit
在 [1, ] 范围内选择一个随机整数:调用它。我从可选择范围内省略了 0,因为它没有效果。测试可能表明,通过包含它,您可以获得更好的随机性;我不确定。limit
change
现在设置 和 .这不会影响平均值,数组的所有元素仍在所需的边界内。num1 <- num1 + change
num2 <- 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
评论
这是我在 Java 中的解决方案。它功能齐全,包含两个生成器:用于未排序分区和已排序分区。您的生成器也在类中实现以进行比较。该类按顺序枚举所有可能的分区(未排序或已排序,具体取决于参数)。我已经为所有这些生成器添加了全面的测试(包括您的测试用例)。
在大多数情况下,实现是可以自我解释的。如果您有任何问题,我会在几天内回答。PermutationPartitionGenerator
CombinationPartitionGenerator
SmithTromblePartitionGenerator
SequentialEnumerator
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 上尝试此操作。
评论
distributionTable
get()
d.t.[n][s]
d.t.[i + 1][s]
d.t.[i][s]
d.t.[i][s]
这是 John McClane 的 PermutationPartitionGenerator 中的算法,在本页的另一个答案中。它有两个阶段,即设置阶段和采样阶段,并在 [, ] 中生成随机变量,其总和为 ,其中数字以随机顺序列出。n
min
max
sum
设置阶段:首先,使用以下公式生成解决方案表(其中位于 [0, ] 中,位于 [0, ] 中):t(y, x)
y
n
x
sum - 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) 而言。y
x
y
抽样阶段:在这里,我们生成一个数字样本。设置为 ,然后对于每个位置,从 0 开始并向后工作到 0:n
s
sum - n * min
i
n - 1
- 设置为 [0, t(i+1, s)) 中的统一随机整数。
v
- 设置为 。
r
min
- 从 中减去 t(i, s)。
v
- 当保持 0 或大于 0 时,从 中减去 t(i, s-1),将 1 加到 ,然后从 中减去 1。
v
v
r
s
- 样本中位置处的数字设置为 。
i
r
编辑:
似乎通过对上述算法进行微不足道的更改,可以让每个随机变量使用单独的范围,而不是对所有变量使用相同的范围:
位置 ∈ [0, ) 处的每个随机变量都有一个最小值 min(i) 和一个最大值 max(i)。i
n
设 = - ∑min(i)。adjsum
sum
设置阶段:首先,使用以下公式生成解决方案表(其中位于 [0, ] 中,位于 [0, ] 中):t(y, x)
y
n
x
adjsum
- 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) (而不是 )。s
adjsum
sum - n * min
r
min
编辑:
对于 John McClane 的 CombinationPartitionGenerator,设置和采样阶段如下。
设置阶段:首先,使用以下公式生成解决方案表(其中 is in [0, ]、 in [0, ]和 is in [0, ]):t(z, y, x)
z
n
y
max - min
x
sum - 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:n
s
sum - n * min
mrange
max - min
i
n - 1
- 设置为 [0, t(i+1, mrange, s)) 中的统一随机整数。
v
- 设置为 min(,
mrange
mrange
s
) - 从 中减去 。
mrange
s
- 设置为 。
r
min + mrange
- 从 中减去 t(, , )。
i
mrange
s
v
- 当保持 0 或大于时,将 1 加到 ,从中减去 1 和从 1 中减去 1,然后从 中减去 t(, , )。
v
s
r
mrange
i
mrange
s
v
- 样本中位置处的数字设置为 。
i
r
正如 OP 所指出的,有效取消排名的能力非常强大。如果我们能够做到这一点,可以通过三个步骤生成均匀分布的分区(重申 OP 在问题中列出的内容):
- 计算长度为 N 的分区的总数 M,使各部分在 [, ] 范围内。
sum
min
max
- 生成整数的均匀分布。
[1, M]
- 将步骤 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 且部分在 [, ] 中的受限分区之间存在一对一映射。sum
min
max
sum - N * (min - 1)
1
max - (min - 1)
举个小例子,考虑 of length 的分区,使得 和 .这将具有与长度限制分区相同的结构,最大部分等于 。50
4
min = 10
max = 15
50 - 4 * (10 - 1) = 14
4
15 - (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
在前面的步骤中,我们有一个最大值,现在我们只考虑一个最大值。6
5
考虑到这一点,取消分区排名与取消标准排列或组合的排名没有什么不同。我们必须能够计算给定部分中的分区数。例如,要计算以上述开头的分区数,我们所做的就是删除第一列中的:10
10
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 = 3
max = 10
n = 7
sum = 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
词典索引在左边,未排名的分区在右边。
评论
如果均匀生成 [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 也应该(接近)整数。
我为 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)
评论
sum
N