1. 程式人生 > >[轉] 如何用matlab生成服從混合高斯分佈的隨機數?

[轉] 如何用matlab生成服從混合高斯分佈的隨機數?

M=10; %產生M行N列的隨機數矩陣
N=8;
miu1=1;%第一個分佈的引數
sigma1=2;%第一個分佈的引數
miu2=6;%第二個分佈的引數
sigma2=1;%第二個分佈的引數


R = 0.2*normrnd(miu1,sigma1,M,N)+0.8*normrnd(miu2,sigma2,M,N);


單點的概率全是0,那你取出來的隨機數算什麼?
若干個隨機數要滿足統計分佈,是要按區間統計的
另外我不知道你要做什麼就是了。
你如果想按一定的概率密度來產生隨機數,你最好用反函式法之類的來弄。

比如產生一個x.^2分佈的隨機數,不過這些要歸一化。

============================================
首先,我知道我的是錯的了。如下圖就可知
M=1000; %產生M行N列的隨機數矩陣
N=1;
miu1=1;%第一個分佈的引數
sigma1=2;%第一個分佈的引數
miu2=6;%第二個分佈的引數
sigma2=1;%第二個分佈的引數


R = 0.2*normrnd(miu1,sigma1,M,N)+0.8*normrnd(miu2,sigma2,M,N);


x=-5:0.001:15;
y1=normpdf(x,miu1,sigma1);
y2=normpdf(x,miu2,sigma2);
subplot(2,2,1);
plot(x,y1);
subplot(2,2,2);
plot(x,y2);
subplot(2,2,3);
y3=0.2*y1+0.8*y2;
plot(x,y3);

subplot(2,2,4)
dx=0.5;
xx=-5:dx:15;
yy=hist(R,xx);
yy=yy/M/dx;
plot(x,y3);
hold on
bar(xx,yy)

=======================================
正確做法,我還沒弄出來,繼續中。。。。


============================================
_____________________新的嘗試
下面的結果我覺得可能可以接受。
思路:基於反變換法
Matlab下面有
p=normpdf(x,miu,sigma)是求出x處的概率密度。
p=normcdf(x,miu,sigma)是求出X<x的累積概率密度(就是從負無窮大到x處的概率密度的積分)
我給定一個區間,這個區間外的概率我認為是0(這一點不夠嚴謹,理論上應當是從負無窮到正無窮)
我這裡取的是-10:15,其間我取了25000個點,求出這些點的累積概率值(兩個的加權和y3),記這個為F(x),根據反變換法,
F(x)=u,其中u是一個0到1的均勻隨機數。只要求出它的解x0,那麼x0就滿足所給定的概率密度分佈。這裡我用的是插值。用
(y3,x)來插值出u所在的位置

宣告,這裡有一些地方不夠嚴謹,嚴謹應當用解析的方法來做反變換。
%%%%%下面是程式


M=1000; %產生M行N列的隨機數矩陣
N=1;
miu1=1;%第一個分佈的引數
sigma1=2;%第一個分佈的引數
miu2=6;%第二個分佈的引數
sigma2=1;%第二個分佈的引數

x=-10:0.001:15;

y1=normpdf(x,miu1,sigma1);
y2=normpdf(x,miu2,sigma2);
y3=0.2*y1+0.8*y2;



y1=normcdf(x,miu1,sigma1);
y2=normcdf(x,miu2,sigma2);
y=0.2*y1+0.8*y2;
u=rand(N,M);
R=interp1(y,x,u,'linear');

dx=0.5;
xx=-10:dx:15;
yy=hist(R,xx);
yy=yy/M/dx;
bar(xx,yy)


hold on;
plot(x,y3,'r*')