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