Menghasilkan angka acak dari "distribusi seragam miring" dari teori matematika

9

Untuk beberapa tujuan, saya perlu menghasilkan angka acak (data) dari distribusi "seragam miring". "Kemiringan" dari distribusi ini dapat bervariasi dalam beberapa interval yang masuk akal, dan kemudian distribusi saya harus berubah dari seragam menjadi segitiga berdasarkan pada kemiringan. Inilah derivasi saya:

masukkan deskripsi gambar di sini

Mari kita membuatnya sederhana dan menghasilkan data dari hingga (biru, merah adalah distribusi seragam). Untuk mendapatkan fungsi kerapatan probabilitas garis biru, saya hanya perlu persamaan garis itu. Jadi:B0B

f(x)=tg(φ)x+Y(0)

dan karena (gambar):

tg(φ)=1/BY(0)B/2Y(0)=1Btg(φ)B2

Kami memiliki itu:

f(x)=tg(φ)x+(1Btg(φ)B2)

Karena adalah PDF, CDF sama dengan:f(x)

F(x)=tg(φ)x22+x(1Btg(φ)B2)

Sekarang mari kita membuat generator data. Idenya adalah, bahwa jika saya akan memperbaiki , angka acak dapat dihitung jika saya akan mendapatkan angka dari dari distribusi seragam seperti dijelaskan di sini . Jadi, jika saya memerlukan 100 angka acak dari distribusi saya dengan fixed , maka untuk setiap dari distribusi seragam ada dari "distribusi miring", dan dapat dihitung sebagai:( 0 , 1 ) x i xφ,Bx(0,1)φ,Bti(0,1)xix

tg(φ)xi22+xi(1Btg(φ)B2)ti=0

Dari teori ini saya membuat kode dalam Python yang terlihat seperti:

import numpy as np
import math
import random
def tan_choice():
    x = random.uniform(-math.pi/3, math.pi/3)
    tan = math.tan(x)
    return tan

def rand_shape_unif(N, B, tg_fi):
    res = []
    n = 0
    while N > n:
        c = random.uniform(0,1)
        a = tg_fi/2
        b = 1/B - (tg_fi*B)/2
        quadratic = np.poly1d([a,b,-c])
        rots = quadratic.roots
        rot = rots[(rots.imag == 0) & (rots.real >= 0) & (rots.real <= B)].real
        rot = float(rot)
        res.append(rot)
        n += 1
    return res

def rand_numb(N_, B_):
    tan_ = tan_choice()
    res = rand_shape_unif(N_, B_, tan_)
    return res

Tetapi angka-angka yang dihasilkan dari rand_numbsangat dekat dengan nol atau ke B (yang saya setel 25). Tidak ada perbedaan, ketika saya menghasilkan 100 angka, semuanya hampir 25 atau semuanya hampir nol. Dalam sekali jalan:

num = rand_numb(100, 25)
numb
Out[140]: 
[0.1063241766836174,
 0.011086243095907753,
 0.05690217839063588,
 0.08551031241199764,
 0.03411227661295121,
 0.10927087752739746,
 0.1173334720516189,
 0.14160616846114774,
 0.020124543145515768,
 0.10794924067959207]

Jadi pasti ada sesuatu yang sangat salah dalam kode saya. Adakah yang bisa membantu saya dengan derivasi atau kode saya? Saya tergila-gila dengan ini sekarang, saya tidak dapat melihat kesalahan apa pun. Saya kira kode R akan memberi saya hasil yang serupa.

Robert
sumber
2
Jika Anda hanya perlu membuat angka acak, Anda tidak perlu menghitung distribusi sama sekali. Cukup lempar anak panah ke gambar Anda dan pertahankan koordinat x mereka, tetapi ketika panah mendarat di segitiga kiri berlabel " ", ubah koordinat x-nya dari ke . Misalnya, berikan nilai apa saja ke dan (parameter nyata yang, ketika diberi nilai antara dan , menghasilkan distribusi Anda) dan setel menjadi jumlah nilai acak yang Anda butuhkan. Ini kode:x B - x - 1 1ϕxBxBtheta11nRx<-runif(n,-1,1);x<-(ifelse(runif(n,-1,1)>theta*x,-x,x)+1)*(B/2)
whuber

Jawaban:

9

Deriviasi Anda baik-baik saja. Perhatikan bahwa untuk mendapatkan kepadatan positif pada , Anda harus membatasi Dalam kode Anda sehingga Anda harus mengambil antara , di situlah kode Anda gagal.(0,B)

B2tanϕ<2.
B=25ϕ±tan12625

Anda bisa (dan harus) menghindari menggunakan solver kuadrat, dan kemudian pilih akar antara 0 dan . Persamaan polinomial kuadrat dalam harus dipecahkan adalah dengan Dengan konstruksi dan ; juga meningkat pada .Bx

F(x)=t
F(x)=12tanϕx2+(1BB2tanϕ)x.
F(0)=0F(B)=1F(0,B)

Dari sini mudah untuk melihat bahwa jika , bagian parabola yang kami minati adalah bagian dari sisi kanan parabola, dan akar untuk dijaga adalah yang tertinggi dari dua akar, yang adalah Sebaliknya, jika , parabola terbalik, dan kami tertarik pada sisi kiri bagian. Akar untuk menjaga adalah yang terendah. Memperhatikan tanda , tampaknya ini adalah root yang sama (yaitu yang dengan ) daripada pada kasus pertama.tanϕ>0

x=1tanϕ(B2tanϕ1B+(B2tanϕ1B)2+2tanϕt.)
tanϕ<0tanϕ+Δ

Berikut adalah beberapa kode R.

phi <- pi/8; B <- 2
f <- function(t) (-(1/B - 0.5*B*tan(phi)) + 
       sqrt( (1/B - 0.5*B*tan(phi))**2 + 2 * tan(phi) * t))/tan(phi)
hist(f(runif(1e6)))

histogram 1

Dan dengan :ϕ<0

phi <- -pi/8
hist(f(runif(1e6)))

masukkan deskripsi gambar di sini

Elvis
sumber
1
Saya membuat kesalahan, karena saya mengatur sudut saya keluar dari batas, saya mengerti. Tetapi penjelasan Anda mengapa saya harus menghindari menggunakan pemecah angka masih berkabut bagi saya. Bisakah Anda mencoba menjelaskannya lebih lanjut? id suka mendapatkannya. F(x)
Robert
@ Robert Saya pikir kode Anda berfungsi dengan baik jika nilai benar. Namun, itu mencegah Anda menangkap masalah potensial (bagaimana jika tidak ada solusi antara 0 dan ? Atau jika kedua solusi itu? Atau jika tidak ada solusi nyata?). Pekerjaan tambahan untuk menghindari penggunaan solver yang sudah jadi sepadan. BϕB
Elvis