用MATLAB软件解首位数问题
杨亚辉
（海南软件职业技术学院，海南 琼海    571400）
【内容摘要】首位数问题是一个关于自然数分布特征的著名问题，若用纯数学的方法求解，则比较困难.本文以MATLAB软件为工具，设计了一个在MATLAB上可执行的程序，借助MATLAB强大的数值计算能力，巧妙的给出了首位数问题的近似解.为了减小误差我们在尽可能大的区间上产生尽可能多的随机自然数，并多次运行求解程序，对数据进行综合处理.最后得出结论：首位数为k(k=1,2,3,4,5,6,7,8,9)的自然数占全体自然数的比例依次为0.3596,0.1289,0.0865,0.0809, 0.0774,0.0734,0.0691,0.0644，0.0597.至此，我们对自然数的分布情况有了一个比较全面的认识.
【关 键 词】首位数问题； MATLAB；随机数；分布；统计;模拟实验
1. 引  言
天文学家西蒙·纽科姆 (1835-1909)在一次查对数表时，偶然发现了这样的现象：对数表开始的几页总要比后面几页磨损的厉害.这说明人们在查对数表时，较多地是使用了以1为首的那几页.于是，纽科姆便产生了这样一个疑问：首位数是1的自然数在全体自然数中占有多大的比例？它是不是要比首位数是其它数字的自然数要多？人们后来把这个问题称为“首位数问题”.
大家可能会认为这个问题是显而易见的.因为除0以外，共有9个数字：matlab生成随机数1,2,3,4,5,6,7,8,9用其中任何一个数字开头的自然数，在全体自然数中的分布是均匀的.这就是说，首位数为１的自然数应该占全体自然数的1/9.可是，事实并不这么简单. 1974年，美国斯坦福大学统计学家珀西·迪亚科尼斯已研究了这个问题，运用高等数学的方法，得出合理平均值，所得到的结论出乎人们的意料：首位数是1的自然数约占全体自然数的1/3. 准确一点说，这个数值应该是lg2约为0.30103.这是怎么一回事呢？
本文从统计学的角度出发，利用MATLAB强大的数值计算能力对首位数问题做出了一个比较全面的解答――验证珀西·迪亚科尼斯的结果，同时计算出其他首位数分别是2,3,4,5,6,7,8,9的自然数在全体自然数中的分布情况.
2.求解首位数问题的算法和程序
要确定基于首位数的自然数的分布规律，一个可行的办法就是用计算机模拟出尽可能多的随机自然数，然后统计出其中分别以k(k=1,2,3,4,5,6,7,8,9)为首位数的自然数的个数，从而确定它们各自在全体自然数中的分布比例.限于计算机系统和MATLAB的性能，在每次模拟中只产生200万个随机自然数.为了减小误差，重复进行了100次模拟实验.
求解首位数问题的算法如下：
第一步：初始化，a=0,k=1,2,3,…,9;
第二步：randn()产生常态变异数r;
第三步：将r转换成自然数；
第四步：max()求该组自然数的最大值；
第五步：length()统计该组自然数的个数n;
第六步：判断首位数是否以k开头，是a加1，否则a不变；
第七步：求出首位数为k的自然数占全体自然数的比例p=a/n;
第八步：输出p；
停  止
MATLAB工具箱提供了大量功能强大的函数.本文主要用到的MATLAB函数有(如表1):
表1
Randn()
产生常态变异数
Round()
整数取整
Jbtest()
正态分布的Jarque-Bera假设检验方法
mean()
求均值
Min()和max()
求最小值和最大值
Length()
统计矩阵中元素的个数
Bar()
二维条形图绘制
Find()
出符合特定条件数在矩阵中的位置
根据该思路编写出具体的求解程序，记为First_figure.m 文件.文件中%后的文字是对当前语句的说明或解释.其具体程序如下：
%% M-file, First_figure.m
for b=1:10          % 重复实验M =10*b次，本文重复实验了100次.
    for i=((b-1)*10+1):10*b  %分行显示,每行显示10个.
        a=0;          % a、b、i、l、m、n、R、x、z都是暂时变量.
        m=0;          % m表示自然数的位数，初始值取0.
        r=randn(U,1);      % U产生常态变异数的个数，U取值越大，产生的结果更精确，本文U=2000000.
        R=Q*r;        % Q是变异数，Q取值越大，产生随机数的随机性越强，本文取Q=1000000.
        l=find(R>0);  % find()函数出该组常态变异数大于0的位置.
        x=round(R(l));% round()函数将矩阵中的元素按最近的整数取整.
        z=max(x);    % z表示随机抽取自然数中的最大数.
        n=length(x);  % n表示随机抽取的自然数的个数.
        k=1;          % 首位数为k(依次取k＝1，2，3，4，5，6，7，8，9).
        while (10^m)<=z
            for i=1:n
                if x(i)>=k*(10^m)&x(i)<(k+1)*(10^m)
                    a=a+1;%统计以首位数为k的自然数的个数.
                end
            end
            m=m+1;
        end
        p=a/n;      % p是首位数为k的自然数占全体自然数比例.
        fprintf('%10.6f',p);
    end
    fprintf(' …\n');% …的作用是为了好以矩阵表示.
end
clear all
用MATLAB7.0中的mean()函数对First_figure.m运行的结果以矩阵表示进行均值统计(得到下表1),且分析发现First_figure.m运行的结果都在均值附近波动.为此，对First_figure.m运行结果用jbtest()进行分布正态性的Jarque-Bera假设检验方法检验，发现首位数为k的自然数在全体自然数中的分布服从正态分布.
用min()和max()函数对First_figure.m文件运行的结果统计出最小和最大值，发现首位数为k的自然数占全体自然数的比例的大致范围，具体如下表2所示：
表2
        所占比例 
首位数为k
均值
最小值
最大值
1
0.3596
0.3587
0.3609
2
0.1289
0.1279
0.1297
3
0.0865
0.0857
0.0872
4
0.0809
0.0804
0.0817
5
0.0774
0.0768
0.0780
6
0.0734
0.0729
0.0740
7
0.0691
0.0683
0.0697
8
0.0644
0.0637
0.0650
9
0.0597
0.0590
0.0602
为直观起见，仅取均值表示首位数为k的自然数的分布规律.如表3所示：
                            表3
k
1
2
3
4
5
6
7
8
9
均值
0.3596
0.1289
0.0865
0.0809
0.0774
0.0734
0.0691
0.0644
0.0597
根据表2中得出的数据，用bar()绘制二维条形图，相关绘图命令如下：
>> k=[1 2 3 4 5 6 7 8 9];
>> p=[0.3596 0.1289 0.0865 0.0809 0.0774 0.0734 0.0691 0.0644 0.0597];
>> bar(k,p);
>> xlabel('首位数为k的自然数');
>> ylabel('首位数为k的自然数占全体自然数的比例');
运行后得到图1.它直观的反映了首位数为k(k=1,2,3,4,5,6,7,8,9)的自然数占全体自然数的比例.
图 1
3.结论
以k＝1为例，我们用MATLAB模拟得出的结论（0.3596）和统计学家珀西·迪亚科尼斯用复杂的数学方法求得结论（0.30103）是一致的.同时我们还额外的求出了以2,3,4,5,6,7,8,9为首位数的自然数的分布比例.
即首位数为k(k=1,2,3,4,5,6,7,8,9)的自然数占全体自然数的比例分别为：0.3596, 0.1289, 0.0865, 0.0809, 0.0774, 0.0734,0.0691,0.0644,0.0597.其变化趋势按照k从1到9的顺序依次减小的.k＝1时，比例高达0.3596，k＝9时，仅为0.0597.而并非是均匀的分布（1/9）.
需要说明的是本文所得到的结果还是比较粗略的.这与本文使用的工具和方法有关： (1)计算机系统和MATLAB软件的性能限制了所取随机数的个数和运算的精度.(2)本文用MATLAB生成的随机数并不是真正的随机数,而是一种“伪随机数”.(3)本文用到的数学方法还是过于单一.要得到更精确的结果还得从以上三个方面改进.
首位数问题的结论在科学技术上发挥了重大的作用.例如，美国西雅图波音航天局数学家梅尔达德·沙沙哈尼在研究用计算机描绘自然景象的问题时，用上了这个结论.美国波音航天局还将这一成果用于飞机模拟器，使飞行员在不离开地面的情况下接受训练，而能得到一种在空中飞行的实感..
【参考文献】
[1] 薛定宇、陈阳泉，高等应用数学问题的MATLAB求解，清华大学出版社，2004年8月
[2] 刘树利、孙云龙、王家玉，计算机数学基础，高等教育出版社，2004年8月 [3] 叶其孝、姜启源、(美)吉奥丹诺，数学建模，机械工业出版社，2005年1月
[4] 高等数学（二）—线性代数与概率统计，复旦大学出版社，2000年1月
[5] 苏金明、刘宏、刘波.，MATLAB高级编程，电子工业出版，2005年8月
Solution First Figure Question
with MATLAB
YANG Ya-hui
  （Hainan Software Profession Institute，Hainan Qinghai 571400）
Abstract： The first figure question is a famous question about the natural number distrib
ution characteristic .It will be difficult to work out by pure mathematical method. But In my paper，it is be solved handily by the MATLAB with powerful ability about value computation. In order to reduce the error, we make as far as possible many random natural numbers in a as far as big interval, and run many times the solution procedure, carries on synthesis processing to the data. Finally draws the conclusion: the natural numbers which its first figure is k (k=1,2,3,4,5,6,7,8,9)  occupy all natural number the proportion is 0.3596,0.1289,0.0865,0.0809,0.0774,0.0734,0.0691,0.0644,0.0597.Heres, we have a quite comprehensive understanding about natural number's distribution.

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。