In [1]:
import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt
In [3]:
p = 10.0
D = 1000.0
sigma = 1.0
def gammapdf( x, k, h ):
return x ** ( k - 1 ) / ( gamma( k ) * h ** k ) * np.exp( - x / h )
e = 2 * np.random.randint( 0, 2, D ) - 1
t1 = np.random.normal( 0, sigma, ( D, p ) )
t2 = np.random.normal( 0, sigma, ( D, p ) )
r = e * np.linalg.norm( t1 - t2, axis = 1 ) ** 2 / ( D * p )
plt.hist( r, bins = 50, normed = True )
# plt.plot( )
x = np.linspace( plt.xlim( )[ 0 ], plt.xlim( )[ 1 ], 1000 )
plt.plot( x, 0.5 * gammapdf( np.abs( x ), p / 2.0, 4.0 * sigma ** 2 / ( D * p ) ) )
plt.show( )
In [28]:
p = 5.0
D = 1000.0
sigma = 1.0
def f( t ):
h = 4.0 * sigma ** 2 / ( p * D )
n = 0.5 * p
return 0.5 * ( ( 1 - h * t ) ** ( -n ) + ( 1 + h * t ) ** (-n) )
def b( t ):
h = 4.0 * sigma ** 2 / ( p * D )
n = 0.5 * p
return np.exp( 0.5 * p * ( p + 1 ) * ( 2 ** ( n + 2 ) + ( 2.0 / 3.0 ) ** ( n + 2 ) ) * h ** 2 * t ** 2.0 / 2.0 )
x = np.linspace( - ( p * D ) / ( 8 * sigma ** 2 ), ( p * D ) / ( 8 * sigma ** 2 ), 1000 )
print f( x ) - b( x )
plt.plot( x, f( x ) )
plt.plot( x, b( x ))
plt.show( )
[ -3.60485424e+18 -3.03850694e+18 -2.56201406e+18 -2.16098393e+18
-1.82335116e+18 -1.53899728e+18 -1.29943374e+18 -1.09753701e+18
-9.27327132e+17 -7.83782405e+17 -6.62684489e+17 -5.60488690e+17
-4.74215401e+17 -4.01359186e+17 -3.39812627e+17 -2.87802506e+17
-2.43836312e+17 -2.06657392e+17 -1.75207325e+17 -1.48594360e+17
-1.26066921e+17 -1.06991363e+17 -9.08332901e+16 -7.71418629e+16
-6.55366070e+16 -5.56963244e+16 -4.73497696e+16 -4.02678051e+16
-3.42568004e+16 -2.91530763e+16 -2.48182260e+16 -2.11351742e+16
-1.80048572e+16 -1.53434239e+16 -1.30798764e+16 -1.11540798e+16
-9.51508353e+15 -8.11970437e+15 -6.93133009e+15 -5.91890949e+15
-5.05609899e+15 -4.32054159e+15 -3.69325729e+15 -3.15812773e+15
-2.70146018e+15 -2.31161872e+15 -1.97871199e+15 -1.69432892e+15
-1.45131476e+15 -1.24358150e+15 -1.06594709e+15 -9.13999154e+14
-7.83979475e+14 -6.72685930e+14 -5.77389322e+14 -4.95762763e+14
-4.25821723e+14 -3.65873095e+14 -3.14471931e+14 -2.70384663e+14
-2.32557833e+14 -2.00091501e+14 -1.72216622e+14 -1.48275789e+14
-1.27706844e+14 -1.10028920e+14 -9.48305534e+13 -8.17595459e+13
-7.05143341e+13 -6.08366232e+13 -5.25051064e+13 -4.53301064e+13
-3.91490012e+13 -3.38223194e+13 -2.92304055e+13 -2.52705704e+13
-2.18546574e+13 -1.89069604e+13 -1.63624441e+13 -1.41652222e+13
-1.22672540e+13 -1.06272300e+13 -9.20961663e+12 -7.98383912e+12
-6.92358064e+12 -6.00618195e+12 -5.21212673e+12 -4.52460011e+12
-3.92911002e+12 -3.41316222e+12 -2.96598139e+12 -2.57827161e+12
-2.24201061e+12 -1.95027302e+12 -1.69707844e+12 -1.47726079e+12
-1.28635600e+12 -1.12050536e+12 -9.76372297e+11 -8.51070844e+11
-7.42103927e+11 -6.47310272e+11 -5.64818640e+11 -4.93008373e+11
-4.30475384e+11 -3.76002819e+11 -3.28535748e+11 -2.87159320e+11
-2.51079903e+11 -2.19608802e+11 -1.92148184e+11 -1.68178927e+11
-1.47250109e+11 -1.28969918e+11 -1.12997796e+11 -9.90376388e+10
-8.68319059e+10 -7.61565291e+10 -6.68164971e+10 -5.86420329e+10
-5.14852805e+10 -4.52174328e+10 -3.97262410e+10 -3.49138544e+10
-3.06949459e+10 -2.69950850e+10 -2.37493256e+10 -2.09009783e+10
-1.84005456e+10 -1.62047949e+10 -1.42759536e+10 -1.25810090e+10
-1.10910986e+10 -9.78098099e+09 -8.62857375e+09 -7.61455220e+09
-6.72199946e+09 -5.93610168e+09 -5.24388238e+09 -4.63397092e+09
-4.09640058e+09 -3.62243238e+09 -3.20440136e+09 -2.83558245e+09
-2.51007330e+09 -2.22269198e+09 -1.96888755e+09 -1.74466197e+09
-1.54650180e+09 -1.37131844e+09 -1.21639592e+09 -1.07934516e+09
-9.58063936e+08 -8.50701874e+08 -7.55629734e+08 -6.71412529e+08
-5.96785957e+08 -5.30635746e+08 -4.71979528e+08 -4.19950955e+08
-3.73785747e+08 -3.32809453e+08 -2.96426703e+08 -2.64111763e+08
-2.35400249e+08 -2.09881831e+08 -1.87193827e+08 -1.67015568e+08
-1.49063439e+08 -1.33086521e+08 -1.18862747e+08 -1.06195522e+08
-9.49107479e+07 -8.48541996e+07 -7.58892112e+07 -6.78946399e+07
-6.07630677e+07 -5.43992167e+07 -4.87185500e+07 -4.36460367e+07
-3.91150622e+07 -3.50664651e+07 -3.14476868e+07 -2.82120194e+07
-2.53179416e+07 -2.27285310e+07 -2.04109440e+07 -1.83359565e+07
-1.64775561e+07 -1.48125824e+07 -1.33204076e+07 -1.19826534e+07
-1.07829411e+07 -9.70666908e+06 -8.74081573e+06 -7.87376515e+06
-7.09515195e+06 -6.39572374e+06 -5.76721900e+06 -5.20225866e+06
-4.69424975e+06 -4.23729980e+06 -3.82614077e+06 -3.45606138e+06
-3.12284695e+06 -2.82272587e+06 -2.55232196e+06 -2.30861207e+06
-2.08888830e+06 -1.89072436e+06 -1.71194566e+06 -1.55060249e+06
-1.40494631e+06 -1.27340845e+06 -1.15458119e+06 -1.04720082e+06
-9.50132592e+05 -8.62357225e+05 -7.82958859e+05 -7.11114321e+05
-6.46083504e+05 -5.87200781e+05 -5.33867324e+05 -4.85544232e+05
-4.41746381e+05 -4.02036925e+05 -3.66022367e+05 -3.33348149e+05
-3.03694697e+05 -2.76773882e+05 -2.52325846e+05 -2.30116153e+05
-2.09933238e+05 -1.91586117e+05 -1.74902330e+05 -1.59726097e+05
-1.45916664e+05 -1.33346808e+05 -1.21901508e+05 -1.11476738e+05
-1.01978387e+05 -9.33212903e+04 -8.54283544e+04 -7.82297696e+04
-7.16623032e+04 -6.56686623e+04 -6.01969197e+04 -5.51999972e+04
-5.06351994e+04 -4.64637941e+04 -4.26506335e+04 -3.91638128e+04
-3.59743623e+04 -3.30559689e+04 -3.03847258e+04 -2.79389052e+04
-2.56987543e+04 -2.36463098e+04 -2.17652312e+04 -2.00406492e+04
-1.84590297e+04 -1.70080500e+04 -1.56764865e+04 -1.44541142e+04
-1.33316149e+04 -1.23004938e+04 -1.13530049e+04 -1.04820826e+04
-9.68128003e+03 -8.94471322e+03 -8.26701014e+03 -7.64326480e+03
-7.06899546e+03 -6.54010669e+03 -6.05285497e+03 -5.60381744e+03
-5.18986344e+03 -4.80812877e+03 -4.45599219e+03 -4.13105411e+03
-3.83111720e+03 -3.55416872e+03 -3.29836447e+03 -3.06201420e+03
-2.84356826e+03 -2.64160551e+03 -2.45482226e+03 -2.28202224e+03
-2.12210737e+03 -1.97406947e+03 -1.83698255e+03 -1.70999592e+03
-1.59232783e+03 -1.48325964e+03 -1.38213055e+03 -1.28833276e+03
-1.20130707e+03 -1.12053884e+03 -1.04555427e+03 -9.75917076e+02
-9.11225374e+02 -8.51108867e+02 -7.95226265e+02 -7.43262917e+02
-6.94928651e+02 -6.49955792e+02 -6.08097341e+02 -5.69125321e+02
-5.32829241e+02 -4.99014709e+02 -4.67502139e+02 -4.38125582e+02
-4.10731645e+02 -3.85178500e+02 -3.61334973e+02 -3.39079711e+02
-3.18300412e+02 -2.98893125e+02 -2.80761595e+02 -2.63816674e+02
-2.47975769e+02 -2.33162339e+02 -2.19305433e+02 -2.06339259e+02
-1.94202798e+02 -1.82839436e+02 -1.72196633e+02 -1.62225618e+02
-1.52881104e+02 -1.44121026e+02 -1.35906305e+02 -1.28200622e+02
-1.20970214e+02 -1.14183689e+02 -1.07811845e+02 -1.01827517e+02
-9.62054224e+01 -9.09220256e+01 -8.59554122e+01 -8.12851708e+01
-7.68922849e+01 -7.27590326e+01 -6.88688938e+01 -6.52064648e+01
-6.17573788e+01 -5.85082324e+01 -5.54465174e+01 -5.25605582e+01
-4.98394531e+01 -4.72730204e+01 -4.48517483e+01 -4.25667483e+01
-4.04097121e+01 -3.83728718e+01 -3.64489627e+01 -3.46311887e+01
-3.29131909e+01 -3.12890171e+01 -2.97530952e+01 -2.83002066e+01
-2.69254632e+01 -2.56242851e+01 -2.43923797e+01 -2.32257230e+01
-2.21205416e+01 -2.10732962e+01 -2.00806662e+01 -1.91395356e+01
-1.82469791e+01 -1.74002502e+01 -1.65967694e+01 -1.58341135e+01
-1.51100054e+01 -1.44223050e+01 -1.37690001e+01 -1.31481986e+01
-1.25581204e+01 -1.19970911e+01 -1.14635346e+01 -1.09559673e+01
-1.04729922e+01 -1.00132938e+01 -9.57563240e+00 -9.15884019e+00
-8.76181630e+00 -8.38352287e+00 -8.02298120e+00 -7.67926813e+00
-7.35151268e+00 -7.03889289e+00 -6.74063289e+00 -6.45600009e+00
-6.18430266e+00 -5.92488706e+00 -5.67713578e+00 -5.44046522e+00
-5.21432366e+00 -4.99818948e+00 -4.79156928e+00 -4.59399635e+00
-4.40502905e+00 -4.22424941e+00 -4.05126172e+00 -3.88569131e+00
-3.72718332e+00 -3.57540157e+00 -3.43002751e+00 -3.29075924e+00
-3.15731056e+00 -3.02941007e+00 -2.90680040e+00 -2.78923738e+00
-2.67648934e+00 -2.56833641e+00 -2.46456988e+00 -2.36499157e+00
-2.26941328e+00 -2.17765626e+00 -2.08955064e+00 -2.00493504e+00
-1.92365605e+00 -1.84556780e+00 -1.77053163e+00 -1.69841564e+00
-1.62909434e+00 -1.56244837e+00 -1.49836410e+00 -1.43673340e+00
-1.37745333e+00 -1.32042585e+00 -1.26555759e+00 -1.21275965e+00
-1.16194727e+00 -1.11303974e+00 -1.06596010e+00 -1.02063501e+00
-9.76994541e-01 -9.34972025e-01 -8.94503872e-01 -8.55529433e-01
-8.17990851e-01 -7.81832926e-01 -7.47002985e-01 -7.13450761e-01
-6.81128277e-01 -6.49989738e-01 -6.19991423e-01 -5.91091592e-01
-5.63250391e-01 -5.36429763e-01 -5.10593366e-01 -4.85706492e-01
-4.61735995e-01 -4.38650217e-01 -4.16418924e-01 -3.95013236e-01
-3.74405576e-01 -3.54569603e-01 -3.35480164e-01 -3.17113240e-01
-2.99445896e-01 -2.82456238e-01 -2.66123366e-01 -2.50427334e-01
-2.35349109e-01 -2.20870535e-01 -2.06974297e-01 -1.93643888e-01
-1.80863574e-01 -1.68618370e-01 -1.56894005e-01 -1.45676898e-01
-1.34954134e-01 -1.24713437e-01 -1.14943146e-01 -1.05632199e-01
-9.67701072e-02 -8.83469376e-02 -8.03532956e-02 -7.27803070e-02
-6.56196025e-02 -5.88633023e-02 -5.25040023e-02 -4.65347605e-02
-4.09490850e-02 -3.57409225e-02 -3.09046474e-02 -2.64350521e-02
-2.23273377e-02 -1.85771060e-02 -1.51803513e-02 -1.21334539e-02
-9.43317344e-03 -7.07664371e-03 -5.06136736e-03 -3.38521175e-03
-2.04640524e-03 -1.04353414e-03 -3.75540327e-04 -4.17193742e-05
-4.17193742e-05 -3.75540327e-04 -1.04353414e-03 -2.04640524e-03
-3.38521175e-03 -5.06136736e-03 -7.07664371e-03 -9.43317344e-03
-1.21334539e-02 -1.51803513e-02 -1.85771060e-02 -2.23273377e-02
-2.64350521e-02 -3.09046474e-02 -3.57409225e-02 -4.09490850e-02
-4.65347605e-02 -5.25040023e-02 -5.88633023e-02 -6.56196025e-02
-7.27803070e-02 -8.03532956e-02 -8.83469376e-02 -9.67701072e-02
-1.05632199e-01 -1.14943146e-01 -1.24713437e-01 -1.34954134e-01
-1.45676898e-01 -1.56894005e-01 -1.68618370e-01 -1.80863574e-01
-1.93643888e-01 -2.06974297e-01 -2.20870535e-01 -2.35349109e-01
-2.50427334e-01 -2.66123366e-01 -2.82456238e-01 -2.99445896e-01
-3.17113240e-01 -3.35480164e-01 -3.54569603e-01 -3.74405576e-01
-3.95013236e-01 -4.16418924e-01 -4.38650217e-01 -4.61735995e-01
-4.85706492e-01 -5.10593366e-01 -5.36429763e-01 -5.63250391e-01
-5.91091592e-01 -6.19991423e-01 -6.49989738e-01 -6.81128277e-01
-7.13450761e-01 -7.47002985e-01 -7.81832926e-01 -8.17990851e-01
-8.55529433e-01 -8.94503872e-01 -9.34972025e-01 -9.76994541e-01
-1.02063501e+00 -1.06596010e+00 -1.11303974e+00 -1.16194727e+00
-1.21275965e+00 -1.26555759e+00 -1.32042585e+00 -1.37745333e+00
-1.43673340e+00 -1.49836410e+00 -1.56244837e+00 -1.62909434e+00
-1.69841564e+00 -1.77053163e+00 -1.84556780e+00 -1.92365605e+00
-2.00493504e+00 -2.08955064e+00 -2.17765626e+00 -2.26941328e+00
-2.36499157e+00 -2.46456988e+00 -2.56833641e+00 -2.67648934e+00
-2.78923738e+00 -2.90680040e+00 -3.02941007e+00 -3.15731056e+00
-3.29075924e+00 -3.43002751e+00 -3.57540157e+00 -3.72718332e+00
-3.88569131e+00 -4.05126172e+00 -4.22424941e+00 -4.40502905e+00
-4.59399635e+00 -4.79156928e+00 -4.99818948e+00 -5.21432366e+00
-5.44046522e+00 -5.67713578e+00 -5.92488706e+00 -6.18430266e+00
-6.45600009e+00 -6.74063289e+00 -7.03889289e+00 -7.35151268e+00
-7.67926813e+00 -8.02298120e+00 -8.38352287e+00 -8.76181630e+00
-9.15884019e+00 -9.57563240e+00 -1.00132938e+01 -1.04729922e+01
-1.09559673e+01 -1.14635346e+01 -1.19970911e+01 -1.25581204e+01
-1.31481986e+01 -1.37690001e+01 -1.44223050e+01 -1.51100054e+01
-1.58341135e+01 -1.65967694e+01 -1.74002502e+01 -1.82469791e+01
-1.91395356e+01 -2.00806662e+01 -2.10732962e+01 -2.21205416e+01
-2.32257230e+01 -2.43923797e+01 -2.56242851e+01 -2.69254632e+01
-2.83002066e+01 -2.97530952e+01 -3.12890171e+01 -3.29131909e+01
-3.46311887e+01 -3.64489627e+01 -3.83728718e+01 -4.04097121e+01
-4.25667483e+01 -4.48517483e+01 -4.72730204e+01 -4.98394531e+01
-5.25605582e+01 -5.54465174e+01 -5.85082324e+01 -6.17573788e+01
-6.52064648e+01 -6.88688938e+01 -7.27590326e+01 -7.68922849e+01
-8.12851708e+01 -8.59554122e+01 -9.09220256e+01 -9.62054224e+01
-1.01827517e+02 -1.07811845e+02 -1.14183689e+02 -1.20970214e+02
-1.28200622e+02 -1.35906305e+02 -1.44121026e+02 -1.52881104e+02
-1.62225618e+02 -1.72196633e+02 -1.82839436e+02 -1.94202798e+02
-2.06339259e+02 -2.19305433e+02 -2.33162339e+02 -2.47975769e+02
-2.63816674e+02 -2.80761595e+02 -2.98893125e+02 -3.18300412e+02
-3.39079711e+02 -3.61334973e+02 -3.85178500e+02 -4.10731645e+02
-4.38125582e+02 -4.67502139e+02 -4.99014709e+02 -5.32829241e+02
-5.69125321e+02 -6.08097341e+02 -6.49955792e+02 -6.94928651e+02
-7.43262917e+02 -7.95226265e+02 -8.51108867e+02 -9.11225374e+02
-9.75917076e+02 -1.04555427e+03 -1.12053884e+03 -1.20130707e+03
-1.28833276e+03 -1.38213055e+03 -1.48325964e+03 -1.59232783e+03
-1.70999592e+03 -1.83698255e+03 -1.97406947e+03 -2.12210737e+03
-2.28202224e+03 -2.45482226e+03 -2.64160551e+03 -2.84356826e+03
-3.06201420e+03 -3.29836447e+03 -3.55416872e+03 -3.83111720e+03
-4.13105411e+03 -4.45599219e+03 -4.80812877e+03 -5.18986344e+03
-5.60381744e+03 -6.05285497e+03 -6.54010669e+03 -7.06899546e+03
-7.64326480e+03 -8.26701014e+03 -8.94471322e+03 -9.68128003e+03
-1.04820826e+04 -1.13530049e+04 -1.23004938e+04 -1.33316149e+04
-1.44541142e+04 -1.56764865e+04 -1.70080500e+04 -1.84590297e+04
-2.00406492e+04 -2.17652312e+04 -2.36463098e+04 -2.56987543e+04
-2.79389052e+04 -3.03847258e+04 -3.30559689e+04 -3.59743623e+04
-3.91638128e+04 -4.26506335e+04 -4.64637941e+04 -5.06351994e+04
-5.51999972e+04 -6.01969197e+04 -6.56686623e+04 -7.16623032e+04
-7.82297696e+04 -8.54283544e+04 -9.33212903e+04 -1.01978387e+05
-1.11476738e+05 -1.21901508e+05 -1.33346808e+05 -1.45916664e+05
-1.59726097e+05 -1.74902330e+05 -1.91586117e+05 -2.09933238e+05
-2.30116153e+05 -2.52325846e+05 -2.76773882e+05 -3.03694697e+05
-3.33348149e+05 -3.66022367e+05 -4.02036925e+05 -4.41746381e+05
-4.85544232e+05 -5.33867324e+05 -5.87200781e+05 -6.46083504e+05
-7.11114321e+05 -7.82958859e+05 -8.62357225e+05 -9.50132592e+05
-1.04720082e+06 -1.15458119e+06 -1.27340845e+06 -1.40494631e+06
-1.55060249e+06 -1.71194566e+06 -1.89072436e+06 -2.08888830e+06
-2.30861207e+06 -2.55232196e+06 -2.82272587e+06 -3.12284695e+06
-3.45606138e+06 -3.82614077e+06 -4.23729980e+06 -4.69424975e+06
-5.20225866e+06 -5.76721900e+06 -6.39572374e+06 -7.09515195e+06
-7.87376515e+06 -8.74081573e+06 -9.70666908e+06 -1.07829411e+07
-1.19826534e+07 -1.33204076e+07 -1.48125824e+07 -1.64775561e+07
-1.83359565e+07 -2.04109440e+07 -2.27285310e+07 -2.53179416e+07
-2.82120194e+07 -3.14476868e+07 -3.50664651e+07 -3.91150622e+07
-4.36460367e+07 -4.87185500e+07 -5.43992167e+07 -6.07630677e+07
-6.78946399e+07 -7.58892112e+07 -8.48541996e+07 -9.49107479e+07
-1.06195522e+08 -1.18862747e+08 -1.33086521e+08 -1.49063439e+08
-1.67015568e+08 -1.87193827e+08 -2.09881831e+08 -2.35400249e+08
-2.64111763e+08 -2.96426703e+08 -3.32809453e+08 -3.73785747e+08
-4.19950955e+08 -4.71979528e+08 -5.30635746e+08 -5.96785957e+08
-6.71412529e+08 -7.55629734e+08 -8.50701874e+08 -9.58063936e+08
-1.07934516e+09 -1.21639592e+09 -1.37131844e+09 -1.54650180e+09
-1.74466197e+09 -1.96888755e+09 -2.22269198e+09 -2.51007330e+09
-2.83558245e+09 -3.20440136e+09 -3.62243238e+09 -4.09640058e+09
-4.63397092e+09 -5.24388238e+09 -5.93610168e+09 -6.72199946e+09
-7.61455220e+09 -8.62857375e+09 -9.78098099e+09 -1.10910986e+10
-1.25810090e+10 -1.42759536e+10 -1.62047949e+10 -1.84005456e+10
-2.09009783e+10 -2.37493256e+10 -2.69950850e+10 -3.06949459e+10
-3.49138544e+10 -3.97262410e+10 -4.52174328e+10 -5.14852805e+10
-5.86420329e+10 -6.68164971e+10 -7.61565291e+10 -8.68319059e+10
-9.90376388e+10 -1.12997796e+11 -1.28969918e+11 -1.47250109e+11
-1.68178927e+11 -1.92148184e+11 -2.19608802e+11 -2.51079903e+11
-2.87159320e+11 -3.28535748e+11 -3.76002819e+11 -4.30475384e+11
-4.93008373e+11 -5.64818640e+11 -6.47310272e+11 -7.42103927e+11
-8.51070844e+11 -9.76372297e+11 -1.12050536e+12 -1.28635600e+12
-1.47726079e+12 -1.69707844e+12 -1.95027302e+12 -2.24201061e+12
-2.57827161e+12 -2.96598139e+12 -3.41316222e+12 -3.92911002e+12
-4.52460011e+12 -5.21212673e+12 -6.00618195e+12 -6.92358064e+12
-7.98383912e+12 -9.20961663e+12 -1.06272300e+13 -1.22672540e+13
-1.41652222e+13 -1.63624441e+13 -1.89069604e+13 -2.18546574e+13
-2.52705704e+13 -2.92304055e+13 -3.38223194e+13 -3.91490012e+13
-4.53301064e+13 -5.25051064e+13 -6.08366232e+13 -7.05143341e+13
-8.17595459e+13 -9.48305534e+13 -1.10028920e+14 -1.27706844e+14
-1.48275789e+14 -1.72216622e+14 -2.00091501e+14 -2.32557833e+14
-2.70384663e+14 -3.14471931e+14 -3.65873095e+14 -4.25821723e+14
-4.95762763e+14 -5.77389322e+14 -6.72685930e+14 -7.83979475e+14
-9.13999154e+14 -1.06594709e+15 -1.24358150e+15 -1.45131476e+15
-1.69432892e+15 -1.97871199e+15 -2.31161872e+15 -2.70146018e+15
-3.15812773e+15 -3.69325729e+15 -4.32054159e+15 -5.05609899e+15
-5.91890949e+15 -6.93133009e+15 -8.11970437e+15 -9.51508353e+15
-1.11540798e+16 -1.30798764e+16 -1.53434239e+16 -1.80048572e+16
-2.11351742e+16 -2.48182260e+16 -2.91530763e+16 -3.42568004e+16
-4.02678051e+16 -4.73497696e+16 -5.56963244e+16 -6.55366070e+16
-7.71418629e+16 -9.08332901e+16 -1.06991363e+17 -1.26066921e+17
-1.48594360e+17 -1.75207325e+17 -2.06657392e+17 -2.43836312e+17
-2.87802506e+17 -3.39812627e+17 -4.01359186e+17 -4.74215401e+17
-5.60488690e+17 -6.62684489e+17 -7.83782405e+17 -9.27327132e+17
-1.09753701e+18 -1.29943374e+18 -1.53899728e+18 -1.82335116e+18
-2.16098393e+18 -2.56201406e+18 -3.03850694e+18 -3.60485424e+18]
In [156]:
def bound1( D, eps, sigma, p ):
q = ( p + 2 ) * 2 ** ( ( p + 6 ) / 2 )
dsg = ( ( 4 * sigma ** 4 * q ) / ( p * D ) )
delta = ( p * D ) / ( 8 * sigma ** 2 )
print 1 - eps / ( delta * dsg )
return np.exp( - ( eps * delta ) / 2 * ( 1 - ( 1 - eps / ( delta * dsg ) ).clip( 0, np.inf ) ) * ( 1 + ( 1 - ( delta * dsg ) / eps ).clip( 0, np.inf ) ) )
def bound2( D, eps, sigma, p ):
return np.exp( -eps ** 2 * D / ( 2 * 2 * sigma ** 2 ) )
In [157]:
p = 1.0
sigma = 1
D = np.logspace( 0, 13, 1000, 10 )
plt.subplot( 1, 1, 1 )
plt.plot( D, bound1( D, 0.01, sigma, p ), label = "theoretical" )
plt.plot( D, bound2( D, 0.01, sigma, p ), label = "kontorovich" )
plt.xlim( [ 1, 1e13 ] )
plt.xscale( 'log' )
plt.yscale( 'log' )
plt.legend( )
plt.show( )
# eps = np.logspace( 0, -10, 100, 10 )
# plt.subplot( 1, 1, 1 )
# # for i in np.logspace( -5, 0, 6, 10 ):
# # plt.plot( D, bound2( D, 0.1, 1, 3.0 ), label = i )
# plt.plot( D, bound1( 1000, eps, 1, 3.0 ), label = "mine" )
# plt.xlim( [ 1, 1e13 ] )
# plt.xscale( 'log' )
# # plt.yscale( 'log' )
# plt.legend( )
# plt.show( )
[ 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074 0.99941074 0.99941074
0.99941074 0.99941074 0.99941074 0.99941074]
In [1]:
import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt
In [38]:
D = 100
p = 5
eps = 0.1
sigma = 2
def bound( D, eps, sigma, p ):
if eps <= 6 * np.sqrt( 2 ):
return np.exp( -eps ** 2 * D * p / ( 96 * np.sqrt( 2 ) ) )
else:
return np.exp( -eps * D * p / 16.0 )
n = 100
k = 10
prob = np.zeros( k )
D = np.logspace( 0, 4, k )
Dt = np.logspace( 0, 7, 100 )
for idx, nb_draw in enumerate( D ):
for i in xrange( 0, n ):
t = sigma * np.random.randn( p, nb_draw )
if np.linalg.norm( ( np.dot( t, t.T ) / nb_draw - sigma ** 2 * np.eye( p ) ), 2 ) >= sigma ** 2 * p * eps:
prob[ idx ] = prob[ idx ] + 1
plt.subplot( 1, 1, 1 )
plt.plot( D, prob / n, marker = '+', color = 'b' )
plt.plot( Dt, bound( Dt, eps, sigma, p ), marker = '+', color = 'k' )
plt.axhline(y=1-1.0/n, color = 'r')
plt.axhline(y=1.0/n, color = 'r')
print prob / n
print bound( Dt, eps, sigma, p )
plt.xscale( 'log' )
plt.show( )
[ 1. 1. 1. 1. 0.37 0.02 0. 0. 0. 0. ]
[ 9.99631783e-001 9.99566692e-001 9.99490097e-001 9.99399968e-001
9.99293912e-001 9.99169119e-001 9.99022282e-001 9.98849509e-001
9.98646226e-001 9.98407053e-001 9.98125665e-001 9.97794626e-001
9.97405196e-001 9.96947105e-001 9.96408287e-001 9.95774572e-001
9.95029326e-001 9.94153025e-001 9.93122773e-001 9.91911727e-001
9.90488444e-001 9.88816125e-001 9.86851735e-001 9.84545015e-001
9.81837345e-001 9.78660463e-001 9.74935035e-001 9.70569065e-001
9.65456177e-001 9.59473757e-001 9.52481037e-001 9.44317169e-001
9.34799412e-001 9.23721600e-001 9.10853127e-001 8.95938786e-001
8.78699902e-001 8.58837345e-001 8.36037161e-001 8.09979701e-001
7.80353294e-001 7.46873545e-001 7.09309285e-001 6.67515884e-001
6.21475937e-001 5.71346101e-001 5.17506910e-001 4.60609732e-001
4.01611721e-001 3.41786269e-001 2.82694171e-001 2.26101329e-001
1.73834586e-001 1.27580172e-001 8.86492292e-002 5.77578900e-002
3.48858572e-002 1.92740239e-002 9.58815430e-003 4.21581594e-003
1.60299559e-003 5.13725301e-004 1.34632179e-004 2.78442887e-005
4.35823887e-006 4.91447518e-007 3.76751425e-008 1.83405635e-009
5.23214310e-011 7.95828123e-013 5.77483416e-015 1.75391894e-017
1.91139360e-020 6.23529699e-024 4.91938179e-028 7.30321207e-033
1.51844499e-038 3.12333374e-045 4.22216281e-053 2.31794730e-062
2.93303456e-073 4.39209692e-086 3.55202556e-101 6.16310599e-119
7.74146782e-140 1.96021883e-164 2.22242594e-193 1.92062991e-227
1.57487728e-267 1.05519080e-314 0.00000000e+000 0.00000000e+000
0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000
0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000]
In [9]:
In [ ]:
Content source: RomainBrault/OVFM
Similar notebooks: