generated from klb2/reproducible-paper-python-template
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rayleigh.py
98 lines (79 loc) · 3.07 KB
/
rayleigh.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
import numpy as np
from scipy import special
def pdf_skg_rate_rayleigh(s, lam_x, lam_y):
_denom = (lam_y + lam_x * (2**s - 1)) ** 2
_factor = np.log(2) * lam_x * lam_y
_part1 = 2**s * np.exp(lam_x * (1 - 2**s))
_part2 = 1 + lam_y + lam_x * (2**s - 1)
_pdf = _factor * _part1 * _part2 / _denom
_pdf = np.where(s < 0, 0.0, _pdf)
_pdf = np.maximum(_pdf, 0)
return _pdf
def logpdf_skg_rate_rayleigh(s, lam_x, lam_y):
_denom = 2 * np.log(lam_y + lam_x * (2**s - 1))
_factor = np.log(np.log(2)) + np.log(lam_x) + np.log(lam_y)
_part1 = s * np.log(2) + (lam_x * (1 - 2**s))
_part2 = np.log(1 + lam_y + lam_x * (2**s - 1))
_pdf = _factor + _part1 + _part2 - _denom
_pdf = np.where(s < 0, -np.infty, _pdf)
# _pdf = np.maximum(_pdf, 0)
return _pdf
def cdf_skg_rate_rayleigh(s, lam_x, lam_y):
_denom = lam_y + lam_x * (2**s - 1)
_enum = lam_y * np.exp(lam_x * (1 - 2**s))
_cdf = 1 - _enum / _denom
_cdf = np.where(s <= 0, 0.0, _cdf)
_cdf = np.clip(_cdf, 0, 1)
return _cdf
def sf_skg_rate_rayleigh(s, lam_x, lam_y):
return 1.0 - cdf_skg_rate_rayleigh(s, lam_x, lam_y)
def expectation_skg_rate_rayleigh(power, lam_bob, lam_eve):
if power <= 0:
return 0
if lam_bob == lam_eve:
expect = (power + np.exp(lam_eve / power) * special.expi(-lam_eve / power)) / (
power * np.log(2)
)
else:
_part1 = np.exp(lam_bob / power) * special.expi(-lam_bob / power)
_part2 = np.exp(lam_eve / power) * special.expi(-lam_eve / power)
_denom = (lam_bob - lam_eve) * np.log(2)
expect = lam_eve * (_part1 - _part2) / _denom
return expect
def expectation_skg_rate_rayleigh_conditional_bob(power, channel_bob, lam_eve):
if power <= 0:
return 0
_part1 = np.exp(lam_eve / power) * special.expi(-lam_eve / power)
_part2 = -np.exp(lam_eve * (channel_bob + 1 / power)) * special.expi(
-(lam_eve + channel_bob * lam_eve * power) / power
)
_part3 = np.log(1 + channel_bob * power)
conditional_expect = (_part1 + _part2 + _part3) / np.log(2)
return conditional_expect
if __name__ == "__main__":
from scipy import stats
import matplotlib.pyplot as plt
num_samples = 40000
power_db = 2
lam_x = 10 ** (-(10 + power_db) / 10.0)
lam_y = 10 ** (-power_db / 10.0)
prob_tx = 0.3
cond_expect = expectation_skg_rate_rayleigh_conditional_bob(10, 5, 1)
assert np.round(cond_expect, 2) == 3.01
x = stats.expon.rvs(scale=1.0 / lam_x, size=num_samples)
y = stats.expon.rvs(scale=1.0 / lam_y, size=num_samples)
skg = np.log2((1 + x + y) / (1 + y))
s = np.linspace(-2, max(skg) * 1.1, 200)
# s = np.linspace(-2, 100, 200)
from ultimate_ruin_prob import calculate_ultimate_ruin_mixed
b, out = calculate_ultimate_ruin_mixed(
lambda s: np.exp(logpdf_skg_rate_rayleigh(s, lam_x, lam_y)),
lambda s: sf_skg_rate_rayleigh(s, lam_x, lam_y),
message_length=5,
prob_tx=prob_tx,
max_b=210,
num_points=9000,
)
plt.figure()
plt.plot(b, out)
plt.show()