求教高手,该怎么解决

求教高手
⑴已知二值随机分布序列{X1, X2, X3, …, Xi, …,}，在时刻i, Xi的取值为1的概率如下：
               /0.5, if 2k  < loglog i ≤ 2k + 1,
        p1(i)=<
               \0.0, if 2k+1< loglog i ≤ 2k + 2
    其中，k=0,1,2,3,4,...，且第一个X1=0;
    求该随机序列的熵率随i的变化情况－注为了体现变化特征
        要求：
        ①产生的序列长度接近3e8次方；
        ②概率p1的取值情况保存在"p1.txt"中，而序列保存在"Xi.txt"中；
        ③给出产生序列的MATLAB代码和对Xi序列进行熵率统计的代码；
        ④由于序列较长，不进行蒙特卡罗法平均，仅按e^(e^k)的长度或其他等长分组方式进行概率统计和熵率变化情况的演示即可。
------解决方案--------------------
第一次运行结果，1~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
line1(i) = 1;
if i <= 2
X(i) = 0;
P1(i) = 0;
continue;
end;
if mod(floor(logi(i)), 2) == 0 && mod(ceil(logi(i)), 2) ~= 0
p1(i) = 0.5;
if randi(i) > 0.5
X(i) = 1;
else
X(i) = 0;
end;
else
p1(i) = 0;
X(i) = 0;
end;
end;
% p1 和 Xi 产生完毕，可以根据需要进行存储。

% 下面计算熵率，
for i = 1: 1: Len
temp = temp - p1(i) * log2(p1(i) + eps);  % log = ln，不确定是以什么为底。
HN(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;


------解决方案--------------------
不知道怎么上传图片。。。