-
Notifications
You must be signed in to change notification settings - Fork 326
/
exo4.jl
38 lines (32 loc) · 818 Bytes
/
exo4.jl
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
klist = [10, 30, 50]
P = 200
ntrials = 200
tmin = 0
tmax = 2.5
q = 50
t = linspace(tmin, tmax, q)
t1 = linspace(tmin, tmax, 1000)
dt = (tmax - tmin)/q
for j in 1:length(klist)
k = klist[j]
# simulation
v = []
for i in 1:ntrials
v = [v; svd(randn(P, k)./sqrt(P))[2].^2]
end
figure(figsize = (10, 10))
subplot(length(klist), 1, j)
h = hist(v, t)[1]
h = h/sum(h)/dt
bar(t[1:end], h, width = 1/20, color = "darkblue", edgecolor="black")
# theoritical law
beta = k/P
a = (1 - sqrt(beta))^2
b = (1 + sqrt(beta))^2
z = sqrt(max(t1 - a, zeros(length(t1))).*max(b - t1, zeros(length(t1))))./(2*pi.*beta.*t1)
plot(t1, z, "r", linewidth = 3)
xlim(tmin, tmax)
ylim(0, maximum(h)*1.05)
title(@sprintf("P = %d, k = %d", P, k))
show()
end