function walk(num_steps)
x = 0
for i in 1:num_steps
x += rand( (-1, +1) )
end
return x
end
walk (generic function with 1 method)
num_steps = 20
walk(num_steps)
0
experiment(num_steps, num_walks) = [walk(num_steps) for i in 1:num_walks]
experiment (generic function with 1 method)
data = experiment(20, 1000)
data' # transpose(data) -- ' = transpose
1×1000 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}: 2 -2 -10 -4 12 2 0 2 2 -2 4 … 4 -4 6 2 -2 2 4 0 0 2 -4
using StatsBase
counts = countmap(data);
using Plots
histogram(data, bins=10) # in Plots.jl package
xlabel!("position x")
using StatsBase
h = fit(Histogram, data)
Histogram{Int64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}} edges: -20.0:5.0:15.0 weights: [1, 8, 109, 296, 438, 121, 27] closed: left isdensity: false
h.weights / sum(h.weights)
7-element Array{Float64,1}: 0.001 0.008 0.109 0.296 0.438 0.121 0.027
histogram(data, bins=30, normed=true) # in Plots.jl package
data
1000-element Array{Int64,1}: 2 -2 -10 -4 12 2 0 2 2 -2 4 0 -6 ⋮ 2 4 -4 6 2 -2 2 4 0 0 2 -4
mean(data)
0.116
sum(data) / length(data)
0.116
histogram(data, normed=true, bins=30)
vline!([mean(data)], lw=3, ls=:dash, leg=false, c=:purple)
scatter(data)
hline!([mean(data)], lw=3, ls=:dash, c=:purple)
m = mean(data) # de-mean: subtract:
centred_data = data .- m;
scatter(centred_data)
hline!([mean(centred_data)], lw=3, ls=:dash, c=:purple)
mean(centred_data)
2.9665159217984184e-16
1e-16
is notation for $10^{-16}$
nextfloat(1.0)
1.0000000000000002
nextfloat(1.0) - 1.0 # machine epsilon for double precision -- Float64
2.220446049250313e-16
abs_data = abs.(centred_data)
1000-element Array{Float64,1}: 1.884 2.116 10.116 4.116 11.884 1.884 0.116 1.884 1.884 2.116 3.884 0.116 6.116 ⋮ 1.884 3.884 4.116 5.884 1.884 2.116 1.884 3.884 0.116 0.116 1.884 4.116
scatter(abs_data)
mean(abs_data)
3.595024
scatter(abs_data, alpha=0.2)
hline!([mean(abs_data)], lw=3, ls=:dash, c=:purple)
m = mean(abs_data)
count(abs_data .< m) / length(abs_data)
0.491
count(abs_data .< 2m) / length(abs_data) # 2m instead of 2*m
0.889
scatter(abs_data, alpha=0.2)
hline!([2*mean(abs_data)], lw=3, ls=:dash, c=:purple)
squared_data = data.^2;
scatter(squared_data, alpha=0.2)
@time data = experiment(50, 10^6);
0.356394 seconds (2 allocations: 7.629 MiB)
using StatsBase
m = mean(data)
variance = mean((data .- m).^2)
σ = √(variance) # \sigma<TAB> # \sqrt<TAB> # standard deviation
7.066092713168997
[(m - σ) < x < (m + σ) for x in data];
count((m - σ) .< data .< (m + σ)) / length(data)
0.677962
data2 = data[1:5_000];
scatter(data2, alpha=0.2)
hline!([m, m+σ, m-σ], lw=3, ls=:dash, leg=false)
hline!([m, m+2σ, m-2σ], lw=3, ls=:dash, leg=false)
count((m - 2σ) .< data .< (m + 2σ)) / length(data)
0.967248
$N$ data points sampled from the same distribution, calculate the mean
num_steps = 20
N = 100
data = experiment(num_steps, N)
m = mean(data)
-0.22
num_steps = 20
N = 1000
data = experiment(num_steps, N)
m = mean(data)
0.164
using Statistics
σ = std(data)
4.190116116896174
$X_1 + X_2 + \cdots + X_N$
standard_error_of_mean = σ / √N
0.1325031058997225
mean_data = [mean(experiment(num_steps, N)) for i in 1:1000];
std(mean_data)
0.14236756134219103
mean_data = [mean(experiment(num_steps, 4N)) for i in 1:1000];
std(mean_data)
0.07322830162689259
histogram(mean_data)
Random walker:
x = 0
x += 1
1
x2 = 0
x2 -= 1
-1
y = 0.0
y += randn() # random number from a gaussian distribution
-0.41883079557282155
histogram(randn(1000)) # ziggurat algorithm for normally distributed random numbers
# See also Distributions.jl
How write generic code? A single walk
function to generate random walk that works with discrete and continuous distributions
discrete_jump() = rand( (-1, +1 ))
continuous_jump() = randn()
continuous_jump (generic function with 1 method)
function walk(N)
x = 0
for i in 1:N
x += jump()
end
return x
end
walk (generic function with 1 method)
# DON'T DO THIS!
function walk(N, walker_type)
if walker_type == "discrete"
...
end
function walk(N, jump)
x = 0 # integer
for i in 1:N
x += jump()
end
return x
end
walk (generic function with 2 methods)
walk(20, discrete_jump)
-4
walk(20, continuous_jump)
6.0904266229323145
using BenchmarkTools
@btime walk(20, discrete_jump)
125.640 ns (0 allocations: 0 bytes)
-4
@btime walk(20, continuous_jump)
155.878 ns (0 allocations: 0 bytes)
1.0104291115561674
function walk(N, jump, x0)
x = x0
for i in 1:N
x += jump()
end
return x
end
walk (generic function with 3 methods)
walk(20, continuous_jump, 0.0)
-4.3678555053328
Julia has objects (of user-defined types), but not object-oriented prog lang
struct MyDiscreteRandomWalker
x::Int # type annotation: x is of type Integer
end
"Box with data inside"
To make an object of this type:
MyDiscreteRandomWalker
MyDiscreteRandomWalker
w = MyDiscreteRandomWalker(10)
MyDiscreteRandomWalker(10)
w
MyDiscreteRandomWalker(10)
w.x # the variable x that lives inside the **object** w
10
Move the walker:
w.x = w.x + 1
setfield! immutable struct of type MyDiscreteRandomWalker cannot be changed Stacktrace: [1] setproperty!(::MyDiscreteRandomWalker, ::Symbol, ::Int64) at ./Base.jl:34 [2] top-level scope at In[91]:1
mutable struct MyDiscreteRandomWalker2
x::Int # type annotation: x is of type Integer
end
w = MyDiscreteRandomWalker2(10)
MyDiscreteRandomWalker2(10)
w.x += 1
11
jump!(w::MyDiscreteRandomWalker2) = w.x += rand( (-1, +1) )
jump! (generic function with 1 method)
function jump!(w::MyDiscreteRandomWalker2) # takes arg of type MyDiscreteRandomWalker2
w.x += rand( (-1, +1) )
return w
end
jump! (generic function with 1 method)
w
MyDiscreteRandomWalker2(11)
jump!(w)
MyDiscreteRandomWalker2(10)
w
MyDiscreteRandomWalker2(10)
function walk!(w::MyDiscreteRandomWalker2, N)
for i in 1:N
jump!(w)
end
end
walk! (generic function with 1 method)
w
MyDiscreteRandomWalker2(10)
w = MyDiscreteRandomWalker2(0)
MyDiscreteRandomWalker2(0)
MyDiscreteRandomWalker2() = MyDiscreteRandomWalker2(0)
MyDiscreteRandomWalker2
methods(MyDiscreteRandomWalker2)
w = MyDiscreteRandomWalker2()
MyDiscreteRandomWalker2(0)
w2 = MyDiscreteRandomWalker2()
MyDiscreteRandomWalker2(0)
walk!(w, 1)
w
MyDiscreteRandomWalker2(-1)
w2
MyDiscreteRandomWalker2(0)
w.x
-1
w2.x
0