个settings

下分别测得的数据

，其中每组数据

是一个随机离散时间序列

，其中

。然而理论合作者那边要求获得的时间序列要满足每个数据点所对应的测量setting必须是从那

个settings中按照给定概率分布随机抽取的，那我必须生成一条随机setting序列

，其中

，然后用我那测得的

组数据按照这个setting序列生成一条新的随机序列：也就是说假设这个setting序列里面有

个是

，那么我就用数据

中的前

个点往对应的位置填充，以此类推，那么生成的随机序列

就满足了理论合作者的要求。

以后，还需要用它来获得两条新序列

和

用来实验结果的曲线图。第一条序列是这样的：对于第

个点，先求得前

个数据点的平均值

，然后求

，其中

是KL散度函数，

是一个固定不变的值。第二条序列是这样的：对于第

个点，同样先求得前

个数据点的平均值

，然后求

是KL散度

在固定

后的逆函数

（这种条件下的逆函数是存在的，因为单调，只不过需要用数值方法求解）。

次

实验，也就是说上面过程重复个

次。

，再计算

和

。当时这个程序写出来是为了仿真一个低维高保真度的情况，

不需要太大（5000左右吧，

取50），跑的时间还是可以接受的（不到一分钟吧），毕竟就画个图给老板看。结果实验数据维度高了保真度低了，

得取到50000左右，结果运行时间直接飙升到十几分钟，但我还要对着跑出来图debug，这来回几下也太费时间了。

个settings上的概率分布为

，那这个随机序列生成函数可以这么写：

function randseq = getRandSeq(len, probs)
cutpoints = cumsum(probs); % 给不同结果划分区间
randseq = rand(1, len);
idx = randseq < cutpoints(1);   % 找到落在第一个区间的所有点的下标
randseq(idx) = 2;  % 给第一个区间上的点赋予对应的值
for k = 2 : max(size(probs))   % 对后面K-1个区间重复上面过程idx = randseq < cutpoints(k) & randseq >= cutpoints(k-1);randseq(idx) = k+1;
end
randseq = randseq - 1;
end

results = zeros(1, N);       % 给生成的随机序列预分配空间
% 生成随机setting序列
setting_seq = getRandSeq(N, probs);
% 用原始数据填充setting序列
for k = 1 : setting_num    idx = (setting_seq == k);    % 在setting序列中找出所有S_k对应的下标results(idx) = raw_data{k}(1:sum(idx));   % 用对应原始数据往对应位置填充
end

序列同样可以向量化：

results = cumsum(results) ./ (1:N);  % 直接在原地更新，不新建数组了

也给它向量化了：

function D = klDiv(p1, p0)
% 计算两个伯努利分布之间的KL散度
% p1和p2要么是两个长度相等的array，要么其中一个是数，另一个是array% 检查输入的范围
if sum(p0 < 0) > 0 || sum(p0 > 1) > 0 || sum(p1 > 1) > 0 || sum(p1 < 0) > 0error("illegal input.")
end
if max(size(p0)) == 1 % p1为array, p0为数值D = zeros(size(p1));  % 预分配内存indices0 = (p1 <= p0);  % 特殊情况0点下标indices1 = (p1 == 1);  % 特殊情况1点下标indices2 = (p1 > p0) & (p1 < 1);   % 正常情况点下标D(indices0) = 0;D(indices1) = log(1 / p0);D(indices2) = p1(indices2) .* log(p1(indices2) / p0) +...(1 - p1(indices2)) .* log((1 - p1(indices2)) / (1 - p0));
elseif max(size(p1)) == 1 % p0为array, p1为数值D = zeros(size(p0));  % 预分配内存if p1 == 1D(1:end) = log(1 ./ p0);elseindices0 = (p0 >= p1); % 特殊情况0点下标indices1 = (p0 < p1); % 正常情况点下标D(indices0) = 0;D(indices1) = p1 * log(p1 ./ p0) +...(1 - p1) * log((1 - p1) ./ (1 - p0));end
elseif max(size(p1)) == max(size(p0)) % p0和p1为等长度的arrayD = zeros(size(p1));  % 预分配内存indices0 = (p1 <= p0);  % 特殊情况0点下标indices1 = (p1 == 1);  % 特殊情况1点下标indices2 = (p1 > p0) & (p1 < 1);   % 正常情况点下标D(indices0) = 0;D(indices1) = log(1 ./ p0(indices1));D(indices2) = p1(indices2) .* log(p1(indices2) ./ p0(indices2)) +...(1 - p1(indices2)) .* log((1 - p1(indices2)) ./ (1 - p0(indices3)));
elseerror("illegal input.")
end
end

序列就可以一口气完成了：

Delta = exp(-klDiv(results, p0) .* (1:N));

序列的计算，因为每个点都需要用数值方法求解（我这用的是二分法，因为函数是单调的），当时我的思路是能不能想出什么近似方法把KL散度的逆函数显式表达出来，这样就可以向量化计算了，结果没想出来。。。后来就直接把二分法这个过程给向量化了：

function x = invKLDiv(ps, y)
threshold = 1e-4;  % 计算精度
left = zeros(size(ps));
right = ps;
x = (left + right) / 2;
res = klDiv(ps, x) - y;
idx = (right - left) > threshold;   % 找出精度仍然不达标的所有下标
idx1 = (res > 0) & idx; % 找出idx中需要更新区间左边界的所有下标
idx2 = (res <= 0) & idx;   % % 找出idx中需要更新区间右边界的所有下标
while sum(idx) > 0   % 循环直到所有点都达到精度要求left(idx1) = x(idx1);right(idx2) = x(idx2);x(idx) = (left(idx) + right(idx)) / 2;res(idx) = klDiv(ps(idx), x(idx)) - y(idx);idx = (right - left) > threshold;idx1 = (res > 0) & idx;idx2 = (res <= 0) & idx;
end
end

序列也可以一口气计算了：

eps_A = 16 * (1 - invKLDiv(results, -log(delta)./ (1:N))) / 3;

次

实验数据后处理给向量化了，很简单，就是把上面生成的随机setting序列的长度从

拓展到

就是了，最后整个过程长这样：

% start post-processing the raw data
one_to_N = 1 : N;
one_to_N_M = repmat(one_to_N, 1, M);
results = zeros(1, N*M);       % test results sequence
% generate a random setting sequence
setting_seq = getRandSeq(N*M, probs);
% fill the results sequence by raw data according to setting_seq
for k = 1 : setting_num    idx = (setting_seq == k);results(idx) = raw_data{k}(1:sum(idx));
end
% compute the significance and infidelity sequence
for k = 1 : Mresults((k-1)*N+1 : k*N) = cumsum(results((k-1)*N+1 : k*N)) ./ one_to_N;
end
Delta = exp(-klDiv(results, p0) .* one_to_N_M);
eps_A = d * (1 - invKLDiv(results, -log(delta) ./ one_to_N_M)) / ((d + 1) * nu);
% store the M i.i.d. exp data
for k = 1 : Mexp_data{4*(k - 1) + 1} = results(k*N);exp_data{4*(k - 1) + 2} = Delta((k-1)*N+1 : k*N);exp_data{4*(k - 1) + 3} = eps_A((k-1)*N+1 : k*N);exp_data{4*(k - 1) + 4} = results((k-1)*N+1 : k*N);
end