实用Matlab产生给定密度函数的随机数时，怎么会产生复数呢？！

秦时明星

求助各位GGJJ还有 老师啊！
使用的是谢中华老师的连续分布随机数产生的代码,
先是 M_file：
function y = crnd(pdffun, pdfdef, m, n)
%生成任意一元连续分布随机数
%   y = crnd(pdffun, pdfdef, m, n)，产生指定一元连续分布的随机数，m行n列。pdffun为密度
%   函数表达式，pdfdef为密度函数定义域。注意：pdfdef只能是有限区间，若密度函数定义域为无限区
%   间，应设成比较大的有限区间，例如[-10000,10000]
%
%   example:
%   pdffun = 'x*(x>=0 & x<1)+(2-x)*(x>=1 & x<2)';
%   x = crnd(pdffun, [0 2], 1000, 1);
%   [fp,xp] = ecdf(x);
%   ecdfhist(fp, xp, 20);
%   hold on
%   fplot(pdffun, [0 2], 'r')
%
%   Copyright 2009 - 2010 xiezhh.
%   $Revision: 1.0.0.0 $  $Date: 2009/12/17 12:10:00 $

fun = vectorize(['(' pdffun ')' '*x']);    % x*f(x)运算向量化
try
    xm = quadl(fun, min(pdfdef), max(pdfdef));    % 计算x的数学期望xm
catch
    xm = mean(pdfdef);    % 计算定义区间的平均值xm
end

pdffun = vectorize(['(' pdffun ')' '*x/x']);    % x*f(x)/x运算向量化
cdfrnd = rand(m*n, 1);    % 产生[0,1]上均匀分布随机数
y = zeros(m*n, 1);        % 产生0向量作为变量y的初值
options = optimset;       % 产生一个控制迭代过程的结构体变量
options.Display = 'off';  % 不显示中间迭代过程
% 通过循环计算指定一元连续分布的随机数
for i = 1:m*n
    funcdf = @(x)[quadl(pdffun, min(pdfdef), x) - cdfrnd(i)];
    y(i) = fsolve(funcdf, xm, options);
end
y=reshape(y,[m,n]);    % 把向量y变为矩阵
然后接下来是我的命令行：
k=-0.345119445;
a=11.67962773;
e=68.61830219;
pdffun=['1/' num2str(a) '/exp((1-(' num2str(k) '))*((-1/' num2str(k) '*log(1-(' num2str(k) ')*(x-' num2str(e) ')/' num2str(a) '))))/(1+exp(-((-1/' num2str(k) '*log(1-(' num2str(k) ')*(x-' num2str(e) ')/' num2str(a) ')))))^2'];
x = crnd(pdffun, [0 600], 10, 1);
运行结果为：
54.3107137607307 - 29.4301471730884i
54.3549771521278 - 29.1269826917345i
54.2471473724783 - 29.6516842154325i
54.3596474076828 - 28.9914163339997i
54.2680896805071 - 29.5886019413901i
54.3596393266516 - 28.9844750403979i
54.3468399374602 - 29.2164034806281i
54.2881735718014 - 29.5203171051013i
54.2334236943346 - 29.6895982304876i
54.3549393560571 - 28.8530002248770i

秦时明星

回复 1 # 秦时明星 的帖子

其中命令行中的密度函数是广义逻辑分布的密度函数：
f(x)=1/a/exp((1-k)*y)/(1+exp(-y))^2
其中，y=-1/k*ln(1-k*(x-e)/a)