numpy.random.normal การแจกแจงที่แตกต่างกัน: การเลือกค่าจากการแจกแจง

ฉันมีการกระจายพลังงานตามกฎกำลัง และฉันต้องการเลือกพลังงานแบบสุ่มตามการกระจาย ฉันลองทำสิ่งนี้ด้วยตนเองโดยใช้ตัวเลขสุ่ม แต่มันก็ไม่มีประสิทธิภาพเกินไปสำหรับสิ่งที่ฉันต้องการทำ ฉันสงสัยว่ามีวิธีใดที่เป็นตัวเลข (หรืออื่นๆ) ที่ทำงานเหมือน numpy.random.normal ยกเว้นว่าแทนที่จะใช้การแจกแจงแบบปกติ อาจมีการระบุการแจกแจงแทน ดังนั้นในใจของฉันตัวอย่างอาจมีลักษณะดังนี้ (คล้ายกับ numpy.random.normal):

import numpy as np

# Energies from within which I want values drawn
eMin = 50.
eMax = 2500.

# Amount of energies to be drawn
n = 10000

photons = []

for i in range(n):

    # Method that I just made up which would work like random.normal,
    # i.e. return an energy on the distribution based on its probability,
    # but take a distribution other than a normal distribution
    photons.append(np.random.distro(eMin, eMax, lambda e: e**(-1.)))

print(photons)

การพิมพ์ photons ควรให้รายการความยาว 10,000 รายการที่มีพลังงานในการแจกแจงนี้แก่ฉัน ถ้าผมทำฮิสโตแกรม มันจะมีค่าถังขยะมากกว่ามากที่พลังงานต่ำกว่า

ฉันไม่แน่ใจว่ามีวิธีดังกล่าวหรือไม่ แต่ดูเหมือนว่าควรจะเป็นเช่นนั้น ฉันหวังว่ามันชัดเจนว่าฉันต้องการทำอะไร

แก้ไข:

ฉันเคยเห็น numpy.random.power แล้ว แต่เลขชี้กำลังของฉันคือ -1 ดังนั้นฉันไม่คิดว่ามันจะใช้ได้


python numpy random distribution
person davly    schedule 07.07.2014    source แหล่งที่มา
comment
คุณต้องการ PDF อะไรกันแน่? การจ่ายพลังงานเป็นกรณีพิเศษของเบต้า คุณสามารถใช้ docs.scipy.org/doc/numpy/reference/generated/ ?   -  person wim    schedule 07.07.2014
comment
@wim ฉันเชื่อว่าฉันต้องการฟังก์ชั่นแบบแยกส่วนซึ่ง f(x)=0 อยู่นอกช่วงพลังงานของฉันและ f(x)=x**a (โดยที่ a สามารถเป็นค่าได้ตั้งแต่ -5 ถึง 5) ภายใน ฉันไม่เห็นว่าเบต้าจะทำงานที่นี่อย่างไร   -  person davly    schedule 08.07.2014
comment
@davly อัปเดตคำตอบของฉันด้วยข้อมูลโค้ดในกรณีที่มีประโยชน์   -  person John Greenall    schedule 08.07.2014

คำตอบ (3)


arrow_upward
1
arrow_downward

การสุ่มตัวอย่างจากไฟล์ PDF โดยพลการนั้นค่อนข้างยากจริงๆ มีหนังสือขนาดใหญ่และหนาแน่นเกี่ยวกับวิธีสุ่มตัวอย่างจากกลุ่มการกระจายมาตรฐานอย่างมีประสิทธิภาพและแม่นยำ

ดูเหมือนว่าคุณอาจจะใช้วิธีผกผันแบบกำหนดเองสำหรับตัวอย่างที่คุณให้ไว้

person Robert Kern    schedule 07.07.2014
comment
ฉันจะใช้วิธีผกผันแบบกำหนดเองได้อย่างไร ฉันไม่เห็นเหมือนที่นี่ - person davly; 08.07.2014
comment
สืบทอดฟังก์ชันผกผันของ CDF ใช้ random_sample() เพื่อรับค่าที่กระจายสม่ำเสมอระหว่าง 0 ถึง 1 ส่งค่าเหล่านี้ผ่าน CDF ผกผันเพื่อรับค่าที่เป็นไปตามการกระจายที่คุณต้องการ ในกรณีของคุณ CDF แบบผกผันคือ lambda u: eMin*(eMax/eMin)**u - person Robert Kern; 09.07.2014

arrow_upward
1
arrow_downward

หากคุณต้องการสุ่มตัวอย่างจากการแจกแจงตามอำเภอใจ คุณต้องมีฟังก์ชันผกผันของฟังก์ชันความหนาแน่นสะสม (ไม่ใช่ pdf)

จากนั้นคุณสุ่มตัวอย่างความน่าจะเป็นอย่างสม่ำเสมอจากช่วง [0,1] และป้อนค่านี้เข้าไปในค่าผกผันของ cdf เพื่อให้ได้ค่าที่สอดคล้องกัน

มักจะไม่สามารถรับ cdf จาก pdf ในเชิงวิเคราะห์ได้ อย่างไรก็ตาม หากคุณยินดีที่จะประมาณการกระจายตัว คุณสามารถทำได้โดยคำนวณ f(x) ในช่วงเวลาปกติของโดเมน จากนั้นทำผลรวมบนเวกเตอร์นี้เพื่อให้ได้ค่าประมาณของ cdf และจากค่านี้ก็จะประมาณค่าผกผัน

ข้อมูลโค้ดคร่าวๆ:

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate

def f(x):
   """
   substitute this function with your arbitrary distribution
   must be positive over domain
   """
   return 1/float(x)


#you should vary inputVals to cover the domain of f (for better accurracy you can
#be clever about spacing of values as well). Here i space them logarithmically
#up to 1 then at regular intervals but you could definitely do better
inputVals = np.hstack([1.**np.arange(-1000000,0,100),range(1,10000)])

#everything else should just work
funcVals = np.array([f(x) for x in inputVals])
cdf = np.zeros(len(funcVals))
diff = np.diff(funcVals)
for i in xrange(1,len(funcVals)):
   cdf[i] = cdf[i-1]+funcVals[i-1]*diff[i-1]
cdf /= cdf[-1]

#you could also improve the approximation by choosing appropriate interpolator
inverseCdf = scipy.interpolate.interp1d(cdf,inputVals)

#grab 10k samples from distribution
samples = [inverseCdf(x) for x in np.random.uniform(0,1,size = 100000)]

plt.hist(samples,bins=500)
plt.show()
person John Greenall    schedule 07.07.2014

arrow_upward
-1
arrow_downward

ทำไมคุณไม่ใช้ eval และใส่การแจกแจงเป็นสตริงล่ะ?

>>> cmd = "numpy.random.normal(500)"
>>> eval(cmd)

คุณสามารถจัดการสตริงได้ตามที่คุณต้องการเพื่อตั้งค่าการแจกแจง

person BigBrownBear00    schedule 07.07.2014
comment
ขออภัย ฉันเข้าใจผิดคำถามของคุณ - person BigBrownBear00; 07.07.2014