2016-04-30 17:28 GMT+02:00 Francesca Senatore <
francesca.senatore1...@gmail.com>:
> Ciao ragazzi,
> vi ringrazio per i vostri consigli. Scusate se rispondo solo ora ma causa
> influenza sono costretta a letto.
>
> Avete ragione! Non credo di essermi spiegata molto bene. Cerco di
> spiegarmi meglio dicendovi cosa esattamente devo fare. Io voglio simulare
> una distribuzione di exposures (a e b sono secondi di osservazione da 5.6e3
> a 2.5e6 secondi). La distribuzione che devo simulare รจ fatta in questo
> modo:
>
> l' 88.5% delle exposures cade nell'intervallo [2.5e6,1e5]
> il 9.5% delle exposure cade nell'intervallo [1e5,2.5e4]
> 0.7% delle exposures cade nell'intervallo [2.5e4,1.1e4]
> 1.4% delle exposures cade nell'intervallo [1.1e4,6.3e3]
>
Allora rv_discrete e random.choice non fanno al tuo caso. Quello che
potresti fare e` invece interpolare la tua cdf per poterla invertire, ed
utilizzarla poi per trasformare i campioni estratti da una distribuzione
uniforme. La procedura si chiama inverse transform sampling, e.g. [1].
Di seguito un esempio.
from scipy.interpolate import interp1d
import numpy as np
n_samples = 1e4
px = [2.5e6, 1e5, 2.5e4, 1.1e4, 6.3e3]
py = [0.885, 0.095, 0.07, 0.014, 0]
py = py / np.sum(py)
cdf_x = np.flipud(px)
cdf_y = np.cumsum(np.flipud(py))
inv_cdf = interp1d(cdf_y, cdf_x, bounds_error=False, assume_sorted=True)
r = np.random.uniform(0, 1, n_samples)
ys = inv_cdf(r)
for b, a, nom in zip(px[:-1], px[1:], py):
obs = np.sum((ys < b) & (ys > a)) / len(ys)
print('observed: %.3f - nominal: %.3f' % (obs, nom))
[1]
http://www.nehalemlabs.net/prototype/blog/2013/12/16/how-to-do-inverse-transformation-sampling-in-scipy-and-numpy/
Ciao,
L
_______________________________________________
Python mailing list
Python@lists.python.it
http://lists.python.it/mailman/listinfo/python
