对数间隔的整数
Posted
技术标签:
【中文标题】对数间隔的整数【英文标题】：logarithmically spaced integers 【发布时间】：2012-09-07 05:56:41 【问题描述】：假设我有一个 10,000 pt 的向量，我想取一个只有 100 个对数间隔点的切片。我想要一个函数来给我索引的整数值。这是一个简单的解决方案，只需使用 around + logspace，然后摆脱重复项。
def genLogSpace( array_size, num ):
lspace = around(logspace(0,log10(array_size),num)).astype(uint64)
return array(sorted(set(lspace.tolist())))-1
ls=genLogspace(1e4,100)
print ls.size
>>84
print ls
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 13, 14, 15, 17, 19, 21, 23, 25, 27, 30,
33, 37, 40, 44, 49, 54, 59, 65, 71, 78, 86,
94, 104, 114, 125, 137, 151, 166, 182, 200, 220, 241,
265, 291, 319, 350, 384, 422, 463, 508, 558, 613, 672,
738, 810, 889, 976, 1071, 1176, 1291, 1416, 1555, 1706, 1873,
2056, 2256, 2476, 2718, 2983, 3274, 3593, 3943, 4328, 4750, 5213,
5721, 6279, 6892, 7564, 8301, 9111, 9999], dtype=uint64)
注意有 16 个重复，所以现在我只有 84 分。
有没有人可以有效地确保输出样本数为 num 的解决方案？对于这个具体的例子，121 和 122 的 num 的输入值给出 100 个输出点。
【问题讨论】：
您的问题是因为logspace 返回均匀分布的样本。
一般情况下，当(num+1) 是 2 的幂时，您只会拥有完全对数间隔的索引。观察上面的结果：前 15 个点实际上是完全线性间隔的。
@chthonicdaemon：你是怎么得到这条规则的（num+1 是2 的力量）？从技术上讲，如果array_size ** (1/(num-1)) 是整数，则可以使用精确的整数对数索引（假设索引从1 开始并以array_size 结束）。
好收获！我是从第一个 /interval/ 等于 1 推断的。当然，第一个间隔可以是任何整数，但其余的间隔必须从那里加倍。
【参考方案1】：
这有点棘手。你不能总是得到对数间隔的数字。如您的示例所示，第一部分是相当线性的。如果你同意，我有一个解决方案。但是对于解决方案，您应该了解为什么会有重复。
对数刻度满足条件：
s[n+1]/s[n] = constant
我们称这个常量为r 为ratio。对于n 范围内1...size 之间的这些数字，您将得到：
1, r, r**2, r**3, ..., r**(n-1)=size
所以这给了你：
r = size ** (1/(n-1))
在您的情况下，n=100 和 size=10000，r 将是 ~1.0974987654930561，这意味着，如果您以 1 开头，您的下一个数字将是 1.0974987654930561，然后四舍五入为 @987654335 @ 再次。因此你的重复。少数人会出现此问题。在一个足够大的数字之后，与 ratio 相乘会得到一个不同的四舍五入整数。
记住这一点，最好的办法是将连续整数相加到某个点，这样与比率的乘法就不再是问题了。然后你可以继续对数缩放。以下函数可以做到这一点：
import numpy as np
def gen_log_space(limit, n):
result = [1]
if n>1: # just a check to avoid ZeroDivisionError
ratio = (float(limit)/result[-1]) ** (1.0/(n-len(result)))
while len(result)<n:
next_value = result[-1]*ratio
if next_value - result[-1] >= 1:
# safe zone. next_value will be a different integer
result.append(next_value)
else:
# problem! same integer. we need to find next_value by artificially incrementing previous value
result.append(result[-1]+1)
# recalculate the ratio so that the remaining values will scale correctly
ratio = (float(limit)/result[-1]) ** (1.0/(n-len(result)))
# round, re-adjust to 0 indexing (i.e. minus 1) and return np.uint64 array
return np.array(list(map(lambda x: round(x)-1, result)), dtype=np.uint64)
Python 3 更新：最后一行曾经是 return np.array(map(lambda x: round(x)-1, result), dtype=np.uint64) 在 Python 2 中
以下是一些使用它的示例：
In [157]: x = gen_log_space(10000, 100)
In [158]: x.size
Out[158]: 100
In [159]: len(set(x))
Out[159]: 100
In [160]: y = gen_log_space(2000, 50)
In [161]: y.size
Out[161]: 50
In [162]: len(set(y))
Out[162]: 50
In [163]: y
Out[163]:
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11,
13, 14, 17, 19, 22, 25, 29, 33, 38, 43, 49,
56, 65, 74, 84, 96, 110, 125, 143, 164, 187, 213,
243, 277, 316, 361, 412, 470, 536, 612, 698, 796, 908,
1035, 1181, 1347, 1537, 1753, 1999], dtype=uint64)
为了向您展示结果的对数程度，这里是x = gen_log_scale(10000, 100) 的输出的半对数图（如您所见，左侧不是真正的对数）：
【讨论】：
【参考方案2】：Avaris's answer 中直接生成对数间隔点的方法绝对是可行的方法。但我认为看看如何选择适当的值传递给logspace 以获得你想要的东西会很有趣。
logspace(0, k, n) 生成的数组中的值是数字 10ik / (n−1) for 0 ≤ i n:
>>> numpy.logspace(0, 2, 10)
array([ 1. , 1.66810054, 2.7825594 , 4.64158883,
7.74263683, 12.91549665, 21.5443469 , 35.93813664,
59.94842503, 100. ])
>>> [10 ** (i * 2 / 9.0) for i in xrange(10)]
[1.0, 1.6681005372000588, 2.7825594022071245, 4.641588833612778,
7.742636826811269, 12.91549665014884, 21.544346900318832,
35.938136638046274, 59.94842503189409, 100.0]
此序列由一个初始段组成，其中值比单位间距更接近（因此当它们四舍五入到最接近的整数时可能会出现重复），然后是值比单位间距更宽的段和没有重复。
>>> ' '.join(':.2f'.format(10 ** (i * 2 / 19.0)) for i in xrange(20))
'1.00 1.27 1.62 2.07 2.64 3.36 4.28 5.46 6.95 8.86 11.29 14.38 18.33 23.36
29.76 37.93 48.33 61.58 78.48 100.00'
>>> [int(0.5 + 10 ** (i * 2 / 19.0)) for i in xrange(20)]
[1, 1, 2, 2, 3, 3, 4, 5, 7, 9, 11, 14, 18, 23, 30, 38, 48, 62, 78, 100]
值之间的间距为s(i) = 10iK − 10(i-1)K，其中K = k / (n - 1)。令 m 为满足 s(m) ≥ 1 的最小值。（示例中为 m = 7上面。）然后当删除重复项时，正好有 ⌊½ + 10(m−1)K⌋ + n - m 个剩余数字。
一点代数发现：
m = ⌈ − log(1 − 10−K) / K log 10 ⌉
让我们检查一下。
from math import ceil, floor, log
def logspace_size(k, n):
"""
Return the number of distinct integers we'll get if we round
`numpy.logspace(0, k, n)` to the nearest integers and remove
duplicates.
>>> logspace_size(4, 100)
84
>>> logspace_size(4, 121)
100
>>> from numpy import around, logspace
>>> all(logspace_size(k, n) == len(set(around(logspace(0, k, n))))
... for k in xrange(1,10) for n in xrange(2,100))
True
"""
K = float(k) / (n - 1)
m = int(ceil(- log(1 - 10 ** -K) / (K * log(10))))
if m < n:
return int(0.5 + 10 ** ((m - 1) * K)) + n - m
else:
return int(0.5 + 10 ** ((n - 1) * K))
文档测试通过了，所以这对我来说看起来不错。所以你需要做的就是找到n 这样logspace_size(4, n) == 100。您可以通过二进制印章或scipy.optimize 方法之一来做到这一点：
>>> f = lambda x, k, n:(logspace_size(k, x) - n)**2
>>> int(round(scipy.optimize.fmin(f, 100, args=(4,100), xtol=0.5, ftol=0.5)[0]))
Optimization terminated successfully.
Current function value: 0.015625
Iterations: 8
Function evaluations: 17
122
【讨论】：
【参考方案3】：我在这里搜索一个简单的方法来在 python 中获取对数间隔系列（以 10 为底）（省略使用 numpy）。但是对于我超简单的需求，您的解决方案太复杂了。
def logarithmic_decade(numbers_per_decade, offset=10):
for n in xrange(numbers_per_decade):
yield offset * 10.0 ** (n / float(numbers_per_decade))
由于它是生成器，因此您必须：
numbers = list(logarithmic_decade(5))
print numbers
[10.0, 15.848931924611136, 25.118864315095802, 39.81071705534972, 63.095734448019336]
for p, n in zip(numbers, numbers[1:] + [100]):
print 'prev = p:.2f, next = n:.2f, next/prev = rt:.4f'.format(p=p, n=n, rt=n / p)
给出以下输出：
prev = 10.00, next = 15.85, next/prev = 1.5849
prev = 15.85, next = 25.12, next/prev = 1.5849
prev = 25.12, next = 39.81, next/prev = 1.5849
prev = 39.81, next = 63.10, next/prev = 1.5849
prev = 63.10, next = 100.00, next/prev = 1.5849
【讨论】：
【参考方案4】：1到1e4之间的单行解：
y = [(lambda x:int(x))(tmp) for tmp in np.logspace(0,4, 10)]
【讨论】：
以上是关于对数间隔的整数的主要内容，如果未能解决你的问题，请参考以下文章
R语言生成对数线性间隔数据序列(Log-linearly Spaced Sequences)