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

2008057

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

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

ChaChing

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

xiezhh

matlab中好像没有现成的函数，给你一个我编的函数，代码如下：

1. function y = crnd(pdffun, pdfdef, m, n)
   
2. %生成任意一元连续分布随机数
   
3. %   y = crnd(pdffun, pdfdef, m, n)，产生指定一元连续分布的随机数，m行n列。pdffun为密度
   
4. %   函数表达式，pdfdef为密度函数定义域。注意：pdfdef只能是有限区间，若密度函数定义域为无限区
   
5. %   间，应设成比较大的有限区间，例如[-10000,10000]
   
6. %
   
7. %   example:
   
8. %   pdffun = 'x*(x>=0 & x<1)+(2-x)*(x>=1 & x<2)';
   
9. %   x = crnd(pdffun, [0 2], 1000, 1);
   
10. %   [fp,xp] = ecdf(x);
    
11. %   ecdfhist(fp, xp, 20);
    
12. %   hold on
    
13. %   fplot(pdffun, [0 2], 'r')
    
14. %
    
15. %   Copyright 2009 - 2010 xiezhh.
    
16. %   $Revision: 1.0.0.0 $  $Date: 2009/12/17 12:10:00 $
    
17. 
    
18. fun = vectorize(['(' pdffun ')' '*x']);    % x*f(x)运算向量化
    
19. try
    
20.     xm = quadl(fun, min(pdfdef), max(pdfdef));    % 计算x的数学期望xm
    
21. catch
    
22.     xm = mean(pdfdef);    % 计算定义区间的平均值xm
    
23. end
    
24. 
    
25. pdffun = vectorize(['(' pdffun ')' '*x/x']);    % x*f(x)/x运算向量化
    
26. cdfrnd = rand(m*n, 1);    % 产生[0,1]上均匀分布随机数
    
27. y = zeros(m*n, 1);        % 产生0向量作为变量y的初值
    
28. options = optimset;       % 产生一个控制迭代过程的结构体变量
    
29. options.Display = 'off';  % 不显示中间迭代过程
    
30. % 通过循环计算指定一元连续分布的随机数
    
31. for i = 1:m*n
    
32.     funcdf = @(x)[quadl(pdffun, min(pdfdef), x) - cdfrnd(i)];
    
33.     y(i) = fsolve(funcdf, xm, options);
    
34. end
    
35. y=reshape(y,[m,n]);    % 把向量y变为矩阵

复制代码

crnd函数能解决你的问题，只是该函数的效率还有待进一步提高，
这里用一个例子说明该函数的用法：

1. pdffun = '6*x*(1-x)';    % 密度函数表达式
   
2. % 调用crnd函数生成100个服从指定一元连续分布的随机数
   
3. x = crnd(pdffun, [0 1], 100, 1);
   
4. [fp,xp] = ecdf(x);    % 计算经验累积概率分布函数值
   
5. ecdfhist(fp,xp,20);   % 绘制频率直方图
   
6. hold on
   
7. fplot(pdffun, [0 1], 'r')    % 绘制真实密度函数曲线
   
8. xlabel('x');       % 为X轴加标签
   
9. ylabel('f(x)');    % 为Y轴加标签
   
10. legend('频率直方图', '密度函数曲线')    % 为图形加标注框

复制代码

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

2008057

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

ChaChing

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

2008057

我要求的概率密度函数是 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)[quadl(pdffun, min(pdfdef), x) - cdfrnd(i)] 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)[quadl(pdffun, min(pdfdef), x) - cdfrnd(i)] at 32
  In optim\private\trustnleqn at 210
  In fsolve at 295
  In crnd at 33
>> 这只是一部分，不明白怎么回事，楼主给的例子倒是能实现，请帮我看看这是怎么回事啊，产生的随机数数量应该没有限制吧，我用1000和10都试过了，结果都一样，郁闷啊。求教， 谢谢了

xiezhh

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

八戒

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

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

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

xiezhh

这样吧

1. 
   
2. >> aa=1;
   
3. >> pdffun =['6*x*(' num2str(aa) '-x)'];
   
4. >> x = crnd(pdffun, [0 1], 100, 1);
   

复制代码

八戒

非常感谢 谢老师 @
用上num2st命令，问题解决了@
您上面编的crnd命令也很好用，很好地生成的分布函数！！
但是运行时都会跳出类似以下的警告，不知为什么：（但可正常得到分布数）
> In inlineeval at 13
  In inline.feval at 34
  In quadl at 64
  In crnd>@(x)[quadl(pdffun, min(pdfdef), x) - cdfrnd(i)] 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

出现了分母为0的情况。

八戒

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

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

xiezhh

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

yanyang0425

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

ljelly

回复 14 # yanyang0425 的帖子

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

你的不能实现是什么，具体错误信息是什么