2024年2月22日发(作者：)
fortran产生随机数方法介绍（附代码）
注意：现在计算机产生的随机数都是伪随机数。
1.0-1之间均匀分布的随机数
random_number(x) 产生一个0到1之间的随机数（x可以是向量），但是每次总是那几个数。
用了random_seed ()后，系统根据日期和时间随机地提供种子，使得随机数更随机了。
program random
implicit none
real :: x
call random_seed () ! 系统根据日期和时间随机地提供种子
call random_number (x) ! 每次的随机数就都不一样了
write(*,*) x
stop
end program random
2.任意区间均匀分布的随机数
function my_random (lbound,ubound)
implicit none
real :: lbound,ubound
real :: len
real :: my_random
real :: t
len=ubound-lbound !计算范围大小
call random_number(t) !t是0-1之间的随机数
my_random=lbound+len*t
return
end
注意：在循环外call random_seed()
3.产生一个随机数数组，只需加一个循环即可
function my_random (lbound,ubound)
implicit none
real :: lbound,ubound
real :: len
integer size
real :: my_random(size) !size代表数组元素的个数
real :: t
integer i
len=ubound-lbound !计算范围大小
do i=1,10
call random_number(t) !t是0-1之间的随机数
my_random(i)=lbound+len*t !把t转换成lbound-ubound间的随机数
end do
return
end
注意：同理在循环外call random_seed()
4.标准正态分布随机数/高斯分布随机数
(1)徐士良的那本程序集里介绍了正态分布随机数产生的原理，不过他的方法只能产生较为简
单的随机数，随机数的质量并不高，特别是随机数的数目较多时。
(2)Box 和 Muller 在 1958 年给出了由均匀分布的随机变量生成正态分布的随机变量的算法。
设 U1, U2 是区间 (0, 1) 上均匀分布的随机变量，且相互独立。令
X1 = sqrt(-2*log(U1)) * cos(2*PI*U2);
X2 = sqrt(-2*log(U1)) * sin(2*PI*U2);
那么 X1, X2 服从 N(0,1) 分布，且相互独立。等于说我们用两个独立的 U(0,1) 随机数得到了
两个独立的 N(0,1)随机数。
值得说明的是，该方法产生的随机数质量很高。嘻嘻，本人亲自验证过。
SAS和蒙特卡罗模拟_随机数
一、为什么选择SAS做蒙特卡罗模拟？
为什么要用SAS做蒙卡？首先，对我来说，我只会用SAS，而且打算用SAS完成我所有的工作。当然，其他一些通用的理由有(Fan, etc.,2002)：
1. 蒙卡是个计算密集的活，而SAS Base、SAS Macro、SAS/IML强大而灵活的编程能力能满足这一要求；
2. 做蒙卡时要用到大量的统计/数学技术，而SAS就内置了大量的统计/数学函数(在 SAS/Stat和SAS/ETS)；用Fortran或C++当然也是非常好的主意，只是他们缺少内置的统计函数，代码就要冗长复杂很多。
二、什么是蒙卡？一个启发性例子
好，开始，什么是蒙卡？了解它背景知识的最好办法当然是wiki-Monte_Carlo_method。蒙特卡罗是位于摩洛哥的一家赌场，二战时，美国Los Alamos国家实验室把它作为核裂变计算机模拟的代码名称。作为模拟方法，蒙卡以前就叫统计抽样(statistical sampling)，我们感兴趣的结果因为输入变量的不确定而不可知，但如果能依概率产生输入变量的样本，我们就可以估计到结果变量的分布。跟蒙卡对应的，还有一种模拟技术叫系统模拟，包括排队、库存等模型，这些模型都跟随时间推移而出现的事件序列有关。下面举个蒙卡的例子，来自Evans, etc.( 2001)的超级简化版。
假设一家企业，利润是其需求量的函数，需求是随机变量。为了简化讨论，假定利润就是需求的两倍。这里输入变量就是不可控的需求，结果变量就是我们感兴趣的利润。假设需求以相同的概率取10、20、30、40、50、60这六种情况。在这样的简化下，我们就可以投一枚均匀的骰子来产生需求的样本，如果点数为1，对应得需求就是10，点数为2，需求就是20，以下类推。这样，我们的模拟过程就是：
1. 投骰子；
2. 根据骰子的点数确定需求量；
3. 根据需求量，求利润。
掷10次骰子，假设我们的模拟结果如下：
重复次数
1
2
3
4
5
6
7
8
9
10
骰子点数
5
3
3
6
1
3
4
5
2
5
需求量
50
30
30
60
10
30
40
50
20
50
利润
100
60
60
120
20
60
40
100
20
50
这样通过模拟需求的样本，我们对结果利润的分布也就有所知晓，比如平均利润可以算出就是63。蒙卡一个重要的步骤就是生成随机数，这里我们是用投骰子来完成。
三、又一个例子：利用蒙特卡罗模拟方法求圆周率∏(pi)
再举个很有名的例子，就是估计圆周率∏的值，来自Ross(2006)。这个试验的思路正好可以帮我们温习一下几何概型的概念。我们知道概率的古典概型，就是把求概率的问题转化为计数：样本空间由n个样本点组成，事件A由k个样本点组成，则事件A的概率就是k/n。考虑到概率和面积在测度上具有某种共性，几何概型的基本想法是把事件跟几何区域相对应，用面积来计算概率，其要点是：
1. 认为样本空间Ω是平面上的某个区域，其面积记为υ(Ω)；
2. 向区域Ω随机投掷一点，该点落入Ω内任何部分区域内的可能性只与这部分区域的面积成比例，而与这部分区域的位置和形状无关；
3. 设事件A是Ω的某个区域，面积为υ(A)，则向区域Ω上随机投掷一点，该点落在区域A的可能性称为事件A的概率，P(A)=υ(A)/υ(Ω)。
扯远了，回到用蒙卡估计圆周率∏的实验思路。假设样本空间是一个边长为2面积为4的正方形，我们感兴趣的事件是正方形内的一个半径为1面积为∏的圆，所以向正方形内随机投掷一点，落在圆里面的概率为∏/4。实验的思路如下：
1. 1）生成随机数——生成n个均匀落在正方形内的点；
2. 2）对落在正方形内的n个点，数一数正好落在圆里面的点的个数，假设为k（另外n-k个点就落在圆外面的正方形区域内）。
3. 3）k/n就可以大致认为是圆的面积与正方形的面积之比，另其等于∏/4，就可以求出圆周率∏的估计值。
又，Forcode提供了利用电子表格Excel求解以上问题的详细过程，有兴趣可以看看。另外，这里也有一个详细的说明，用Mathematica实现。
四、随机数很重要（以及下期预告）
在第一个例子中，我强调了骰子要是均匀的。那里骰子就是我们的随机数生成器，骰子的均匀程度就是我们随机数的“随机”成分。在上面的10次试验中，投出了3次点数“3”和3次点数“5”，然后点数“1”、“2”、“4”、“6”各出现一次——因为只是重复了10次，这样我们对骰子的均匀程度不好评估，但如果重复无数次——假如是十万次，如果还出现类似比例的结果，那么就有理由怀疑这粒骰子的均匀程度了。如果骰子灌了铅，比如投出点数“6”的可能性从六分之一上升到6分之三，那么根据这粒骰子来做模拟，我们就有高估需求以及利润的危险。
下次就讲讲随机数的生成原理，当然结合SAS来讲，或者也会提一下Excel。
五、其他软件包
在商业领域，基于Excel的插件比较兴盛，比如Crystal Ball应用就较为普遍，有试用版可以下载。其他竞争产品有@Risk、Risk Solver等。现在据说Risk
Simulator也开始兴盛起来，大有后来居上之势。他的开发者Johnathan Mun还在Wiley Finance系列有一本书，Modeling Risk: Applying Monte Carlo
Simulation, Real Options Analysis, Forecasting, and Optimization
Techniques。
S语言（R或者S-Plus）由于其面对对象的特性，加之丰富的内置函数和诸多用户提供的库，使得R或者S-Plus也是蒙卡研究的得力工具——或者比SAS更有前景。这个我不熟，略之。
Random Numbers
最简单、最基本、最重要的随机变量是在[0,1]上均匀分布的随机变量。一般地，我们把[0,1]上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都是借助于随机数来实现的。以下谈的都是所谓“伪随机数”(Pseudo Random
Numbers)。产生随机数，可以通过物理方法取得（很久很久以前，兰德公司就曾以随机脉冲源做信息源，利用电子旋转轮来产生随机数表），但当今最为普遍的乃是在计算机上利用数学方法产生随机数。这种随机数根据特定的迭代公式计算出来，初值确定后，序列就可以预测出来，所以不能算是真正的随机数（就成为“伪随机数”）。不过，在应用中，只要产生的伪随机数列能通过一系列统计检验，就可以把它们当成“真”随机数来用。
现在大多软件包内置的随机数产生程序，都是使用同余法(Congruential Random
Numbers Generators)。“同余”是数论中的概念。
0.预备知识：同余
捡回小学一年级的东西："4/2=2"读作“4除以2等于2”，或者，“2除4等于2”。还有求模的符号mod(number,divisor)，其中，number是被除数（在上式中，为4），divisor是除数（上式中的2）。这样的约定对SAS和Excel都通用，如mod(4,13)=4，mod(13,4)=1。
现在我们可以开讲“同余”了。设m是正整数，用m去除整数a、b，得到的余数相同，则称“a与b关于模m同余”。上面的定义可以读写成，对整数a、b和正整数m，若mod(a,m)=mod(b,m)，则称“a与b关于模m有相同的余数（同余）”，记做a≡b(mod m)（这就是同余式）。举个例子，mod(13,4)=1，mod(1,4)=1，则读成13和1关于模4同余，记做13≡1(mod 4)。当然，同余具有对称性，上式还可以写成1=13(mod 4)。
a≡b(mod m)的一个充要条件是a=b+mt，t是整数，比如13=1+3*4。a=b+mt可以写成(a-b)/m=t，即m能整除(a-b)。
这些就够了。更多基础性的介绍，可以参考《同余（数据基础）》。
1.乘同余法
同余法是一大类方法的统称，包括加同余法、乘同余法等。因为这些方法中的迭代公式都可以写成上面我们见过的同余式形式，故统称同余法。常用的就是下面的乘同余法(Multiplicative Congruential Generator.)。符号不好敲，做些约定，如R(i)就是R加一个下标i。
乘同余法随机数生成器的同余形式如下：R(i+1)=a*R(i) (mod m)。这个迭代式可以写成更直观的形式，R(i+1)=mod[a*R(i),m]，其中初值R(0)称为随机数种子。因为mod(x,m)总是等于0到m-1的一个整数，所以最后把R(i+1)这个随机序列都除以m，就可以得到在[0,1]上均匀分布的随机数。下面用电子表格演示一遍，假设随机数种子R(0)=1，a=4，m=13：
1
2
3
4
A
i
0
1
2
B
R(i)
1
=D2
=D3
C
a*R(i)
=B2*4
=B3*4
=B4*4
D E
mod[a*R(i),m] mod[a*R(i),m]/m
=mod(C2,13) =D2/13
=mod(C3,13) =D3/13
=mod(C4,13) =D4/13
上面显示的是公式（Excel中，公式与计算结果的转换，用快捷键ctrl+~实现），结果看起来是这样的，E列就是我们想要的在[0,1]上均匀分布的随机数：
1
2
3
4
A
i
0
1
2
B
R(i)
1
4
3
C
a*R(i)
4
16
12
D E
mod[a*R(i),m] mod[a*R(i),m]/m
4 0.307692308
3 0.230769231
12 0.923076923
上表中，a*R(1)=16，R(2)=3，m=13，一个同余式就是R(2)=a*R(1) (mod m)，3=16(mod 13)，或者说，16-3能够被13整除——同余法就是这意思了。说“乘同余法”，要点在于a*R(i)是乘法形式。在SAS系统中，a=397,204,094，m =
2^31-1=2147483647（一个素数），随机种子R(0)可以取1到m-1之间的任何整数。
2.伪随机数的检验
现在内置于各大软件包的随机数生成器都经过了彻底的检验，我们当然可以安心地使用这个黑箱。或则我们也可以多了解些。一个好的“伪”随机数列，应该看起来就跟从[0,1]上均匀分布的随机数列中随机抽取出来的一样。两个比较直观的检验有：
1. 均匀性检验——所有的数是否均匀地分布在[0,1]区间；
2. 独立性（或不相关）检验——序列中的数不存在序列相关，每个数都跟其前后出现的数独立无关。一个例子是如0.1、0.2、0.3、0.4、0.5，这里每一个相继的数都正好比它前面的一个数大0.1，这样这个数列就似乎有了某种“格局”。
具体的检验方法就略之了，更多请参见徐钟济（1985）。
3.生成其他概率分布的随机数
前面提到过，我们把[0,1]上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都是借助于随机数来实现的。对其他连续型分布，（累积）分布函数为F(x)，r是在[0,1]上均匀分布的随机变量，另r=F(x)，解出的x就是该连续型分布的随机抽样。由于x可以写成r的反函数，该变换就称作“逆变换”(Inverse Transformation method)。对离散型分布，思路类似，不过繁琐些，具体可以见徐钟济（1985）、Ross（2006）和Evans等（2001）。大多数相关软件包都会提供各种分布的随机数生成函数，知道有这回事就可以。最重要的，是我们陈述一类问题时，知道需要用哪种概率分布来描述Evans等（2001）：
常用的连续分布
1. 均匀分布(Uniform Distribution)描述在某最小值和最大值之间所有结果等可能的随机变量的特性。
2. 正态分布(Normal Distribution)是对称的，具有中位数等于均值的性质，就是我们熟悉的钟形形状。各种类型的误差常常是正态分布的。最后，作为中心极限定理的推论，大量具有任意分布的随机变量的平均数也呈正态分布。
3. 三角形分布(Triangular Distribution)由三个参数来定义，最大值、最小值和最可能值。临近最可能值的结果比那些位于端点的结果有更大的出现机会。
4. 指数分布(Exponential Distribution)常用来构建在时间上随机重现的事件的模型，如顾客到达服务系统的时间间隔、机器元件失效前的工作时间等。它的主要性质是“无记忆性”(Memoryless) ，即当前时间对未来结果没有影响。例如，不论机器原先已经运转了多长时间，它继续运转直至出现故障的时间间隔总有同样的分布。
5. 对数正态分布(Lognormal Distribution)：若随机变量x的自然对数是正态的，则x就服从对数正态分布。对数正态分布的常见例子是股票价格。
常用的离散分布
1. 贝努利分布(Bernoulli Distribution)描述的是只有两个以常数概率出现的可能结果的随机变量的特性；
2. 二项分布(Binomial Distribution)给出每次试验成功概率为p的n次独立重复贝努利实验的模型；
3. 泊松分布(Poisson Distribution)用于建立某种度量单位内发生次数的模型的一种离散分布，如某时段内某事件发生的次数。
其他有用的分布如伽玛分布(Gamma Distribution)、威布尔分布(Weibull
Distribution)、贝塔分布(Beta Distribution)略之。
4.下期预告
上次落了些东西。说到蒙卡与C++，在数量金融领域，不能不提的一本书就是Mark
Joshi的C++ Design Patterns and Derivatives Pricing (Cambridge Uni. Press,
2004)，面试中一定要说自己读过的，这样人家以为你至少会用C++做一个蒙特卡罗模拟。这个系列只讲SAS，下次先讲如何用SAS生成随机数，然后就是具体的项目。
SAS随机数函数及CALL子程序
**************************************************************************************************************************
《SAS和蒙特卡罗模拟（1）：开篇》——简介，通过例子建立起蒙卡的直观概念；参考软件包及书目
《SAS和蒙特卡罗模拟（2）：随机数基础》——随机数入门；参考书目
**************************************************************************************************************************************************
一、SAS随机数函数和CALL子程序
SAS系统产生随机数，两种方式，利用SAS函数(Functions)或者CALL子程序(CALL Routines)，它们的语法格式是（具体的区别容后讨论）：
方 代
式 码 明
函数 var=name(seed,
CALL子程序 call 同上,记得seed=0, ±1,±2, , ± (2**31-2)
name(seed,
SAS可用的随机数函数列表如下（可以参见SAS Help and Documentation-SAS
Products-Base SAS-SAS Language Dictionary-Functions and CALL
Routines-Functions and CALL Routines by Category）：
分 SAS函数
布
二项分布(Binomial) ranBin(seed,n,p)
指数分布ranExp(seed)
(Exponential)
正态分布(Normal) ranNor(seed),or
normal(seed)
说明
n为独立实验的次数，p为成功概率
ranNor和normal是同质的，但normal没有相对应的CALL子程序
泊松分布(Poisson) ranPoi(seed,m) m为均值
均匀分布(Uniform) ranUni(seed),or ranUni和uniform是uniform(seed) 同质的，但uniform没有相对应的CALL子程序
柯西分布(Cauchy) ranCau(seed)
伽玛分布(Gamma) ranGam(seed,a) a>0为形状参数
由分布律表格决定的ranTbl(seed,p1,p2,...pn) ∑p=1
离散分布(tabled
probability
distribution)
三角分布ranTri(seed,h) h为斜边（最可能值）
(Triangular)
上面的随机函数，除了normal和uniform，都可以由直接相应的CALL子程序调用。
二、SAS随机数函数和CALL子程序：比较
用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。这话费解，一个例子先，创建两个随机数变量，各包含3个记录，其中x1的种子为123，x2的种子为456：
data ranuni(drop=i);
retain seed1 123 seed2 456;
do i=1 to 3;
x1=ranuni(seed1);
x2=ranuni(seed2);
output;
end;
run;
proc print data=ranuni;run;
结果为：
Obs seed1 seed2 x1 x2
1 123 456 0.75040 0.32091
2 123 456 0.17839 0.90603
3 123 456 0.35712 0.22111
好像没什么异样。我们把上面的x1增加为6个记录：
data ranuni2(drop=i);
retain seed1 123;
do i=1 to 6;
x1=ranuni(seed1);
output;
end;
run;
proc print data=ranuni2;run;
结果如下，把上下用红色标出的数字对照看一看：
Obs seed1 x1
1
2
3
4
5
6
123
123
123
123
123
123
0.75040
0.32091
0.17839
0.90603
0.35712
0.22111
什么意思？在第一段代码中，x2的种子根本不起作用，把x2的记录安插到x1里，看起来就是用种子123产生的随机数列加长了而已。x2的第一个值并不是由种子456产生的，而是产生第一个x1后的下一个值，x1、x2其实属于同一个随机数列，尽管x2的种子被指定为456，而x1的被指定为123。现在就可以重复上面的一句话：用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。
用CALL子程序就能够同时产生独立的随机数列。
data ranuni3(drop=i);
retain seed3 123 seed4 456;
do i=1 to 3;
call ranuni(seed3,x3);
call ranuni(seed4,x4);
output;
end;
run;
proc print data=ranuni3;run;
结果如下：
Obs seed3 seed4 x3
x4
1 1611463328 736440516 0.75040 0.34293
2 689153326 774069794 0.32091 0.36045
3 383088854 686944750 0.17839 0.31988
以上x3就是x1。x1和x3的初始种子都是123，但x3那个结果显示的种子是当前种子值。要在SAS随机数函数语句中显示随机种子的当前值，可以由以下代码给出：
data ranuni4(drop=i);
retain seed1 123;
do i=1 to 3;
x1=ranuni(seed1);
seed=x1*(2**31-1);
output;
end;
run;
proc print data=ranuni4;
var seed1 seed x1;
run;
结果如下，可以跟上面由CALL子程序得出的结果对照：
Obs seed1 seed x1
1 123 1611463328 0.75040
2 123 689153326 0.32091
3 123 383088854 0.17839
---------参考资料---------
1. Xitao Fan, etc..Monte Carlo Studies: A Guide for Quantitative
Researchers. SAS Institute Inc.,2002
2. 朱世武《SAS编程技术与金融数据处理》，北京：清华大学出版社，2003
Matlab() 随机数生成方法：
第一种方法是用 random 语句，其一般形式为
y = random('分布的英文名',A1,A2,A3,m,n)，
表示生成 m 行 n 列的 m × n 个参数为 ( A1 , A2 , A3 ) 的该分布的随机数。例如:
(1) R = random('Normal',0,1,2,4): 生成期望为 0,标准差为 1 的(2 行 4 列)2× 4 个正态随机数
(2) R = random('Poisson',1:6,1,6): 依次生成参数为 1 到 6 的(1 行 6 列)6 个 Poisson 随机数
第二种方法是针对特殊的分布的语句：
一． 几何分布随机数 （下面的 P，m 都可以是矩阵）
R = geornd(P) （生成参数为 P 的几何随机数）
R = geornd(P,m) （生成参数为 P 的 × m 个几何随机数）
１
R = geornd(P,m,n) （生成参数为 P 的 m 行 n 列的 m × n 个几何随机数）
例如
(1) R = geornd(1./ 2.^(1:6)) ( 生成参数依次为 1/2,1/2^2,到 1/2^6 的 6 个几何随机数)
(2) R = geornd(0.01,[1 5]) (生成参数为 0.01 的（１行５列）5 个几何随机数).
二．Beta 分布随机数
R = betarnd(A,B) （生成参数为 A,B 的 Beta 随机数）
R = betarnd(A,B,m) （生成 × m 个数为 A,B 的 Beta 随机数）
１
R = betarnd(A,B,m,n) （生成 m 行 n 列的 m × n 个数为 A,B 的 Beta 随机数）.
三．正态随机数
R = normrnd(MU，SIGMA) （生成均值为 MU，标准差为 SIGMA 的正态随机数）
R = normrnd(MU，SIGMA,m) （生成 × m 个正态随机数）
１
R = normrnd(MU，SIGMA,m,n) （生成 m 行 n 列的 m × n 个正态随机数）
例如
(1) R = normrnd(0,1,[1 5]) 生成 5 个正态(0,1) 随机数
(2) R = normrnd([1 2 3;4 5 6],0.1,2,3) 生成期望依次为[1,2,3;4,5,6], 方差为 0.1 的 2× 3 个正态
随机数．
四．二项随机数：类似地有
R = binornd(N,P) R = binornd(N,P,m) R = binornd(N,p,m,n)
例如
n = 10:10:60; r1 = binornd(n,1./n) 或 r2 = binornd(n,1./n,[1 6]) （都生成参数分别为
1 1
), L, ( 60, ) 的６个二项随机数．
(10,
10 60
五．自由度为 V 的 χ 2 随机数：
R = chi2rnd(V) R = chi2rnd(V R = chi2rnd(V
,m) ,m,n)
六．期望为 MU 的指数随机数（即 Exp 随机数）：
1
MU
R = exprnd(MU) R = exprnd(MU,m) R = exprnd(MU,m,n)
七．自由度为 V1， V2 的 F 分布随机数：
R = frnd(V1，V2) R = frnd(V1， V2,m) R = frnd(V1，V2,m,n)
八． Γ ( A, λ ) 随机数：
R = gamrnd（A,lambda） R = gamrnd（A,lambda,m) R = gamrnd（A,lambda,m,n)
九．超几何分布随机数：
R = hygernd(N,K,M) R = hygernd(N,K,M,m) R = hygernd(N,K,M,m,n)
十．对数正态分布随机数
R = lognrnd(MU，SIGMA) R = lognrnd(MU，SIGMA,m) R = lognrnd(MU，SIGMA,m,n)
十一．负二项随机数：
R = nbinrnd(r,p) R = nbinrnd(r,p,m) R = nbinrnd(r,p,m,n)
十二．Poisson 随机数：
R = poissrnd(lambda) R = poissrnd(lambda,m) R = poissrnd(lambda,m,n)
例如，以下 3 种表达有相同的含义：lambda = 2; R = poissrnd(lambda,1,10)
（或 R = poissrnd(lambda,[1 10]) 或 R = poissrnd(lambda(ones(1,10)))
十三．Rayleigh 随机数：
R = raylrnd(B) R = raylrnd(B,m) R = raylrnd(B,m,n)
十四．V 个自由度的 t 分布的随机数：
R = trnd(V) R = trnd(V,m) R = trnd(V,m,n)
42
十五．离散的均匀随机数：
R = unidrnd(N) R = unidrnd(N,m) R = unidrnd(N,m,n)
十六．[A,B] 上均匀随机数
R = unifrnd(A,B) R = unifrnd(A,B,m) R = unifrnd(A,B,m,n)
例如 unifrnd(0,1:6)与 unifrnd(0,1:6,[1 6]) 都依次生成[0,1] 到[0,6]的６个均匀随机数．:
十七．Weibull 随机数
R = weibrnd(A,B) R = weibrnd(A,B,m) R = weibrnd(A,B,m,n)
发布者：admin，转转请注明出处：http://www.yc00.com/news/1708580544a1575610.html
评论列表（0条）