User Tools

Site Tools


random
snippet.juliarepl
julia> pkgchk( [ "julia" => v"1.0.2", "Distributions" => v"0.16.4" ] )

Statistical Distributions and Random Variates

An array of random numbers can be created using the rand or randn functions.

The rand function can also be called as rand(T, n) to generate n random numbers of type T.

Base julia has facilities to generate uniform numbers (rand) and Gaussian (Normal) random numbers (randn). Other distributions rely on the Distributions package.

snippet.juliarepl
julia> using Distributions

Drawing Random Variables

rand( 1 ) draws a vector with 1 element, while rand() draws a scalar. Thus, y= 3.0 + rand( 1 ) results in a vector, not a scalar.

The following is a selection of random numbers drawn from some common distributions. For a more complete list of available distributions, refer to Distributions.

Setting the Seed

snippet.juliarepl
julia> using Random;  Random.seed!(0);

Uniform Floats

snippet.juliarepl
julia> using Random; Random.seed!(0);

julia> rand(3,2)            ## multiple int args are array dimensions
3×2 Array{Float64,2}:
 0.823648  0.177329
 0.910357  0.27888
 0.164566  0.203477

Uniform Integers

snippet.juliarepl
julia> using Random; Random.seed!(0); 

julia> rand(Int, 3)         ## Types as argument to rand mean over the whole range
3-element Array{Int64,1}:
 1339893410598768192
 1575814717733606317
 2119082727832020007

Uniform Ranged Integers (e.g., Bernoullis)

snippet.juliarepl
julia> using Random, Distributions; Random.seed!(0);

julia> rand( 1:10, 3 )		## 3 draws from 1 to 10, equally likely
3-element Array{Int64,1}:
  1
  3
 10

julia> Random.seed!(0);

julia> rand( 0:1, 10)		## 10 draws from 0 to 1, equally likely
10-element Array{Int64,1}:
 0
 0
 1
 1
 0
 0
 0
 1
 1
 0

julia> Random.seed!(0);

julia> rand( Bernoulli( 0.15 ), 10 )    ## 10 draws, each is 1 with prob 0.15
10-element Array{Int64,1}:
 0
 0
 0
 0
 0
 0
 1
 1
 0
 0
  • Note: If speed is of the essence, there are faster methods than rand(). For example, to draw 1,000 uniform numbers between 1 and 4, you can replace [ rand(1:4) for i=1:1000 ] with rand( Random.SamplerRangeFast(1:4), 1000) or with 1 .+ Int.(floor.(rand(1000)*4)).

Binomial Integers

snippet.juliarepl
julia> using Random, Distributions;  Random.seed!(0);

julia> rand( Binomial( 10, 0.2 ), 10 )   ## each draw has 0.2 prob of 1, summed up 10 times
10-element Array{Int64,1}:
 3
 4
 1
 1
 1
 1
 0
 0
 1
 5

Poissons

snippet.juliarepl
julia> using Random, Distributions;  Random.seed!(0);

julia> rand( Poisson( 2 ), 20 )
20-element Array{Int64,1}:
 3
 4
 1
 1
 1
 1
 0
 0
 1
 5
 2
 2
 1
 4
 1
 2
 2
 4
 5
 3

Random Orthogonal Matrix

snippet.juliarepl
julia> using Random, LinearAlgebra;  Random.seed!(0);

julia> QR= qr(randn(3,3));        ## QR is already orthogonal

julia> Q= QR.Q * Diagonal(sign.(diag(QR.R)))        ## orthogonal *and* uniformly distributed
3×3 Array{Float64,2}:
  0.602120.4664330.647991
  0.7345     0.641769   0.2205510.312988   0.6087470.729017

julia> Q'*Q                ## show that we have orthogonal vectors
3×3 Array{Float64,2}:
  1.01.94289e–16  0.01.94289e–16   1.0          1.66533e–16
  0.0           1.66533e–16  1.0

Shuffling an Array

Shuffled values from 1:N can be used as indices to rearrange N vector elements.

snippet.juliarepl
julia> using Random, Distributions;  Random.seed!(0);

julia> Random.seed!(0); randperm(5)
5-element Array{Int64,1}:
 5
 4
 3
 1
 2

julia> Random.seed!(0); shuffle( collect(1:5) )		## a different way of shuffling
5-element Array{Int64,1}:
 5
 4
 2
 3
 1

julia> Random.seed!(0); sample( collect(10:20), 3 )	## Pick 3 random numbers from 10:20 array
3-element Array{Int64,1}:
 10
 12
 19

StatsBase offers many sampling variants.

Drawing "Perfect-Density" Random Variables

Independent random draws do not have empirical densities that perfectly match their theoretical distributions. For example, a million fair coin flips may have 500,021 heads, rather than 500,000 heads. Sometimes, this type of true randomness is desirable. Often, it is not and simply adds extra noise where it is not desired. Users may want to have random draws on each draw, but have exactly 500,000 heads and 500,000 tails. The randomness is in when what value shows up.

The following draws from a “perfect density” distribution, guaranteed to match the theoretical density function, with randomness in when each value is drawn:

snippet.juliarepl
julia> using Random, Distributions;  Random.seed!(0); 

julia> randPerfect(d::UnivariateDistribution, n::Integer)= shuffle!(quantile.(d, (1:n)/(n+1)));

julia> randPerfect( Normal(0,1), 5 )	## Example: draws from a normal distribution
5-element Array{Float64,1}:
  0.9674215661017014
  0.430727299295457330.4307272992954576
  0.00.9674215661017014

Fitting Parameters of Statistical Distributions To Data

snippet.juliarepl
julia> using Random, Distributions;  Random.seed!(0);

julia> fit_mle( Normal, rand( Normal( 10 , 5 ) , 10_000_000 ) )    ## Fitting the correct distribution -> recover params
Normal{Float64}(μ=9.99913328006548, σ=5.000592281736906)

julia> fit_mle( Normal, rand(1_000_000) )         ## Fitting the wrong distribution--a normal fitted to uniforms
Normal{Float64}(μ=0.5001626330884961, σ=0.288518549254692)

Testing for Statistical Distributions

See unistats for tests whether a set of random draws come from a particular distributin.

The Gaussian (Normal) Distribution

Random Numbers

snippet.juliarepl
julia> using Random;  Random.seed!(0); 

julia> randn(3)
3-element Array{Float64,1}:
  0.6791074260357777
  0.82841348290003590.3530074003005963

In distributions,

snippet.juliarepl
julia> using Random, Distributions;  Random.seed!(0);

julia> rand( Normal(0,1), 4 )
4-element Array{Float64,1}:
  0.6791074260357777
  0.82841348290003590.35300740030059630.13485387193052173

See also the example of perfect-density draws above, which also generate Normal variates.

Multivariate Gaussian (Normal) Distribution with Covariance

snippet.juliarepl
julia> using Random, Distributions; Random.seed!(0);

julia> gen3dim= MvNormal( 3, 100.0 )    ## generate 3-dim normal SD=100
ZeroMeanIsoNormal(
dim: 3
μ: [0.0, 0.0, 0.0]
Σ: [10000.0 0.0 0.0; 0.0 10000.0 0.0; 0.0 0.0 10000.0]
)


julia> rand( gen3dim, 5 )
3×5 Array{Float64,2}:
  67.910713.4854    6.49475  157.433    39.7482
  82.8413   58.661710.901768.8907   81.16335.3007   29.733651.42176.280434.6355

julia> Random.seed!(0);

julia> x= 10.0*randn(1000); y= 10.0*randn(1000)+x;  covmat= cov(hcat( x, y )); x= y= nothing

julia> gen2dim= MvNormal( [10.0, 20.0], covmat )	## specific mean and covariance matrix
FullNormal(
dim: 2
μ: [10.0, 20.0]
Σ: [106.834 98.4571; 98.4571 190.721]
)


julia> rand( gen2dim, 5 )
2×5 Array{Float64,2}:
 11.4939  23.7481  17.13091.470978.87016
 30.0597  33.7182  14.2002  13.7076    5.54033

Probability Density Function (PDF)

snippet.juliarepl
julia> using Random, Distributions; Random.seed!(0);

julia> pdf.( Normal(0,1), [2.0, 0, 2.0 ] )  ### use logpdf for tiny values
3-element Array{Float64,1}:
 0.05399096651318806
 0.3989422804014327
 0.05399096651318806

Log Likelihood

snippet.juliarepl
julia> using Random, Distributions; Random.seed!(0);

julia> loglikelihood( Normal(0,1), rand( Normal(0,1), 1_000 ) )1452.9989715615147

Cumulative Density Funtion (CDF)

snippet.juliarepl
julia> using Distributions; Random.seed!(0);

julia> cdf.( Normal(0,1), [2.0, 0, 2.0 ] )  ### use logcdf for tiny values
3-element Array{Float64,1}:
 0.022750131948179205
 0.5
 0.9772498680518208

Inverse CDF

snippet.juliarepl
julia> using Distributions; Random.seed!(0);

julia>  quantile.( Normal(0,1), [0.025,0.05,0.5,0.95,0.975] )
5-element Array{Float64,1}:
 –1.95996398454005921.6448536269514724
  0.0
  1.6448536269514717
  1.9599639845400576

Histograms

snippet.juliarepl
julia> using Random, Distributions, StatsBase; Random.seed!(0);

julia> x= rand( Normal(0,1), 10_000);

julia> fit( Histogram, x , closed= :left, nbins=10  )
Histogram{Int64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
edges:
  –4.0:1.0:5.0
weights: [10, 224, 1398, 3426, 3373, 1342, 213, 13, 1]
closed: left
isdensity: false

julia> fit( Histogram, x , closed= :right, nbins=10  )
Histogram{Int64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
edges:
  –4.0:1.0:5.0
weights: [10, 224, 1398, 3426, 3373, 1342, 213, 13, 1]
closed: right
isdensity: false

Julia's plotting makes it easy to plot histograms:

snippet.juliaeval
[download only julia statements]
using Plots
using Distributions, StatsBase; n=1_000; rn= quantile.( Normal(0.0,1.0), [1:n;]/(n+1));
histogram( rn.^2, legend=false );	# a chi-squared
savefig( "plotting/histogram1.png" );

  • StatsBase contains an ecdf function that is the direct empirical cdf, whereas EmpiricalUnivariateDistribution produces an distribution function for input into cdf, pdf, rand, etc.

Backmatter

Useful Packages on Julia Repository

Notes

  • EmpiricalUnivariateDistribution is broken.
snippet.julianoeval
[download only julia statements]
julia> using Random, Distributions, Plots; Random.seed!(0);
 
julia> xr= -4.0:0.01:4.0; plot( xr, cdf.( Normal(0,1), xr ), legend= false) ## plots the standard normal first
 
julia> x=  max.(randn(10000), (rand(10000) .- 0.5));
 
julia> empcdfun= EmpiricalUnivariateDistribution( x );
 
julia> cdf.( empcdfun, [-0.5,-0.25,0.0,0.5, 10.0] )
5-element Array{Float64,1}:
 0.0
 0.0998
 0.2453
 0.6972
 1.0
 
julia> rand( empcdfun, 4 )
4-element Array{Float64,1}:
 0.43878137920825677
 0.17368409076394875
 0.17510611936949205
 1.0969671028632133
 
julia> plot!( x -> cdf(empcdfun, x) )

References

random.txt · Last modified: 2018/12/05 20:51 (external edit)