第一次运行结果,4e4
细节图：
说明：
图中红色线条为 loglogN,以e为底,注意绘图时其横坐标为对数,纵坐标为均匀坐标.
蓝色线条为熵率.
绿色线条为概率值p1.
计算出的熵率最大值为 0.5 左右,与参考答案中的 1 值差 2 倍,估计是熵率的计算方法在理解上有偏差.
第二次运行结果,4e4+1~8e4
说明：
将此幅图拼接到第一幅后面,就可以得到答案中完整的图.
但是,运算量较大,上述运算量为4e4,若要达到1e25,则需要重复执行2.5e20,计算量相当可观.不过,可以从前期结果分析出熵率的变化趋势和概率值的变化趋势的关系.
MATLAB代码：
close all;
clear all;
clc;
% 由于所要求长度过长,无法一次实现,可以采用数据拼接的方法.
temp = 0;% 计算熵时的临时变量.
Len = 4e4;% 总长度.
i = 1: 1: Len;
logi = log(log(i));
randi = rand(1, Len);
for i = 1: 1: Len
x05line1(i) = 1;
x05if i <= 2
x05x05X(i) = 0;
x05x05P1(i) = 0;
x05x05continue;
x05end;
x05if mod(floor(logi(i)), 2) == 0 mod(ceil(logi(i)), 2) ~= 0
x05x05p1(i) = 0.5;
x05x05if randi(i) > 0.5
x05x05x05X(i) = 1;
x05x05else
x05x05x05X(i) = 0;
x05x05end;
x05else
x05x05p1(i) = 0;
x05x05X(i) = 0;
x05end;
end;
% p1 和 Xi 产生完毕,可以根据需要进行存储.
% 下面计算熵率,
for i = 1: 1: Len
x05temp = temp - p1(i) * log2(p1(i) + eps);% log = ln,不确定是以什么为底.
x05HN(i) = temp / i;% * 2 ?
end;
% 画图,注意参考答案的熵率横坐标为对数形式.
figure;
x = 1: 1: Len;
plot(x, p1, 'g');hold on;% p1 值
plot(x, logi, 'r');hold on;% log(log(N))
plot(x, line1, 'k-');hold on;% 渐近线
plot(x, HN, 'b');% 熵率
set(gca, 'XScale', 'log');
axis([1e0 1e25 0 4.1]);
hold off;grid off;