0728 統算 隨機分配生成 (離散部分)

[markoo (上行下笑阿..)]

同學，補上之前離散型分配 變數生成
那其他的就同學自行揣摩補上了:D
Bernoulli(p) :
P{X=1} = p,
P{X=0} = 1-p,

  X=0        X=1
|_______|_____________|
0  1-p         p      1
所以經由uniform(0,1)的落點可以得到一個樣本值X1=0或1, 直到生成Xn結束

double *rand_bernoulli(int n, double p, int seed) //Bernoulli(p)
{
        int i;
        double *ans;
        double *unif;
        ans = (double *) malloc (n*sizeof(double));
        unif = rand_uniform(n,0,1,seed);
        for (i=0;i<n;i++)
        {
                if (unif[i]<=p) ans[i]=1;
                else ans[i]=0;
        }
        free(unif);
        return ans;
}

Binomial(m,p) 可以直接由m個Bernoulli(p)加起來

double *rand_binomial(int n,int m, double p, int seed) //Binomial(m,p)
{
        int i,j;
        double *ans;
        double *bernoulli;
        ans = (double *) malloc (n*sizeof(double));
        bernoulli = rand_bernoulli(n*m,p,seed);
        for (i=0;i<n;i++)
        {
                ans[i]=0;
                for (j=0;j<m;j++) ans[i]+=bernoulli[j+i*m];
        }
        free(bernoulli);
        return ans;
}

這邊Geometric值由1,2,...,infinity.(X代表失敗次數)
  X=1      X=2        X=3 ....                        X=infinity
|________|_________|_________|_____|____|__|_||||||||||
0  p      (1-p)p    (1-p)^2*p                         1

一樣由uniform(0,1)去判斷他落在哪個區域就可以得到一個樣本

double *rand_geometric(int n, double p, int seed) // Geometric(p)
{
        int i,k;
        int check=0;
        double *unif;
        double *ans;
        double pleft = 0;
        double pright= 0;
        double temp;
        ans = (double *) malloc (n*sizeof(double));
        unif = rand_uniform(n,0,1,seed);
        for (i=0;i<n;i++)
        {
                pright =0;
                temp = p;
                k=0;
                do
                {
                        pleft =pright;
                        pright = pleft+ temp;
                        k++;
                        temp = (1-p)*temp;
                        //這邊temp是P{X=k+1}與P{X=k}之間的比例
                } while (unif[i]<= pleft || unif[i]>pright);
                ans[i]= k;
        }
        return ans;
}


--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 118.166.54.115