《量化投资:以MATLAB为工具》

MATLAB技术论坛

 找回密码
 注册账号
查看: 28411|回复: 22
收起左侧

[提问] 已知某概率密度函数,如何产生一服从该分布的随机数

  [复制链接]
发表于 2011-4-15 10:29:59 | 显示全部楼层 |阅读模式
在做毕设的过程中遇到了一个问题,一概率密度函数为:
R=0.4*(1/(q*sqrt(2*pi)))*exp(-(x-u1)^2/2*q^2)+0.6*(1/(2*b))*exp(-abs(x-u2)/b);
(其中,前者是正态分布,后者是拉普拉斯分布,q、u1、u2、b均已知)
要产生一个服从该分布的随机数。
我查了好多资料了,一直都编不出来。哪位大侠能帮帮忙,万分感谢!
发表于 2012-4-30 11:29:53 | 显示全部楼层
我自己想了个办法,给你一个具体的例子吧
要求生成满足概率分布为f(x)=sin(x)^2/(x)^2*1/pi的随机数列

while i<n    % n为你要生成的随机数个数 %
      x1=2*(rand(1)-0.5)*8;  % 这个地方根据自己函数需要设定,但生成的 x1 必须在定义域内均匀分布  %
      yy=sin(x1)^2/(x1)^2*1/pi;  %  这是一个具体的概率密度函数的例子,你可以代入自己的概率密度函数 %
      if  subs(yy)>rand(1);   %  这里是这个办法的关键 %
          x(i,1)=x1;
          i=i+1;
      else
          i=i;
      end
end

得到的 x 即是一个满足给定的概率密度分布的一维数列
回复 支持 1 反对 0

使用道具 举报

发表于 2015-12-17 14:31:29 | 显示全部楼层
本帖最后由 yoki 于 2015-12-17 14:39 编辑

请问一下,我利用楼上所说的crnd函数,想生成下面概率密度函数y对应的随机数,可是随机数矩阵x的最大值达到274,超过[0.5 100]定义域范围,是为什么呢?编写的程序如下:
clc;
clear all;
y='0.1286*exp(-((x-99.71)/0.8997)^2) + 0.006678*exp(-((x-78.63)/35.77)^2)';
x=crnd(y,[0.5,100],1000,1);
[fp,xp] = ecdf(x);    % 计算经验累积概率分布函数值
ecdfhist(fp,xp,20);   % 绘制频率直方图
hold on
fplot(y,[0.5,100], 'r')    % 绘制真实密度函数曲线
xlabel('x');       % 为X轴加标签
ylabel('f(x)');    % 为Y轴加标签
回复 支持 反对

使用道具 举报

发表于 2016-5-22 23:23:17 | 显示全部楼层
xiezhh 发表于 2011-4-15 10:39
看看《MATLAB统计分析与应用:40个案例分析》,里面有相关案例。

谢老师您好,想请教您一个问题,就是按照您书上的例子可以获得一元函数的随机数,现在我已知一个复数的概率密度函数 fun=||v(x)||^2/64 其中x是一个复数 ,即是一个平面内的点坐标,这种情况我应该如何根据这个密度函数求得随机数呀?望谢老师能够回复!万分感激!!
回复 支持 反对

使用道具 举报

发表于 2011-4-15 10:39:36 | 显示全部楼层
看看《MATLAB统计分析与应用:40个案例分析》,里面有相关案例。
 楼主| 发表于 2011-4-15 10:50:19 | 显示全部楼层
 楼主| 发表于 2011-4-17 11:18:13 | 显示全部楼层
谢老师:
我想用您说的这个程序来求解问题。您编写的程序的M文件就是下面的这个crnd函数。我将我自己编的概率密度函数放入后出现了好多问题,不知道该如何解决,不知道您能不能帮忙看一下,非常感谢!

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               %try模块
    xm = quadl(fun, min(pdfdef), max(pdfdef));    % 计算x的数学期望xm,对fun进行积分,后面两个参数分别为上下限
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变为矩阵

在commend window中,有如下程序:

>> x=-10:0.1:10;
u1=0;
u2=0;
q=5;
b=q/sqrt(2);
pdffun='0.4*(1/(q*sqrt(2*pi)))*exp(-(x-u1).^2/2*q^2)+0.6*(1/(2*b))*exp(-abs(x-u2)/b)';
% 调用crnd函数生成100个服从指定一元连续分布的随机数
x = crnd(pdffun, [-10 10], 100, 1);
[fp,xp] = ecdf(x);    % 计算经验累积概率分布函数值
ecdfhist(fp,xp,20);   % 绘制频率直方图
hold on
fplot(pdffun, [-10 10], 'r')    % 绘制真实密度函数曲线
xlabel('x');       % 为X轴加标签
ylabel('f(x)');    % 为Y轴加标签
legend('频率直方图', '密度函数曲线')    % 为图形加标注框

运行后出现了很多问题:

??? Error using ==> inline.feval
Not enough inputs to inline function.

Error in ==> quadl at 64
y = feval(f,x,varargin{:}); y = y(:).';

Error in ==> crnd>@(x)[quadl(pdffun, min(pdfdef), x) - cdfrnd(i)] at 32
    funcdf = @(x)[quadl(pdffun, min(pdfdef), x) - cdfrnd(i)];

Error in ==> fsolve at 180
        fuser = feval(funfcn{3},x,varargin{:});

Error in ==> crnd at 33
    y(i) = fsolve(funcdf, xm, options);
发表于 2011-4-18 22:34:54 | 显示全部楼层
把参数的值直接代入密度函数表达式中即可。
a = 1; b = 2; fun = 'a+b*x' 和fun = '1+2*x'不是一回事。
 楼主| 发表于 2011-4-19 09:59:40 | 显示全部楼层
谢谢老师,这个问题解决了!
 楼主| 发表于 2011-4-25 15:27:55 | 显示全部楼层
回复 xiezhh 的帖子

       老师,我还有一个问题。在您的程序中,有这样一句:
fun = vectorize(['(' pdffun ')' '*x']);    % x*f(x)运算向量化
       它是否要求矩阵x和f(x)的大小相同?
       我用这个产生随机数的函数时,概率密度函数是由另两个概率密度函数做卷积得到的,我编写的程序如下:
x=-15:0.1:15;
y=0.1355*exp(-(x-0.103).^2/2.776)+0.2408*exp(-abs(x-0.103)/0.8027);
f1=0.0892*exp(x.^2/40);
crnd('conv(f1,y)', [-15 15], 1, 1);
这样一来,矩阵x和f(x)的大小不同。是否在这种情况下就不能用您编写的那个程序crnd.m了呢?
发表于 2011-4-25 20:32:53 | 显示全部楼层
crnd函数支持密度函数表达式作为第一个输入,不支持密度函数值作为输入。
 楼主| 发表于 2011-4-26 09:27:36 | 显示全部楼层
本帖最后由 chengyameng123 于 2011-4-26 09:44 编辑

谢谢老师,我明白啦!
那像我这样就不能用这个函数了,是吗?
不知道老师能不能推荐一下其他函数来解决这个问题:
x=-15:0.1:15;
y=0.1355*exp(-(x-0.103).^2/2.776)+0.2408*exp(-abs(x-0.103)/0.8027);
f1=0.0892*exp(x.^2/40);
z=conv(f1,y);

要产生一个服从z所示分布的随机数。
老师回答了我的好多问题了,真的很谢谢你!
发表于 2011-6-14 16:11:23 | 显示全部楼层
正好在找这个,学习一下
您需要登录后才可以回帖 登录 | 注册账号

本版积分规则

QQ|网站地图|MATLAB技术论坛|Simulink仿真论坛 ( 蜀ICP备19014457号 

GMT+8, 2020-11-27 11:17 , Processed in 0.059241 second(s), 13 queries , Gzip On, MemCached On.

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表