diff --git a/docs/Project.toml b/docs/Project.toml
index 766e5c98b..e286bda59 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,7 +1,11 @@
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 
 [compat]
 Documenter = "1"
 GR = "0.72.1, 0.73"
+
+[sources]
+Distributions = {path = ".."}
diff --git a/docs/src/fit.md b/docs/src/fit.md
index b482f995e..68bdfb4bb 100644
--- a/docs/src/fit.md
+++ b/docs/src/fit.md
@@ -15,7 +15,7 @@ This statement fits a distribution of type `D` to a given dataset `x`, where `x`
     `Exponential{Float32}`.  However, in the latter case the type parameter of
     the distribution will be ignored:
 
-    ```julia
+    ```jldoctest; setup = :(using Distributions)
     julia> fit(Cauchy{Float32}, collect(-4:4))
     Cauchy{Float64}(μ=0.0, σ=2.0)
     ```
diff --git a/docs/src/starting.md b/docs/src/starting.md
index 4ef4ff856..f5c86a8ff 100644
--- a/docs/src/starting.md
+++ b/docs/src/starting.md
@@ -11,49 +11,60 @@ We start by drawing 100 observations from a standard-normal random variable.
 
 The first step is to set up the environment:
 
-```julia
+```jldoctest getting-started
 julia> using Random, Distributions
 
-julia> Random.seed!(123) # Setting the seed
+julia> Random.seed!(123); # Setting the seed
 ```
 
 Then, we create a standard-normal distribution `d` and obtain samples using `rand`:
 
-```julia
+```jldoctest getting-started
 julia> d = Normal()
-Normal(μ=0.0, σ=1.0)
+Normal{Float64}(μ=0.0, σ=1.0)
 ```
 
 The object `d` represents a probability distribution, in our case the standard-normal distribution.
 One can query its properties such as the mean:
 
-```julia
+```jldoctest getting-started
 julia> mean(d)
 0.0
 ```
 
 We can also draw samples from `d` with `rand`.
-```julia
+```julia-repl
 julia> x = rand(d, 100)
-100-element Array{Float64,1}:
-  0.376264
- -0.405272
- ...
+100-element Vector{Float64}:
+  0.8082879284649668
+ -1.1220725081141734
+ -1.1046361023292959
+ -0.4169926351649334
+  0.28758798062385577
+  0.2298186980518676
+ -0.4217686643996927
+  ⋮
+  0.4350014776438522
+  0.8402951127287839
+ -1.088218513936287
+  0.7037583257923017
+  0.14332589323751366
+  0.14837536667608195
 ```
 
 You can easily obtain the `pdf`, `cdf`, `quantile`, and many other functions for a distribution. For instance, the median (50th percentile) and the 95th percentile for the standard-normal distribution are given by:
 
-```julia
+```jldoctest getting-started
 julia> quantile.(Normal(), [0.5, 0.95])
-2-element Array{Float64,1}:
+2-element Vector{Float64}:
  0.0
- 1.64485
+ 1.6448536269514717
 ```
 
 The normal distribution is parameterized by its mean and standard deviation. To draw random samples from a normal distribution with mean 1 and standard deviation 2, you write:
 
-```julia
-julia> rand(Normal(1, 2), 100)
+```jldoctest getting-started
+julia> rand(Normal(1, 2), 100);
 ```
 
 ## Using Other Distributions
@@ -83,11 +94,9 @@ julia> truncated(Normal(mu, sigma), l, u)
 
 To find out which parameters are appropriate for a given distribution `D`, you can use `fieldnames(D)`:
 
-```julia
+```jldoctest getting-started
 julia> fieldnames(Cauchy)
-2-element Array{Symbol,1}:
- :μ
- :β
+(:μ, :σ)
 ```
 
 This tells you that a Cauchy distribution is initialized with location `μ` and scale `β`.
@@ -96,9 +105,9 @@ This tells you that a Cauchy distribution is initialized with location `μ` and
 
 It is often useful to approximate an empirical distribution with a theoretical distribution. As an example, we can use the array `x` we created above and ask which normal distribution best describes it:
 
-```julia
+```julia-repl
 julia> fit(Normal, x)
-Normal(μ=0.036692077201688635, σ=1.1228280164716382)
+Normal{Float64}(μ=-0.04827714875398303, σ=0.9256810813636542)
 ```
 
 Since `x` is a random draw from `Normal`, it's easy to check that the fitted values are sensible. Indeed, the estimates [0.04, 1.12] are close to the true values of [0.0, 1.0] that we used to generate `x`.