Scrambled Halton Sequence(扰乱霍尔顿序列) 基于Python的代码实现

霍尔顿序列) 基于Python的代码实现"/>

Scrambled Halton Sequence(扰乱霍尔顿序列) 基于Python的代码实现

基于Python生成Scrambled Halton sequence

对生成思路的解释可以参考下面这篇文章,不过该作者给出的是普通Halton sequence,且基于java实现的.
Halton Sequence 原理和代码实现

普通Halton sequence在高维上会面临抽样相关性的问题,为了克服这种问题,我们就需要Scrambled Halton sequence,这种方法与普通Halton sequence的区别在于对断点进行了重新排序,这里我们采用Braaten, E.;Weller, G.(1979)所给出的排列组合.

具体文章参见:

[1]Braaten, E.;Weller, G…Improved low-discrepancy sequence for multidimensional quasi-Monte Carlo integration[J].J. Comput. Phys…

代码

这个版本的代码还有比较大的问题,因为采用进制转换的方法,导致无法在使用11作为其中一个base的情况下,计算超过10个元素.

一个可能的解决办法是用 S t + e i / p i S_{t}+e_{i}/p^i St​+ei​/pi这种方法来生成,难度应该也不大,为了偷懒我就不写了.

更多解释可以参考下面这本教材的第9章:

Kenneth E. Train, Discrete Choice Methods with Simulation, second edition, Cambridge University Press.
教材链接

# Author: LSY
# 由于有11作为质数,该程序只能在每个维度上计算不超过10个数
# 原因在于10在11进制下被转换为A,而A无法被替换为相应的数字进行计算
# 类似的还有11在12进制下被转换为B,等等.
import numpy as npdict_pe = {2: [0, 1],3: [0, 2, 1],5: [0, 3, 1, 4, 2],7: [0, 4, 2, 6, 1, 5, 3],11: [0, 5, 8, 2, 10, 3, 6, 1, 9, 7, 4],13: [0, 6, 10, 2, 8, 4, 12, 1, 9, 5, 11, 3, 7],17: [0, 8, 13, 3, 11, 5, 16, 1, 10, 7, 14, 4, 12, 2, 15, 6, 9],19: [0, 9, 14, 3, 17, 5, 11, 1, 15, 7, 12, 4, 18, 8, 2, 16, 10, 5, 13],23: [0, 11, 17, 4, 20, 7, 13, 2, 22, 9, 15, 5, 18, 1, 14, 10, 21, 6, 16, 3, 19, 8, 12],29: [0, 15, 7, 24, 11, 20, 2, 27, 9, 18, 4, 22, 13, 26, 5, 16, 10, 23, 1, 19, 28, 6, 14, 17, 3, 25, 12, 8,21]}  # Braatan所给出的Permutationdef scr_halton(counts):dic_result = {1: [], 2: [], 3: [], 4: [], 5: [], 6: [], 7: [], 8: [], 9: [], 10: []}  # 用于盛放结果dimension = 1for key in dict_pe:  # 依次计算十个维度list_result = []for count in range(1, counts+1):  # 每个维度计算counts次num = np.base_repr(count, base=key)  # 把第count个10进制数转为key进制print(count, "的", key, "进制转换结果为:", num)str_nums = str(num)inv_num = list(str_nums[::-1])  # 转换后的key进制取反print(count, "取反的结果:", list(inv_num))result = 0j = 1for i in inv_num:rep_i = dict_pe[key][int(i)]  # 取反后的数替换为Braatan所给出的相应数print(i, "应该被替换为:", rep_i)result += rep_i / (key ** j)  # 循环计算并得出相应结果j += 1result = round(result, 4)  # 最终结果保留4位小数list_result.append(result)dic_result[dimension] = list_resultdimension += 1return dic_result