2008057 发表于 2010-1-28 15:17

已知概率密度函数的随机数?

用MATLAB生成相应分布的随机数，怎么办啊，已知概率密度函数为：exp(a+bx+cx.^2+dx.^3+ex.^5)可以吗，用哪个函数啊，请高手指教了，小妹不胜感激，不要再给我删了好不好啊

[ 本帖最后由 ChaChing 于 2010-1-29 00:16 编辑 ]

ChaChing 发表于 2010-1-29 00:14

回复 楼主 2008057 的帖子

基本上, 若不违反版内规定(重覆或相似...), 个人是不可能删的! 请放心

xiezhh 发表于 2010-1-29 11:00

matlab中好像没有现成的函数，给你一个我编的函数，代码如下：
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, , 1000, 1);
%    = ecdf(x);
%   ecdfhist(fp, xp, 20);
%   hold on
%   fplot(pdffun, , '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);    % 产生上均匀分布随机数
y = zeros(m*n, 1);      % 产生0向量作为变量y的初值
options = optimset;       % 产生一个控制迭代过程的结构体变量
options.Display = 'off';% 不显示中间迭代过程
% 通过循环计算指定一元连续分布的随机数
for i = 1:m*n
    funcdf = @(x);
    y(i) = fsolve(funcdf, xm, options);
end
y=reshape(y,);    % 把向量y变为矩阵
crnd函数能解决你的问题，只是该函数的效率还有待进一步提高，
这里用一个例子说明该函数的用法：
pdffun = '6*x*(1-x)';    % 密度函数表达式
% 调用crnd函数生成100个服从指定一元连续分布的随机数
x = crnd(pdffun, , 100, 1);
= ecdf(x);    % 计算经验累积概率分布函数值
ecdfhist(fp,xp,20);   % 绘制频率直方图
hold on
fplot(pdffun, , 'r')    % 绘制真实密度函数曲线
xlabel('x');       % 为X轴加标签
ylabel('f(x)');    % 为Y轴加标签
legend('频率直方图', '密度函数曲线')    % 为图形加标注框

这是我的新书《MATLAB统计分析与应用39个案例分析》中的一个例子，
该书即将出版，敬请关注！
效果图：

2008057 发表于 2010-1-29 20:11

回复 板凳 xiezhh 的帖子

多谢了，我按照你说的方法试了，总是显示错误：说crnd是未经定义的命令/函数，怎么回事啊。

ChaChing 发表于 2010-1-29 23:19

回复 地板 2008057 的帖子

LZ未将crnd.m存在有效路径(path)!?
Ref: 常见的程序出错问题整理 (eight), 3F
http://forum.vibunion.com/forum/thread-46001-1-1.html

2008057 发表于 2010-1-30 09:36

我要求的概率密度函数是 exp(134200.*x-9824.*(x.^2)+52.46.*(x.^3)-0.4468.*(x.^4)+0.08184.*(x.^5))；x的区间是【69.84，70.24】；    用上面的方法试，总是显示一大串出错信息：Warning: Infinite or Not-a-Number function value encountered.
> In quadl at 101
In crnd>@(x) at 32
In sfdnls at 90
In optim\private\trustnleqn at 108
In fsolve at 295
In crnd at 33
Warning: Infinite or Not-a-Number function value encountered.
> In quadl at 101
In crnd>@(x) at 32
In optim\private\trustnleqn at 210
In fsolve at 295
In crnd at 33
>> 这只是一部分，不明白怎么回事，楼主给的例子倒是能实现，请帮我看看这是怎么回事啊，产生的随机数数量应该没有限制吧，我用1000和10都试过了，结果都一样，郁闷啊。求教， 谢谢了

xiezhh 发表于 2010-1-30 16:58

你的函数：exp(134200.*x-9824.*(x.^2)+52.46.*(x.^3)-0.4468.*(x.^4)+0.08184.*(x.^5))根本就不是密度函数，最起码在定义区间上积分不是1.

八戒 发表于 2010-4-21 22:14

谢老师，您好，
我想请问一下，您的函数应用时，
概率密度函数中是否不能有‘已经被赋值的变量’
如：
aa=1;
pdffun = '6*x*(aa-x)';    % 密度函数表达式;

我运行的结果是，因为aa存在，所以出错中止；
请问谢老师是否有什么方法解决此问题？
先谢过了！

[ 本帖最后由 八戒 于 2010-4-21 22:16 编辑 ]

xiezhh 发表于 2010-4-22 10:23

这样吧

>> aa=1;
>> pdffun =['6*x*(' num2str(aa) '-x)'];
>> x = crnd(pdffun, , 100, 1);

八戒 发表于 2010-4-22 15:59

非常感谢 谢老师 @
用上num2st命令，问题解决了@
您上面编的crnd命令也很好用，很好地生成的分布函数！！
但是运行时都会跳出类似以下的警告，不知为什么：（但可正常得到分布数）
> In inlineeval at 13
In inline.feval at 34
In quadl at 64
In crnd>@(x) at 32
In sfdnls at 90
In optim\private\trustnleqn at 268
In fsolve at 295
In crnd at 33
In xiangweicha at 101
Warning: Divide by zero.

xiezhh 发表于 2010-4-24 08:57

出现了分母为0的情况。

八戒 发表于 2010-4-24 20:39

但是可以得出分布数，而且和概率密度的趋势线吻合的还不错，说明也是能正确地按概率密度进行分布~

这个出现分母为0应该不严重吧？
有什么方法解决此问题吗？

xiezhh 发表于 2010-4-25 17:07

影响不大，可以不用理会。加个warning off 命令可以不显示警告信息。

yanyang0425 发表于 2011-1-13 16:26

我怎么连老师给的那个例子也实现不了呢？总是提示18行有错误

ljelly 发表于 2011-1-17 15:40

回复 14 # yanyang0425 的帖子

谢老师给出的代码不要放在一起
上面是M-函数
下面是M-script，可以直接运行，并调用了上面的M函数
要存成两个文件，用后才调用前者

你的不能实现是什么，具体错误信息是什么
页: [1] 2
查看完整版本: 已知概率密度函数的随机数?