The following is an excerpt from my upcoming book, Learn Julia, now available in early access from Manning. I imagine you holding this book in your hand at your favourite bookstore, and contemplating whether to buy it. “Why should I care about Julia?”, you wonder. “There are hundreds of languages around. The statistical community uses […]
Tag Archives: Julia
#MonthOfJulia Day 26: Statistics
JuliaStats is a meta-project which consolidates various packages related to statistics and machine learning in Julia. Well worth taking a look if you plan on working in this domain.
julia> x = rand(10); julia> mean(x) 0.5287191472784906 julia> std(x) 0.2885446536178459
Julia already has some builtin support for statistical operations, so additional packages are not strictly necessary. However they do increase the scope and ease of possible operations (as we’ll see below).Julia already has some builtin support for statistical operations. Let’s kick off by loading all the packages that we’ll be looking at today.
julia> using StatsBase, StatsFuns, StreamStats
StatsBase
The documentation for StatsBase can be found here. As the package name implies, it provides support for basic statistical operations in Julia.
High level summary statistics are generated by summarystats().
julia> summarystats(x) Summary Stats: Mean: 0.528719 Minimum: 0.064803 1st Quartile: 0.317819 Median: 0.529662 3rd Quartile: 0.649787 Maximum: 0.974760
Weighted versions of the mean, variance and standard deviation are implemented. There’re also geometric and harmonic means.
julia> w = WeightVec(rand(1:10, 10)); # A weight vector. julia> mean(x, w) # Weighted mean. 0.48819933297961043 julia> var(x, w) # Weighted variance. 0.08303843715334995 julia> std(x, w) # Weighted standard deviation. 0.2881639067498738 julia> skewness(x, w) 0.11688162715805048 julia> kurtosis(x, w) -0.9210456851144664 julia> mean_and_std(x, w) (0.48819933297961043,0.2881639067498738)
There’s a weighted median as well as functions for calculating quantiles.
julia> median(x) # Median.
0.5296622773635412
julia> median(x, w) # Weighted median.
0.5729104703595038
julia> quantile(x)
5-element Array{Float64,1}:
0.0648032
0.317819
0.529662
0.649787
0.97476
julia> nquantile(x, 8)
9-element Array{Float64,1}:
0.0648032
0.256172
0.317819
0.465001
0.529662
0.60472
0.649787
0.893513
0.97476
julia> iqr(x) # Inter-quartile range.
0.3319677541313941
Sampling from a population is also catered for, with a range of algorithms which can be applied to the sampling procedure.
julia> sample(['a':'z'], 5) # Sampling (with replacement).
5-element Array{Char,1}:
'w'
'x'
'e'
'e'
'o'
julia> wsample(['T', 'F'], [5, 1], 10) # Weighted sampling (with replacement).
10-element Array{Char,1}:
'F'
'T'
'T'
'T'
'F'
'T'
'T'
'T'
'T'
'T'
There’s also functionality for empirical estimation of distributions from histograms and a range of other interesting and useful goodies.
StatsFuns
The StatsFuns package provides constants and functions for statistical computing. The constants are by no means essential but certainly very handy. Take, for example, twoπ and sqrt2.
There are some mildly exotic mathematical functions available like logistic, logit and softmax.
julia> logistic(-5)
0.0066928509242848554
julia> logistic(5)
0.9933071490757153
julia> logit(0.25)
-1.0986122886681098
julia> logit(0.75)
1.0986122886681096
julia> softmax([1, 3, 2, 5, 3])
5-element Array{Float64,1}:
0.0136809
0.101089
0.0371886
0.746952
0.101089
Finally there is a suite of functions relating to various statistical distributions. The functions for the Normal distribution are illustrated below, but there’re functions for Beta and Binomial distribution, the Gamma and Hypergeometric distribution and many others. The function naming convention is consistent across all distributions.
julia> normpdf(0); # PDF julia> normlogpdf(0); # log PDF julia> normcdf(0); # CDF julia> normccdf(0); # Complementary CDF julia> normlogcdf(0); # log CDF julia> normlogccdf(0); # log Complementary CDF julia> norminvcdf(0.5); # inverse-CDF julia> norminvccdf(0.99); # inverse-Complementary CDF julia> norminvlogcdf(-0.693147180559945); # inverse-log CDF julia> norminvlogccdf(-0.693147180559945); # inverse-log Complementary CDF
StreamStats
Finally, the StreamStats package supports calculating online statistics for a stream of data which is being continuously updated.
julia> average = StreamStats.Mean()
Online Mean
* Mean: 0.000000
* N: 0
julia> variance = StreamStats.Var()
Online Variance
* Variance: NaN
* N: 0
julia> for x in rand(10)
update!(average, x)
update!(variance, x)
@printf("x = %3.f: mean = %.3f | variance = %.3fn", x, state(average),
state(variance))
end
x = 0.928564: mean = 0.929 | variance = NaN
x = 0.087779: mean = 0.508 | variance = 0.353
x = 0.253300: mean = 0.423 | variance = 0.198
x = 0.778306: mean = 0.512 | variance = 0.164
x = 0.566764: mean = 0.523 | variance = 0.123
x = 0.812629: mean = 0.571 | variance = 0.113
x = 0.760074: mean = 0.598 | variance = 0.099
x = 0.328495: mean = 0.564 | variance = 0.094
x = 0.303542: mean = 0.535 | variance = 0.090
x = 0.492716: mean = 0.531 | variance = 0.080
In addition to the mean and variance illustrated above, the package also supports online versions of min() and max(), and can be used to generate incremental confidence intervals for Bernoulli and Poisson processes.
That’s it for today. Check out the full code on github and watch the video below.
The post #MonthOfJulia Day 26: Statistics appeared first on Exegetic Analytics.
#MonthOfJulia Day 25: Interfacing with Other Languages
Julia has native support for calling C and FORTRAN functions. There are also add on packages which provide interfaces to C++, R and Python. We’ll have a brief look at the support for C and R here. Further details on these and the other supported languages can be found on github.
Why would you want to call other languages from within Julia? Here are a couple of reasons:
- to access functionality which is not implemented in Julia;
- to exploit some efficiency associated with another language.
The second reason should apply relatively seldom because, as we saw some time ago, Julia provides performance which rivals native C or FORTRAN code.
C
C functions are called via ccall(), where the name of the C function and the library it lives in are passed as a tuple in the first argument, followed by the return type of the function and the types of the function arguments, and finally the arguments themselves. It’s a bit klunky, but it works!
julia> ccall((:sqrt, "libm"), Float64, (Float64,), 64.0) 8.0
It makes sense to wrap a call like that in a native Julia function.
julia> csqrt(x) = ccall((:sqrt, "libm"), Float64, (Float64,), x); julia> csqrt(64.0) 8.0
This function will not be vectorised by default (just try call csqrt() on a vector!), but it’s a simple matter to produce a vectorised version using the @vectorize_1arg macro.
julia> @vectorize_1arg Real csqrt;
julia> methods(csqrt)
# 4 methods for generic function "csqrt":
csqrt{T<:Real}(::AbstractArray{T<:Real,1}) at operators.jl:359
csqrt{T<:Real}(::AbstractArray{T<:Real,2}) at operators.jl:360
csqrt{T<:Real}(::AbstractArray{T<:Real,N}) at operators.jl:362
csqrt(x) at none:6
Note that a few extra specialised methods have been introduced and now calling csqrt() on a vector works perfectly.
julia> csqrt([1, 4, 9, 16])
4-element Array{Float64,1}:
1.0
2.0
3.0
4.0
R
I’ll freely admit that I don’t dabble in C too often these days. R, on the other hand, is a daily workhorse. So being able to import R functionality into Julia is very appealing. The first thing that we need to do is load up a few packages, the most important of which is RCall. There’s great documentation for the package here.
julia> using RCall julia> using DataArrays, DataFrames
We immediately have access to R’s builtin data sets and we can display them using rprint().
julia> rprint(:HairEyeColor)
, , Sex = Male
Eye
Hair Brown Blue Hazel Green
Black 32 11 10 3
Brown 53 50 25 15
Red 10 10 7 7
Blond 3 30 5 8
, , Sex = Female
Eye
Hair Brown Blue Hazel Green
Black 36 9 5 2
Brown 66 34 29 14
Red 16 7 7 7
Blond 4 64 5 8
We can also copy those data across from R to Julia.
julia> airquality = DataFrame(:airquality); julia> head(airquality) 6x6 DataFrame | Row | Ozone | Solar.R | Wind | Temp | Month | Day | |-----|-------|---------|------|------|-------|-----| | 1 | 41 | 190 | 7.4 | 67 | 5 | 1 | | 2 | 36 | 118 | 8.0 | 72 | 5 | 2 | | 3 | 12 | 149 | 12.6 | 74 | 5 | 3 | | 4 | 18 | 313 | 11.5 | 62 | 5 | 4 | | 5 | NA | NA | 14.3 | 56 | 5 | 5 | | 6 | 28 | NA | 14.9 | 66 | 5 | 6 |
rcopy() provides a high-level interface to function calls in R.
julia> rcopy("runif(3)")
3-element Array{Float64,1}:
0.752226
0.683104
0.290194
However, for some complex objects there is no simple way to translate between R and Julia, and in these cases rcopy() fails. We can see in the case below that the object of class lm returned by lm() does not diffuse intact across the R-Julia membrane.
julia> "fit <- lm(bwt ~ ., data = MASS::birthwt)" |> rcopy ERROR: `rcopy` has no method matching rcopy(::LangSxp) in rcopy at no file in map_to! at abstractarray.jl:1311 in map_to! at abstractarray.jl:1320 in map at abstractarray.jl:1331 in rcopy at /home/colliera/.julia/v0.3/RCall/src/sexp.jl:131 in rcopy at /home/colliera/.julia/v0.3/RCall/src/iface.jl:35 in |> at operators.jl:178
But the call to lm() was successful and we can still look at the results.
julia> rprint(:fit)
Call:
lm(formula = bwt ~ ., data = MASS::birthwt)
Coefficients:
(Intercept) low age lwt race
3612.51 -1131.22 -6.25 1.05 -100.90
smoke ptl ht ui ftv
-174.12 81.34 -181.95 -336.78 -7.58
You can use R to generate plots with either the base functionality or that provided by libraries like ggplot2 or lattice.
julia> reval("plot(1:10)"); # Will pop up a graphics window...
julia> reval("library(ggplot2)");
julia> rprint("ggplot(MASS::birthwt, aes(x = age, y = bwt)) + geom_point() + theme_classic()")
julia> reval("dev.off()") # ... and close the window.
Watch the videos below for some other perspectives on multi-language programming with Julia. Also check out the complete code for today (including examples with C++, FORTRAN and Python) on github.
The post #MonthOfJulia Day 25: Interfacing with Other Languages appeared first on Exegetic Analytics.

