update docs

Examples/Optimize_Example.jl

@@ -8,7 +8,7 @@ function sample_prior()
     return [rand(Uniform(-5, 5), 2)]
  end
 
- function rastrigin(data, x)
+ function rastrigin(_, x)
      A = 10.0
      n = length(x)
      y = A * n

docs/Project.toml

@@ -4,5 +4,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"

docs/make.jl

@@ -15,7 +15,8 @@ makedocs(
     modules = [DifferentialEvolutionMCMC],
     pages = ["home" => "index.md",
     "examples" => ["Binomial Model" => "binomial.md",
-                "Gaussian Model" => "gaussian.md"],
+                "Gaussian Model" => "gaussian.md",
+                "Optimization" => "optimization.md"],
     "api" => "api.md"]
 )
 

docs/src/binomial.md

@@ -123,7 +123,7 @@ de = DE(;sample_prior, bounds, burnin=1000, Np=3, σ=.01)
 ```
 
 ## Estimate Parameter
-The code block below runs the sampler for $2000$ trials with each group of particles running on a separate thread. The progress bar is also set to display. 
+The code block below runs the sampler for $2000$ iterations with each group of particles running on a separate thread. The progress bar is also set to display. 
 ```@example binomial_setup
 n_iter = 2000
 chains = sample(model, de, MCMCThreads(), n_iter, progress=true)

docs/src/gaussian.md

@@ -124,7 +124,7 @@ de = DE(;sample_prior, bounds, burnin=1000, Np=6)
 ```
 
 ## Estimate Parameter
-The code block below runs the sampler for $2000$ trials with each group of particles running on a separate thread. The progress bar is also set to display. 
+The code block below runs the sampler for $2000$ iterations with each group of particles running on a separate thread. The progress bar is also set to display. 
 ```@example gaussian_setup
 n_iter = 2000
 chains = sample(model, de, MCMCThreads(), n_iter, progress=true)

docs/src/index.md

@@ -1,15 +1,16 @@
 # DifferentialEvolutionMCMC.jl
 
-Documentation for DifferentialEvolutionMCMC is under construction.
-
 ```@setup de_animation 
 using DifferentialEvolutionMCMC
 using DifferentialEvolutionMCMC: sample_init
 using DifferentialEvolutionMCMC: crossover!
 using Distributions
 using StatsPlots
+#using PyPlot
 using Random
 
+#pyplot()
+
 Random.seed!(81872)
 
 data = rand(Normal(0.0, 1.0), 5)
@@ -64,26 +65,26 @@ animation = @animate for i in 1:n_iter
 end
 gif(animation, "de_animation.gif", fps = 4)
 ```
-This Julia package provides tools for  performing Bayesian parameter estimation using Differential Evolution MCMC (DEMCMC). You can find several examples with the navigation panel on the left. 
+Welcome to DifferentialEvolutionMCMC.jl. With this package, you can perform Bayesian parameter estimation using Differential Evolution MCMC (DEMCMC), and perform optimization using the basic DE algorithm  Please see the navigation panel on the left for information pertaining to the API and runnable examples. 
 
 ## How Does it Work?
 ### Intuition 
 
-The basic idea behind DEMCMC is that a group of $P$ interacting particles traverse the parameter space and share information about the joint posterior distribution of model parameters. Across many iterations, the samples obtained from the particles will approximate the posterior distribution. The image below illustrates how the particles sample from the posterior distribution of $\mu$ and $\sigma$ of a simple Gaussian model. In this example, five observations were sampled from a Gaussian distribution with $\mu=0$ and $\sigma=1.0$.  The DEMCMC sampler consists of four color coded groups of particles which operate semi-independently from each other. Note that the dashed lines represent the maximum likelihood estimates. The particles cluster near the maximum likelihood estimates because the true parameters are close the center of the prior distributions. 
+The basic idea behind DEMCMC is that a group of interacting particles traverse the parameter space and share information about the joint posterior distribution of model parameters. Across many iterations, the samples obtained from the particles will approximate the posterior distribution. The image below illustrates how the particles sample from the posterior distribution of $\mu$ and $\sigma$ of a simple Gaussian model. In this example, five observations were sampled from a Gaussian distribution with $\mu=0$ and $\sigma=1.0$. The DEMCMC sampler consists of four color coded groups of particles which operate semi-independently from each other. Note that the dashed lines represent the maximum likelihood estimates. The particles cluster near the maximum likelihood estimates because the true parameters are close the center of the prior distributions. 
 
 ![](de_animation.gif)
 
 ### Technical Description 
 
-This section provides a more technical explanation of the basic algorithm. Please see the references below for more details. More formally, a particle $p \in [1,2,\dots, P]$ is a vector of parameters in parameter space defined as:
+This section provides a more technical explanation of the basic algorithm. Please see the references below for more details. More formally, a particle $p \in [1,2,\dots, P]$ is a vector of $n$ parameters in a $\mathbb{R}^n$ parameter space defined as:
 
 $\Theta_p = [\theta_{p,1},\theta_{p,2},\dots \theta_{p,n}].$ 
 
 On each iteration $i$, a new position for each particle $p$ is proposed by adding the weighted difference of two randomly selected particles $j,k$ to particle $p$ along with a small amount of noise. Formally, the proposal is given by:
 
 $\Theta_p^\prime = \Theta_p + \gamma (\Theta_j - \Theta_k) + b,$
 
-where $b \sim \mathrm{uniform}(-\epsilon, \epsilon)$. DEMCMC uses the difference between randomly selected particles to leverage approximate derivatives in the proposal process. The proposal is accepted according to the Metropolis-Hastings rule whereby the proposal is always accepted if its likelihood is greater than the current position, but is accepted proportionally to the ratio of likelihoods otherwise.
+where $b \sim \mathrm{uniform}(-\epsilon, \epsilon)$. DEMCMC uses the difference between randomly selected particles to leverage approximate derivatives in the proposal process. The proposal is accepted according to the Metropolis-Hastings rule whereby the proposal is always accepted if its log likelihood is greater than that of the current position, but is accepted proportionally to the ratio of log likelihoods otherwise.
 # References
 
 Ter Braak, C. J. (2006). A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces. Statistics and Computing, 16, 239-249.

docs/src/optimization.md

@@ -0,0 +1,116 @@
+```@setup optimization_setup
+using DifferentialEvolutionMCMC
+using DifferentialEvolutionMCMC: minimize!
+using Random 
+Random.seed!(6845)
+
+function sample_prior()
+    return [rand(Uniform(-5, 5), 2)]
+ end
+
+ function rastrigin(_, x)
+     A = 10.0
+     n = length(x)
+     y = A * n
+     for  i in 1:n
+         y +=  + x[i]^2 - A * cos(2 * π * x[i])
+     end
+     return y 
+ end
+
+bounds = ((-5.0,5.0),)
+names = (:x,)
+
+model = DEModel(; 
+    sample_prior, 
+    loglike = rastrigin, 
+    data = nothing,
+    names)
+
+de = DE(;
+    sample_prior,
+    bounds,
+    Np = 10, 
+    n_groups = 1, 
+    update_particle! = minimize!,
+    evaluate_fitness! = evaluate_fun!)
+```
+# Optimization Example
+
+The purpose of this example is to demonstrate how to optimize a function as opposed to perform Bayesian parameter estimation.
+## Load Packages
+Our first step is to load the required packages.
+
+```@example optimization_setup
+using DifferentialEvolutionMCMC
+using DifferentialEvolutionMCMC: minimize!
+using Distributions
+using Random 
+Random.seed!(6845)
+```
+
+## Initialize the DE Algorithm
+Each particle in DifferentialEvolution must be seeded with an initial value. To do so, we define a function that returns the initial value. Even though we are not performing Bayesian parameter estimation, we want to use prior information to intelligently seed the DE algorithm. In our case, we will return a vector contained in a vector. 
+```@example optimization_setup
+function sample_prior()
+    return [rand(Uniform(-5, 5), 2)]
+ end
+```
+
+## Objective Function
+In this example, we will find the minimum of the [rastrigin](https://en.wikipedia.org/wiki/Rastrigin_function), which challenging due to its "egg carton-like" surface containing several local minima. 
+
+![](https://upload.wikimedia.org/wikipedia/commons/8/8b/Rastrigin_function.png)
+
+The minimum of the rastrigin function is the zero vector for input `x`. The code block below defines the rastrigin function for an input vector of an arbitrary length. Since we will not be passing data to the function, we will set the first argument to `_` and assign it a value of `nothing` below.
+
+```@example optimization_setup
+function rastrigin(_, x)
+    A = 10.0
+    n = length(x)
+    y = A * n
+    for  i in 1:n
+        y +=  + x[i]^2 - A * cos(2 * π * x[i])
+    end
+    return y 
+end
+```
+## Define Bounds and Parameter Names
+We must define the lower and upper bounds of each parameter. The bounds are a `Tuple` of tuples, where the $i^{\mathrm{th}}$ inner tuple contains the lower and upper bounds of the $i^{\mathrm{th}}$ parameter. In the `rastrigin` function, `x` is a vector with an unspecified length. The bounds below applies to each element in `x`.
+```example optimization_setup
+bounds = ((-5.0,5.0),)
+```
+We can also pass a `Tuple` of names for the argument `x`.
+
+```@example optimization_setup
+names = (:x,)
+```
+## Define Model Object
+Now that we have defined the necessary components of our model, we will organize them into a `DEModel` object as shown below. Notice that `data = nothing` because we do not need data for this optimization problem.
+
+```@example optimization_setup
+model = DEModel(; 
+    sample_prior, 
+    loglike = rastrigin, 
+    data = nothing,
+    names)
+```
+## Define Sampler
+Next, we will create a sampler with the constructor `DE`. Here, we will pass the `sample_prior` function and the variable `bounds`, which constains the lower and upper bounds of each parameter. In addition, we will specify $1,000$ burnin samples, `Np=10` particles for `n_groups=1` groups. In the last two keyword arguments, we use `minimize!` for the `update_particle!` function because we want to minimize `rastrigin`, and `evaluate_fitness! = evaluate_fun!` because we do not want to incorporate the prior log likelihood into our evaluation of `rastrigin`.
+```@example optimization_setup
+de = DE(;
+    sample_prior,
+    bounds,
+    Np = 10, 
+    n_groups = 1, 
+    update_particle! = minimize!,
+    evaluate_fitness! = evaluate_fun!)
+```
+
+## Optimize the Function
+The code block below runs the optimizer for $2000$ iterations with the progress bar. 
+```@example optimization_setup
+n_iter = 2000
+particles = optimize(model, de, n_iter, progress=true);
+results = get_optimal(de, model, particles)
+```