[轉] 如何用matlab生成服從混合高斯分佈的隨機數？
阿新 • • 發佈：2019-02-18
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*')
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*')