defmodule Chi2fit.MonteCarlo do
# Copyright 2016-2021 Pieter Rijken
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import Chi2fit.Statistics
alias Chi2fit.Matrix, as: M
@doc """
Basic Monte Carlo simulation to repeatedly run a simulation multiple times.
## Options
`:collect_all?` - If true, collects data from each individual simulation run and returns this an the third element of the result tuple
"""
@spec mc(iterations :: pos_integer, fun :: ((pos_integer) -> float), options :: Keyword.t) :: {avg :: float, sd :: float, tries :: [float]} | {avg :: float, sd :: float}
def mc(iterations,fun,options \\ []) do
all? = options[:collect_all?] || false
tries = 1..iterations |> Enum.map(fn _ -> fun.() end)
avg = moment tries, 1
sd = :math.sqrt momentc(tries,2,avg)
if all?, do: {avg,sd,tries}, else: {avg,sd}
end
def total_mc(result, fun, mode \\ :use_bounds, iterations \\ 1000) do
{_, cov, parameters, _} = result
ranges = case mode do
:use_bounds ->
# Pick up the error in the paramater value
errors = cov |> M.diagonal |> Enum.map(fn x -> x|>abs|>:math.sqrt end)
Enum.zip(parameters, errors) |> Enum.map(fn {par, err} -> [par - err, par + err] end)
:use_ranges ->
{_, _, _, parranges} = result
parranges |> Tuple.to_list |> tl
end
outcomes = ranges
|> List.foldr([], fn
[left, right], [] -> [[left],[right]]
[left, right], acc -> Enum.flat_map(acc, & [[left|List.wrap(&1)], [right|List.wrap(&1)]])
end)
|> Enum.map(& mc(iterations, fun.(&1)))
|> Enum.map(& elem(&1,0))
{Enum.min(outcomes), Enum.max(outcomes)}
end
@doc """
Forecasts how many time periods are needed to complete `size` items
Related functions: `forecast_duration/2` and `forecast_items/2`.
"""
@spec forecast(fun :: (() -> non_neg_integer),size :: pos_integer, tries :: pos_integer, update :: (() -> number)) :: number
def forecast(fun, size, tries \\ 0,update \\ fn -> 1 end)
def forecast(fun, size, tries, update) when size>0 do
forecast(fun, size-fun.(),tries+update.(),update)
end
def forecast(_fun,_size,tries,_update), do: tries
@doc """
Returns a function for forecasting the duration to complete a number of items.
This function is a wrapper for `forecast/4`.
## Arguments
`data` - either a data set to base the forecasting on, or a function that returns (random) numbers
`size` - the number of items to complete
"""
@spec forecast_duration(data :: [number] | (()->number), size :: pos_integer) :: (() -> number)
def forecast_duration(data, size) when is_list(data) do
fn -> forecast(fn -> Enum.random(data) end, size) end
end
def forecast_duration(fun, size) when is_function(fun,0) do
fn -> forecast(fun, size) end
end
@doc """
Returns a function for forecasting the number of completed items in a number periods.
This function is a wrapper for `forecast/4`.
## Arguments
`data` - either a data set to base the forecasting on, or a function that returns (random) numbers
`periods` - the number of periods to forecast the number of completed items for
"""
@spec forecast_items(data :: [number] | (()->number), periods :: pos_integer) :: (() -> number)
def forecast_items(data, periods) when is_list(data) do
fn -> forecast(fn -> 1 end, periods, 0, fn -> Enum.random(data) end) end
end
def forecast_items(fun, periods) when is_function(fun,0) do
fn -> forecast(fn -> 1 end, periods, 0, fun) end
end
end