Mengonversi Distribusi Seragam ke Distribusi Normal

Bagaimana cara mengubah distribusi seragam (seperti yang dihasilkan oleh kebanyakan generator bilangan acak, misalnya antara 0,0 dan 1,0) menjadi distribusi normal? Bagaimana jika saya menginginkan mean dan deviasi standar yang saya pilih?

The Ziggurat algoritma ini cukup efisien untuk ini, meskipun Box-Muller transformasi lebih mudah untuk menerapkan dari awal (dan tidak gila lambat).

Ada banyak metode:

• Jangan tidak menggunakan Box Muller. Apalagi jika Anda menggambar banyak angka gaussian. Box Muller menghasilkan hasil yang dijepit antara -6 dan 6 (dengan asumsi presisi ganda. Hal-hal memburuk dengan pelampung.). Dan itu benar-benar kurang efisien dibandingkan metode lain yang tersedia.
• Ziggurat baik-baik saja, tetapi membutuhkan pencarian tabel (dan beberapa penyesuaian khusus platform karena masalah ukuran cache)
• Rasio seragam adalah favorit saya, hanya beberapa penambahan / perkalian dan log 1/50 dari waktu (misalnya lihat di sana ).
• Membalik CDF itu efisien (dan diabaikan, mengapa?), Anda memiliki implementasi yang cepat tersedia jika Anda mencari di google. Ini wajib untuk nomor Quasi-Random.

Mengubah distribusi fungsi apa pun ke fungsi lain melibatkan penggunaan kebalikan dari fungsi yang Anda inginkan.

Dengan kata lain, jika Anda membidik fungsi probabilitas tertentu p (x), Anda mendapatkan distribusi dengan mengintegrasikannya -> d (x) = integral (p (x)) dan menggunakan inversnya: Inv (d (x)) . Sekarang gunakan fungsi probabilitas acak (yang memiliki distribusi seragam) dan berikan nilai hasil melalui fungsi Inv (d (x)). Anda harus mendapatkan nilai acak yang diberikan dengan distribusi sesuai dengan fungsi yang Anda pilih.

Ini adalah pendekatan matematika umum - dengan menggunakannya Anda sekarang dapat memilih fungsi probabilitas atau distribusi yang Anda miliki selama memiliki pendekatan invers atau invers yang baik.

Semoga ini membantu dan terima kasih atas komentar kecil tentang menggunakan distribusi dan bukan probabilitas itu sendiri.

Berikut adalah implementasi javascript menggunakan bentuk kutub dari transformasi Box-Muller.

/*
 * Returns member of set with a given mean and standard deviation
 * mean: mean
 * standard deviation: std_dev 
 */
function createMemberInNormalDistribution(mean,std_dev){
    return mean + (gaussRandom()*std_dev);
}

/*
 * Returns random number in normal distribution centering on 0.
 * ~95% of numbers returned should fall between -2 and 2
 * ie within two standard deviations
 */
function gaussRandom() {
    var u = 2*Math.random()-1;
    var v = 2*Math.random()-1;
    var r = u*u + v*v;
    /*if outside interval [0,1] start over*/
    if(r == 0 || r >= 1) return gaussRandom();

    var c = Math.sqrt(-2*Math.log(r)/r);
    return u*c;

    /* todo: optimize this algorithm by caching (v*c) 
     * and returning next time gaussRandom() is called.
     * left out for simplicity */
}

Gunakan teorema batas pusat entri wikipedia mathworld entry untuk keuntungan Anda.

Hasilkan n dari bilangan terdistribusi seragam, jumlahkan, kurangi n * 0,5 dan Anda memiliki output dari distribusi yang kira-kira normal dengan mean sama dengan 0 dan varians sama dengan (1/12) * (1/sqrt(N))(lihat wikipedia tentang distribusi seragam untuk yang terakhir)

n = 10 memberi Anda sesuatu yang setengah layak dengan cepat. Jika Anda menginginkan sesuatu yang lebih dari setengah yang layak, gunakan solusi tylers (seperti yang dicatat dalam entri wikipedia pada distribusi normal )

Saya akan menggunakan Box-Muller. Dua hal tentang ini:

1. Anda berakhir dengan dua nilai per iterasi
   Biasanya, Anda menyimpan satu nilai ke cache dan mengembalikan yang lain. Pada panggilan berikutnya untuk sampel, Anda mengembalikan nilai yang di-cache.
2. Box-Muller memberikan Z-score
   Anda harus menskalakan Z-score dengan standar deviasi dan menambahkan mean untuk mendapatkan nilai penuh dalam distribusi normal.

Dimana R1, R2 adalah nomor seragam acak:

DISTRIBUSI NORMAL, dengan SD 1: sqrt (-2 * log (R1)) * cos (2 * pi * R2)

Ini tepat ... tidak perlu melakukan semua putaran lambat itu!

Tampaknya luar biasa bahwa saya dapat menambahkan sesuatu ke dalamnya setelah delapan tahun, tetapi untuk kasus Java saya ingin mengarahkan pembaca ke metode Random.nextGaussian () , yang menghasilkan distribusi Gaussian dengan mean 0.0 dan deviasi standar 1.0 untuk Anda.

Penjumlahan dan / atau perkalian sederhana akan mengubah mean dan deviasi standar sesuai kebutuhan Anda.

Standar Python perpustakaan modul random memiliki apa yang Anda inginkan:

normalvariate (mu, sigma)
Distribusi normal. mu adalah mean, dan sigma adalah standar deviasi.

Untuk algoritme itu sendiri, lihat fungsi di random.py di pustaka Python.

The pengguna masuk sini

Ini adalah implementasi JavaScript saya dari Algorithm P ( metode Polar untuk deviasi normal ) dari Bagian 3.4.1 dari buku Donald Knuth The Art of Computer Programming :

function normal_random(mean,stddev)
{
    var V1
    var V2
    var S
    do{
        var U1 = Math.random() // return uniform distributed in [0,1[
        var U2 = Math.random()
        V1 = 2*U1-1
        V2 = 2*U2-1
        S = V1*V1+V2*V2
    }while(S >= 1)
    if(S===0) return 0
    return mean+stddev*(V1*Math.sqrt(-2*Math.log(S)/S))
}

Aku hal yang Anda harus mencoba ini di EXCEL: =norminv(rand();0;1). Ini akan menghasilkan angka acak yang seharusnya terdistribusi normal dengan mean nol dan menyatukan varians. "0" dapat diberikan dengan nilai apa pun, sehingga angka-angka tersebut akan menjadi rata-rata yang diinginkan, dan dengan mengubah "1", Anda akan mendapatkan varians yang sama dengan kuadrat input Anda.

Misalnya: =norminv(rand();50;3)akan menghasilkan bilangan yang terdistribusi normal dengan MEAN = 50 VARIANCE = 9.

T Bagaimana saya dapat mengubah distribusi seragam (seperti yang dihasilkan oleh kebanyakan generator bilangan acak, misalnya antara 0,0 dan 1,0) menjadi distribusi normal?

1. Untuk implementasi perangkat lunak saya tahu beberapa nama generator acak yang memberi Anda urutan acak seragam semu di [0,1] (Mersenne Twister, Linear Congruate Generator). Sebut saja U (x)

2. Ada bidang matematika yang disebut teori probabilitas. Hal pertama: Jika Anda ingin memodelkan rv dengan distribusi integral F maka Anda dapat mencoba mengevaluasi F ^ -1 (U (x)). Dalam teori pr. Terbukti bahwa rv tersebut memiliki distribusi integral F.

3. Langkah 2 dapat diterapkan untuk menghasilkan rv ~ F tanpa menggunakan metode penghitungan apa pun ketika F ^ -1 dapat diturunkan secara analitik tanpa masalah. (mis. exp.distribution)

4. Untuk memodelkan distribusi normal, Anda dapat menghitung y1 * cos (y2), dengan y1 ~ seragam dalam [0,2pi]. dan y2 adalah distribusi relei.

T: Bagaimana jika saya menginginkan mean dan deviasi standar yang saya pilih?

Anda dapat menghitung sigma * N (0,1) + m.

Dapat dibuktikan bahwa pergeseran dan penskalaan tersebut mengarah ke N (m, sigma)

Ini adalah implementasi Matlab menggunakan bentuk kutub dari transformasi Box-Muller :

Fungsi randn_box_muller.m:

function [values] = randn_box_muller(n, mean, std_dev)
    if nargin == 1
       mean = 0;
       std_dev = 1;
    end

    r = gaussRandomN(n);
    values = r.*std_dev - mean;
end

function [values] = gaussRandomN(n)
    [u, v, r] = gaussRandomNValid(n);

    c = sqrt(-2*log(r)./r);
    values = u.*c;
end

function [u, v, r] = gaussRandomNValid(n)
    r = zeros(n, 1);
    u = zeros(n, 1);
    v = zeros(n, 1);

    filter = r==0 | r>=1;

    % if outside interval [0,1] start over
    while n ~= 0
        u(filter) = 2*rand(n, 1)-1;
        v(filter) = 2*rand(n, 1)-1;
        r(filter) = u(filter).*u(filter) + v(filter).*v(filter);

        filter = r==0 | r>=1;
        n = size(r(filter),1);
    end
end

Dan memohon histfit(randn_box_muller(10000000),100);ini adalah hasilnya:

Jelas ini sangat tidak efisien dibandingkan dengan randn bawaan Matlab .

Saya memiliki kode berikut yang mungkin bisa membantu:

set.seed(123)
n <- 1000
u <- runif(n) #creates U
x <- -log(u)
y <- runif(n, max=u*sqrt((2*exp(1))/pi)) #create Y
z <- ifelse (y < dnorm(x)/2, -x, NA)
z <- ifelse ((y > dnorm(x)/2) & (y < dnorm(x)), x, z)
z <- z[!is.na(z)]

Juga lebih mudah menggunakan fungsi rnorm () yang diimplementasikan karena lebih cepat daripada menulis generator bilangan acak untuk distribusi normal. Lihat kode berikut sebagai bukti

n <- length(z)
t0 <- Sys.time()
z <- rnorm(n)
t1 <- Sys.time()
t1-t0

function distRandom(){
  do{
    x=random(DISTRIBUTION_DOMAIN);
  }while(random(DISTRIBUTION_RANGE)>=distributionFunction(x));
  return x;
}