commit 00464bace532f7c893d82a69f366d0b7689f5665 Author: François Pelletier Date: Thu Jan 2 17:11:50 2020 -0500 first commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f10a7aa --- /dev/null +++ b/.gitignore @@ -0,0 +1,12 @@ +pom.xml +pom.xml.asc +*jar +/lib/ +/classes/ +/target/ +/checkouts/ +.lein-deps-sum +.lein-repl-history +.lein-plugins/ +.lein-failures +/bin/ diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..3260e4f --- /dev/null +++ b/LICENSE @@ -0,0 +1,204 @@ +Eclipse Public License - v 1.0 + +THE ACCOMPANYING PROGRAM IS PROVIDED UNDER THE TERMS OF THIS ECLIPSE PUBLIC +LICENSE ("AGREEMENT"). ANY USE, REPRODUCTION OR DISTRIBUTION OF THE PROGRAM +CONSTITUTES RECIPIENT'S ACCEPTANCE OF THIS AGREEMENT. + +1. DEFINITIONS + +"Contribution" means: + +a) in the case of the initial Contributor, the initial code and documentation + distributed under this Agreement, and +b) in the case of each subsequent Contributor: + i) changes to the Program, and + ii) additions to the Program; + + where such changes and/or additions to the Program originate from and are + distributed by that particular Contributor. A Contribution 'originates' + from a Contributor if it was added to the Program by such Contributor + itself or anyone acting on such Contributor's behalf. Contributions do not + include additions to the Program which: (i) are separate modules of + software distributed in conjunction with the Program under their own + license agreement, and (ii) are not derivative works of the Program. + +"Contributor" means any person or entity that distributes the Program. + +"Licensed Patents" mean patent claims licensable by a Contributor which are +necessarily infringed by the use or sale of its Contribution alone or when +combined with the Program. + +"Program" means the Contributions distributed in accordance with this +Agreement. + +"Recipient" means anyone who receives the Program under this Agreement, +including all Contributors. + +2. GRANT OF RIGHTS + a) Subject to the terms of this Agreement, each Contributor hereby grants + Recipient a non-exclusive, worldwide, royalty-free copyright license to + reproduce, prepare derivative works of, publicly display, publicly + perform, distribute and sublicense the Contribution of such Contributor, + if any, and such derivative works, in source code and object code form. + b) Subject to the terms of this Agreement, each Contributor hereby grants + Recipient a non-exclusive, worldwide, royalty-free patent license under + Licensed Patents to make, use, sell, offer to sell, import and otherwise + transfer the Contribution of such Contributor, if any, in source code and + object code form. This patent license shall apply to the combination of + the Contribution and the Program if, at the time the Contribution is + added by the Contributor, such addition of the Contribution causes such + combination to be covered by the Licensed Patents. The patent license + shall not apply to any other combinations which include the Contribution. + No hardware per se is licensed hereunder. + c) Recipient understands that although each Contributor grants the licenses + to its Contributions set forth herein, no assurances are provided by any + Contributor that the Program does not infringe the patent or other + intellectual property rights of any other entity. Each Contributor + disclaims any liability to Recipient for claims brought by any other + entity based on infringement of intellectual property rights or + otherwise. As a condition to exercising the rights and licenses granted + hereunder, each Recipient hereby assumes sole responsibility to secure + any other intellectual property rights needed, if any. For example, if a + third party patent license is required to allow Recipient to distribute + the Program, it is Recipient's responsibility to acquire that license + before distributing the Program. + d) Each Contributor represents that to its knowledge it has sufficient + copyright rights in its Contribution, if any, to grant the copyright + license set forth in this Agreement. + +3. REQUIREMENTS + +A Contributor may choose to distribute the Program in object code form under +its own license agreement, provided that: + + a) it complies with the terms and conditions of this Agreement; and + b) its license agreement: + i) effectively disclaims on behalf of all Contributors all warranties + and conditions, express and implied, including warranties or + conditions of title and non-infringement, and implied warranties or + conditions of merchantability and fitness for a particular purpose; + ii) effectively excludes on behalf of all Contributors all liability for + damages, including direct, indirect, special, incidental and + consequential damages, such as lost profits; + iii) states that any provisions which differ from this Agreement are + offered by that Contributor alone and not by any other party; and + iv) states that source code for the Program is available from such + Contributor, and informs licensees how to obtain it in a reasonable + manner on or through a medium customarily used for software exchange. + +When the Program is made available in source code form: + + a) it must be made available under this Agreement; and + b) a copy of this Agreement must be included with each copy of the Program. + Contributors may not remove or alter any copyright notices contained + within the Program. + +Each Contributor must identify itself as the originator of its Contribution, +if +any, in a manner that reasonably allows subsequent Recipients to identify the +originator of the Contribution. + +4. COMMERCIAL DISTRIBUTION + +Commercial distributors of software may accept certain responsibilities with +respect to end users, business partners and the like. While this license is +intended to facilitate the commercial use of the Program, the Contributor who +includes the Program in a commercial product offering should do so in a manner +which does not create potential liability for other Contributors. Therefore, +if a Contributor includes the Program in a commercial product offering, such +Contributor ("Commercial Contributor") hereby agrees to defend and indemnify +every other Contributor ("Indemnified Contributor") against any losses, +damages and costs (collectively "Losses") arising from claims, lawsuits and +other legal actions brought by a third party against the Indemnified +Contributor to the extent caused by the acts or omissions of such Commercial +Contributor in connection with its distribution of the Program in a commercial +product offering. The obligations in this section do not apply to any claims +or Losses relating to any actual or alleged intellectual property +infringement. In order to qualify, an Indemnified Contributor must: +a) promptly notify the Commercial Contributor in writing of such claim, and +b) allow the Commercial Contributor to control, and cooperate with the +Commercial Contributor in, the defense and any related settlement +negotiations. The Indemnified Contributor may participate in any such claim at +its own expense. + +For example, a Contributor might include the Program in a commercial product +offering, Product X. That Contributor is then a Commercial Contributor. If +that Commercial Contributor then makes performance claims, or offers +warranties related to Product X, those performance claims and warranties are +such Commercial Contributor's responsibility alone. Under this section, the +Commercial Contributor would have to defend claims against the other +Contributors related to those performance claims and warranties, and if a +court requires any other Contributor to pay any damages as a result, the +Commercial Contributor must pay those damages. + +5. NO WARRANTY + +EXCEPT AS EXPRESSLY SET FORTH IN THIS AGREEMENT, THE PROGRAM IS PROVIDED ON AN +"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, EITHER EXPRESS OR +IMPLIED INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OR CONDITIONS OF TITLE, +NON-INFRINGEMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Each +Recipient is solely responsible for determining the appropriateness of using +and distributing the Program and assumes all risks associated with its +exercise of rights under this Agreement , including but not limited to the +risks and costs of program errors, compliance with applicable laws, damage to +or loss of data, programs or equipment, and unavailability or interruption of +operations. + +6. DISCLAIMER OF LIABILITY + +EXCEPT AS EXPRESSLY SET FORTH IN THIS AGREEMENT, NEITHER RECIPIENT NOR ANY +CONTRIBUTORS SHALL HAVE ANY LIABILITY FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING WITHOUT LIMITATION +LOST PROFITS), HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OR DISTRIBUTION OF THE PROGRAM OR THE +EXERCISE OF ANY RIGHTS GRANTED HEREUNDER, EVEN IF ADVISED OF THE POSSIBILITY +OF SUCH DAMAGES. + +7. GENERAL + +If any provision of this Agreement is invalid or unenforceable under +applicable law, it shall not affect the validity or enforceability of the +remainder of the terms of this Agreement, and without further action by the +parties hereto, such provision shall be reformed to the minimum extent +necessary to make such provision valid and enforceable. + +If Recipient institutes patent litigation against any entity (including a +cross-claim or counterclaim in a lawsuit) alleging that the Program itself +(excluding combinations of the Program with other software or hardware) +infringes such Recipient's patent(s), then such Recipient's rights granted +under Section 2(b) shall terminate as of the date such litigation is filed. + +All Recipient's rights under this Agreement shall terminate if it fails to +comply with any of the material terms or conditions of this Agreement and does +not cure such failure in a reasonable period of time after becoming aware of +such noncompliance. If all Recipient's rights under this Agreement terminate, +Recipient agrees to cease use and distribution of the Program as soon as +reasonably practicable. However, Recipient's obligations under this Agreement +and any licenses granted by Recipient relating to the Program shall continue +and survive. + +Everyone is permitted to copy and distribute copies of this Agreement, but in +order to avoid inconsistency the Agreement is copyrighted and may only be +modified in the following manner. The Agreement Steward reserves the right to +publish new versions (including revisions) of this Agreement from time to +time. No one other than the Agreement Steward has the right to modify this +Agreement. The Eclipse Foundation is the initial Agreement Steward. The +Eclipse Foundation may assign the responsibility to serve as the Agreement +Steward to a suitable separate entity. Each new version of the Agreement will +be given a distinguishing version number. The Program (including +Contributions) may always be distributed subject to the version of the +Agreement under which it was received. In addition, after a new version of the +Agreement is published, Contributor may elect to distribute the Program +(including its Contributions) under the new version. Except as expressly +stated in Sections 2(a) and 2(b) above, Recipient receives no rights or +licenses to the intellectual property of any Contributor under this Agreement, +whether expressly, by implication, estoppel or otherwise. All rights in the +Program not expressly granted under this Agreement are reserved. + +This Agreement is governed by the laws of the State of New York and the +intellectual property laws of the United States of America. No party to this +Agreement will bring a legal action under this Agreement more than one year +after the cause of action arose. Each party waives its rights to a jury trial in +any resulting litigation. + diff --git a/README b/README new file mode 100644 index 0000000..9512e80 --- /dev/null +++ b/README @@ -0,0 +1,13 @@ +# actuarial + +Probability distributions for actuarial models + +## Usage + +Use as a function library + +## License + +Copyright (C) 2014 François Pelletier + +Distributed under the Eclipse Public License, the same as Clojure. diff --git a/project.clj b/project.clj new file mode 100644 index 0000000..9845065 --- /dev/null +++ b/project.clj @@ -0,0 +1,7 @@ +(defproject actuarial "1.0.0-SNAPSHOT" + :description "Actuarial functions for Clojure using Incanter and Weka" + :dependencies [[org.clojure/clojure "1.6.0"] + [weka "3.6.2"] + [incanter/incanter-core "1.4.1"] + [ring "1.3.2"] + [org.clojure/data.xml "0.0.8"]]) \ No newline at end of file diff --git a/src/actuarial/basic.clj b/src/actuarial/basic.clj new file mode 100644 index 0000000..4bb4573 --- /dev/null +++ b/src/actuarial/basic.clj @@ -0,0 +1,24 @@ +(ns actuarial.basic) +(import [weka.core Statistics]) + +(defn normalcdf "Normal CDF" [a] (Statistics/normalProbability a)) +(defn normalvar "Normal VAR" [a] (Statistics/normalInverse a)) +(def infinity Double/POSITIVE_INFINITY) + +(defn binomial-coefficient [n k] + (let [rprod (fn [a b] (reduce * (range a (inc b))))] + (/ (rprod (- n k -1) n) (rprod 1 k)))) + +;; from https://gist.github.com/the-max0r/d5affa7c3c8a481733e5 +(defn rootNP + "Newton Raphson univariate root" + [f x] + (def step 1E-4) + (def eps 1E-8) + (let [m (/ (- (f (+ x step)) (f x)) step) + b (- (f x) (* m x)) + x-0 (- (/ b m))] + (if (<= (Math/abs (- x x-0)) eps) + x-0 + (recur f x-0)) + )) \ No newline at end of file diff --git a/src/actuarial/core.clj b/src/actuarial/core.clj new file mode 100644 index 0000000..98682eb --- /dev/null +++ b/src/actuarial/core.clj @@ -0,0 +1,5 @@ +(ns actuarial.core (:require [actuarial.distributions] + [actuarial.mixed-distributions] + [actuarial.basic] + [actuarial.measures] + [actuarial.http])) diff --git a/src/actuarial/distributions.clj b/src/actuarial/distributions.clj new file mode 100644 index 0000000..1cb3f7f --- /dev/null +++ b/src/actuarial/distributions.clj @@ -0,0 +1,1270 @@ +(ns actuarial.distributions (:use actuarial.basic)) +(import [weka.core Statistics]) +(use '(incanter core optimize)) + +;; Methods + +(defmulti pdf "Probability density function" :Distribution) +(defmulti pmf "Probability mass function" :Distribution) +(defmulti survival "Survival function" :Distribution) +(defmulti cdf "Cumulative distribution function" :Distribution) +(defmulti rawmoments "Raw moments function" :Distribution) +(defmulti mean "mean function" :Distribution) +(defmulti variance "Variance function" :Distribution) +(defmulti standard-deviation "Standard deviation function" :Distribution) +(defmulti mgf "Moment generating function" :Distribution) +(defmulti pgf "Probability generating function" :Distribution) +(defmulti truncated-mean "Truncated mean function" :Distribution) +(defmulti value-at-risk "Value at risk function" :Distribution) +(defmulti tail-value-at-risk "Tail value at risk function" :Distribution) +(defmulti stop-loss "Stop-loss function" :Distribution) +(defmulti mean-excess "mean excess function" :Distribution) +(defmulti limited-mean "Limited mean function" :Distribution) +(defmulti standardize "Standardize a value to zero mean and unit variance" :Distribution) +(defmulti describe "Describe a distribution" :Distribution) + +;; ## Default methods + +(defmethod survival :default [f x] + (- 1 (cdf f x))) + +(defmethod stop-loss :default [p d] + (- (truncated-mean p d) + (* d (survival p d)))) + +(defmethod limited-mean :default [f d] + (- (mean f) (stop-loss f d))) + +(defmethod mean-excess :default [f d] + (/ (stop-loss f d) (survival f d))) + +(defmethod mean :default [f] (rawmoments f 1)) + +(defmethod variance :default [f] + (- (rawmoments f 2) (pow (rawmoments f 1) 2))) + +(defmethod standard-deviation :default [f] + (sqrt (variance f))) + +(defmethod value-at-risk :default [f kappa] + {:pre [(<= 0 kappa 1)]} + (case kappa + 0 (get-in f [:Domain :min]) + 1 (get-in f [:Domain :max]) + (rootNP #(- (cdf f %) kappa) (mean f)))) + +(defmethod tail-value-at-risk :default [f kappa] + {:pre [(<= 0 kappa 1)]} + (let [d (value-at-risk f kappa)] + (case kappa + 0 (get-in f [:Domain :min]) + 1 (get-in f [:Domain :max]) + (+ d (/ (stop-loss f d) (- 1 kappa)))))) + +(defmethod standardize :default [f x] + (/ (- x (mean f)) (standard-deviation f))) + +(defmethod mgf :default [f t] + (+ 1 (reduce + (map #(/ (* (pow t %) + (rawmoments f %)) + (factorial %)) + (range 1 100))))) + +;; Describe a distribution + +(defn in-domain [f x] + (let [min (get-in f [:Domain :min]) + max (get-in f [:Domain :max])] + (<= min x max) + )) + +(defn describe [f x d t k] + {:distribution f + :x x + :deductible d + :degree t + :kappa k + :pdf (try (pdf f x) (catch Exception e (println "Error in pdf!" e) 0)) + :cdf (try (cdf f x) (catch Exception e (println "Error in cdf!" e) 0)) + :survival (try (survival f x) (catch Exception e (println "Error in survival!" e) 0)) + :mean (try (mean f) (catch Exception e (println "Error in mean!" e) 0)) + :variance (try (variance f) (catch Exception e (println "Error in variance!" e) 0)) + :truncated-mean (try (truncated-mean f d) (catch Exception e (println "Error in truncated-mean!" e) 0)) + ;:value-at-risk (try (map #(value-at-risk f %) k) (catch Exception e (println "Error in value-at-risk !" e) 0)) + ;:tail-value-at-risk (try (map #(tail-value-at-risk f %) k) (catch Exception e (println "Error in tail-value-at-risk!" e) 0)) + :stop-loss (try (stop-loss f d) (catch Exception e (println "Error in stop-loss!" e) 0)) + :mean-excess (try (mean-excess f d) (catch Exception e (println "Error in mean-excess!" e) 0)) + :limited-mean (try (limited-mean f d) (catch Exception e (println "Error in limited-mean!" e) 0)) + }) + +;; ## Discrete distributions + +(defn arithmetic-dist + [h & prob] + {:pre [(> h 0)(= 1.0 (reduce + (vec prob)))]} + (let [n (dec (count (vec prob)))] + {:Distribution ::Arithmetic + :Parameters {:h h :prob prob :n n} + :Domain {:min 0 :max n} + :Type ::Discrete})) + +(defn poisson-dist + [lambda] + {:pre (> lambda 0)} + {:Distribution ::Poisson + :Parameters {:lambda lambda} + :Domain {:min 0 :max infinity } + :Type ::Discrete}) + +(defn bernoulli-dist + [q] + {:pre [(<= 0 q 1)]} + {:Distribution ::Bernoulli + :Parameters {:q q} + :Domain {:min 0 :max 1} + :Type ::Discrete}) + +(defn binomial-dist + [n q] + {:pre [(integer? n) (<= 0 q 1)]} + (if (= n 1) + (bernoulli-dist q) + ({:Distribution ::Binomial + :Parameters {:n n :q q} + :Domain {:min 0 :max n} + :Type ::Discrete}))) + +(defn geometric-dist + [q] + {:pre [(< 0 q 1)]} + {:Distribution ::Geometric + :Parameters {:q q} + :Domain {:min 0 :max infinity} + :Type ::Discrete}) + +(defn negative-binomial-dist + [r q] + {:pre [(<= 0 q 1)]} + (if (= r 1) + (geometric-dist q) + ({:Distribution ::NegativeBinomial + :Parameters {:r r :q q} + :Domain {:min 0 :max infinity} + :Type ::Discrete}))) + +(defn negative-binomial-type2 + [r betap] + {:pre [(> r 0) (> betap 0)]} + (negative-binomial-dist + r (/ 1 (inc betap)))) + +(defn logarithmic-dist + [nu] + {:pre [(< 0 nu 1)]} + {:Distribution ::Logarithmic + :Parameters {:nu nu} + :Domain {:min 1 :max infinity} + :Type ::Discrete}) + +;; ## Continuous distributions +(defn uniform-dist + ([a b] {:pre [(< a b)]} + {:Distribution ::Uniform + :Parameters {:a a :b b} + :Domain {:min :max} + :Type ::Continuous}) + ([b] (uniform-dist 0 b)) + ([] (uniform-dist 0 1))) + +(defn exponential-dist + [betap] + {:pre [(> betap 0)]} + {:Distribution ::Exponential + :Parameters{:beta betap} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn beta-dist + [alpha betap] + {:pre [(> alpha 0) (> betap 0)]} + {:Distribution ::Beta + :Parameters {:alpha alpha :beta betap} + :Domain {:min 0 :max 1} + :Type ::Continuous}) + +(defn gamma-dist + [alpha betap] + {:pre [(> alpha 0) (> betap 0)]} + {:Distribution ::Gamma + :Parameters {:alpha alpha :beta betap :erlang false} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn erlang-dist + [n betap] + {:pre [(> n 0) (> betap 0) (integer? n)]} + (assoc-in (gamma-dist n betap) [:Parameters :erlang] true)) + +(defn weibull-dist + [tau betap] + {:pre [(> tau 0) (> betap 0)]} + {:Distribution ::Weibull + :Parameters {:tau tau :beta betap} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn generalized-erlang-dist + [& betas] + {:pre [(every? #(> % 0) betas) (apply distinct? betas)]} + {:Distribution ::generalized-erlang + :Parameters {:betas betas :n (count (vec betas))} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn gaussian-dist + ([mu sigma] + {:pre [(> sigma 0)]} + {:Distribution ::Gaussian + :Parameters {:mu mu :sigma sigma} + :Domain {:min (- infinity) :max infinity} + :Type ::Continuous}) + ([] (gaussian-dist 0 1))) + +(defn lognormal-dist + ([mu sigma] + {:pre [(> sigma 0)]} + {:Distribution ::Lognormal :Parameters {:mu mu :sigma sigma} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + ([] (lognormal-dist 0 1))) + +(defn inverse-gaussian-dist + [mu betap] + {:pre [(>= betap 0)]} + {:Distribution ::InverseGaussian + :Parameters {:mu mu :beta betap} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn pareto-dist + [alpha lambda] + {:pre [(> alpha 0) (> lambda 0)]} + {:Distribution ::Pareto + :Parameters {:alpha alpha :lambda lambda} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn generalized-pareto-dist + [alpha lambda tau] + {:pre [(> alpha 0) (> lambda 0) (> tau 0)]} + {:Distribution ::GeneralizedPareto + :Parameters {:alpha alpha :lambda lambda :tau tau} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn burr-dist + [alpha lambda tau] + {:pre [(> alpha 0) (> lambda 0) (> tau 0)]} + {:Distribution ::Burr + :Parameters {:alpha alpha :lambda lambda :tau tau} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn log-logistic-dist + [lambda tau] + {:pre (> tau 0)} + {:Distribution ::LogLogistic + :Parameters {:lambda lambda :tau tau} + :Domain {:min 0 :max infinity} + :Type ::Continuous}) + +(defn student-dist + [nu] + {:pre (> nu 0)} + {:Distribution ::Student + :Parameters {:nu nu} + :Domain {:min (- infinity) :max infinity} + :Type ::Continuous}) + +;; map of all distributions +(def distributions-map + {"arithmetic" arithmetic-dist + "poisson" poisson-dist + "uniform" uniform-dist + "exponential" exponential-dist + "beta" beta-dist + "gamma" gamma-dist + "erlang" erlang-dist + "weibull" weibull-dist + "generalized-erlang" generalized-erlang-dist + "gaussian" gaussian-dist + "lognormal" lognormal-dist + "inverse-gaussian" inverse-gaussian-dist + "pareto" pareto-dist + "generalized-pareto" generalized-pareto-dist + "burr" burr-dist + "log-logistic" log-logistic-dist + "student" student-dist + }) + +;; # Discrete distributions + +;; ## Arithmetic distribution + +(defmethod pmf ::Arithmetic [f x] + {:pre (in-domain f x)} + (let [prob (get-in f [:Parameters :prob]) + h (get-in f [:Parameters :h])] + (nth prob (/ x h)))) + +(defmethod cdf ::Arithmetic [f x] + {:pre (in-domain f x)} + (let [prob (get-in f [:Parameters :prob]) + h (get-in f [:Parameters :h])] + (reduce + (take (inc (/ x h)) + (map #(pmf f (* % h)) + (iterate inc 0)))))) + +(defmethod mean ::Arithmetic [f] + (let [prob (get-in f [:Parameters :prob]) + h (get-in f [:Parameters :h]) + n (get-in f [:Parameters :n])] + (reduce + (map #(* % h (nth prob %)) + (range (inc n)))))) + +(defmethod variance ::Arithmetic [f] + (let [prob (get-in f [:Parameters :prob]) + h (get-in f [:Parameters :h]) + n (get-in f [:Parameters :n])] + (reduce + (map #(* (pow (- (* % h) + (mean f)) 2) + (nth prob %)) + (range (inc n)))))) + +(defmethod mgf ::Arithmetic [f t] + (let [prob (get-in f [:Parameters :prob]) + h (get-in f [:Parameters :h]) + n (get-in f [:Parameters :n])] + (reduce + (map #(* (exp (* t h %)) + (nth prob %)) + (range (inc n)))))) + +(defmethod pgf ::Arithmetic [f t] + (let [prob (get-in f [:Parameters :prob]) + h (get-in f [:Parameters :h]) + n (get-in f [:Parameters :n])] + (reduce + (map #(* (pow t (* h %)) + (nth prob %)) + (range (inc n)))))) + +(defmethod truncated-mean ::Arithmetic [f d] + (let [prob (get-in f [:Parameters :prob]) + h (get-in f [:Parameters :h])] + (reduce + (map #(* % h (nth prob %)) + (range (inc d)))))) + +(defmethod value-at-risk ::Arithmetic [f kappa] + (first (filter #(> (cdf f %) kappa) (iterate inc 0)))) + +(defmethod tail-value-at-risk ::Arithmetic [f kappa] + (let [h (get-in f [:Parameters :h]) + k0 (value-at-risk f kappa)] + (/ (+ (- (mean f) + (reduce + (map #(* % h (pmf f (* % h))) + (range 0 (inc k0))))) + (* k0 h (- (cdf f (* k0 h)) kappa))) + (- 1 kappa)))) + +;; ## Poisson distribution + +(defmethod pmf ::Poisson [p k] + {:pre (in-domain p k)} + (let [lambda (get-in p [:Parameters :lambda])] + (/ (* (pow lambda k) + (exp (- lambda))) + (factorial k)))) + +(defmethod cdf ::Poisson [p k] + {:pre (in-domain p k)} + (reduce + (map #(pmf %) (range (inc k))))) + +(defmethod mean ::Poisson [p] + (let [lambda (get-in p [:Parameters :lambda])] + lambda)) + +(defmethod variance ::Poisson [p] + (let [lambda (get-in p [:Parameters :lambda])] + lambda)) + +(defmethod pgf ::Poisson [p t] + (let [lambda (get-in p [:Parameters :lambda])] + (exp (* lambda (- t 1))))) + +(defmethod mgf ::Poisson [p t] + (let [lambda (get-in p [:Parameters :lambda])] + (pgf p (exp t)))) + + +;; ## Bernoulli distribution + +(defmethod pmf ::Bernoulli [b k] + (let [q (get-in b [:Parameters :q])] + (* (pow q k) + (pow (- 1 q) (- 1 k))) + )) + +(defmethod mean ::Bernoulli [b] + (let [q (get-in b [:Parameters :q])] + q)) + +(defmethod variance ::Bernoulli [b] + (let [q (get-in b [:Parameters :q])] + (* q (- 1 q)))) + +(defmethod mgf ::Bernoulli [b t] + (let [q (get-in b [:Parameters :q])] + (+ (* q (exp t)) (- 1 q)))) + +(defmethod pgf ::Bernoulli [b t] + (let [q (get-in b [:Parameters :q])] + (+ (* q t) (- 1 q)))) + +;; ## Binomial distribution + +(defmethod pmf ::Binomial [b k] + (let [n (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (* (binomial-coefficient n k) + (pow q k) + (pow (- 1 q) (- n k))) + )) + +(defmethod mean ::Binomial [b] + (let [n (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (* n q))) + +(defmethod variance ::Binomial [b] + (let [n (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (* n q (- 1 q)))) + +(defmethod mgf ::Binomial [b t] + (let [n (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (pow (+ (* q (exp t)) (- 1 q)) n))) + +(defmethod pgf ::Binomial [b t] + (let [n (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (pow (+ (* q t) (- 1 q)) n))) + +;; ## Negative binomial distribution + +(defmethod pmf ::NegativeBinomial [b k] + (let [r (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (* (binomial-coefficient (+ r k 1) k ) + (pow q r) + (pow (- 1 q) k)))) + +(defmethod mean ::NegativeBinomial [b] + (let [r (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (/ (* r (- 1 q)) q))) + +(defmethod variance ::NegativeBinomial [b] + (let [r (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (/ (* r (- 1 q)) (pow q 2))) + ) + +(defmethod mgf ::NegativeBinomial [b t] + (let [r (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (pow (/ q (- 1 (* (- 1 q) (exp t)))) r) + )) + +(defmethod pgf ::NegativeBinomial [b t] + (let [r (get-in b [:Parameters :n]) + q (get-in b [:Parameters :q])] + (pow (/ q (- 1 (* (- 1 q) t))) r) + )) + +;; ## Geometric distribution + +(defmethod pmf ::Geometric [b k] + (let [q (get-in b [:Parameters :q])] + (* q (- 1 q)))) + +(defmethod cdf ::Geometric [b k] + (let [q (get-in b [:Parameters :q])] + (- 1 (pow (- 1 q) (inc k))))) + +(defmethod mean ::Geometric [b] + (let [q (get-in b [:Parameters :q])] + (/ (- 1 q) q))) + +(defmethod variance ::Geometric [b] + (let [q (get-in b [:Parameters :q])] + (/ (- 1 q) (pow q 2)))) + +(defmethod mgf ::Geometric [b t] + (let [q (get-in b [:Parameters :q])] + (/ q (- 1 (* (- 1 q) (exp t)))))) + +(defmethod pgf ::Geometric [b t] + (let [q (get-in b [:Parameters :q])] + (/ q (- 1 (* (- 1 q) t))))) + +;; ## Logarithmic distribution + +(defmethod pmf ::Logarithmic [b k] + (let [nu (get-in b [:Parameters :nu])] + (/ (- (pow nu k)) (* (log (- 1 nu)) k)))) + +(defmethod mean ::Logarithmic [b] + (let [nu (get-in b [:Parameters :nu])] + (/ (- nu) (* (log (- 1 nu)) (- 1 nu))))) + +(defmethod variance ::Logarithmic [b] + (let [nu (get-in b [:Parameters :nu])] + (/ (+ nu (log (- 1 nu))) + (* (pow (- 1 nu) 2) + (pow (log (- 1 nu)) 2))))) + +(defmethod mgf ::Logarithmic [b t] + (let [nu (get-in b [:Parameters :nu])] + (/ (log (- 1 (* nu (exp t)))) + (log (- 1 nu))))) + +(defmethod pgf ::Logarithmic [b t] + (let [nu (get-in b [:Parameters :nu])] + (/ (log (- 1 (* nu t))) + (log (- 1 nu))))) + +;; # Continuous distributions + +;; ## Uniform distribution + +(defmethod pdf ::Uniform [u x] + {:pre (in-domain u x)} + (let [a (get-in u [:Parameters :a]) + b (get-in u [:Parameters :b])] + (/ 1 (- b a)))) + +(defmethod cdf ::Uniform [u x] + {:pre (in-domain u x)} + (let [a (get-in u [:Parameters :a]) + b (get-in u [:Parameters :b])] + (min (max 0 (/ (- x a) (- b a))) 1))) + +(defmethod rawmoments ::Uniform [u k] + (let [a (get-in u [:Parameters :a]) + b (get-in u [:Parameters :b])] + (/ (- (pow b (inc k)) (pow a (inc k))) + (* (inc k) (- b a))))) + +(defmethod mgf ::Uniform [u t] + (let [a (get-in u [:Parameters :a]) + b (get-in u [:Parameters :b])] + (/ (- (Math/exp (* b t)) + (Math/exp (* a t))) + (* (- b a) t)))) + +(defmethod truncated-mean ::Uniform [u d] + {:pre (< d (get-in u [:Parameters :b]))} + (let [a (get-in u [:Parameters :a]) + b (get-in u [:Parameters :b])] + (/ (- (pow d 2) (pow a 2)) + (* 2 (- b a))))) + +(defmethod value-at-risk ::Uniform [u kappa] + {:pre [(<= 0 kappa 1)]} + (let [a (get-in u [:Parameters :a]) + b (get-in u [:Parameters :b])] + (+ a (* (- b a) kappa)))) + +;; ## Exponential distribution + +(defmethod pdf ::Exponential [e x] + {:pre (in-domain e x)} + (let [betap (get-in e [:Parameters :beta])] + (* betap (exp (- (* betap x)))))) + +(defmethod cdf ::Exponential [e x] + {:pre (in-domain e x)} + (let [betap (get-in e [:Parameters :beta])] + (- 1 (exp (- (* betap x)))))) + +(defmethod mean ::Exponential [e] + (let [betap (get-in e [:Parameters :beta])] + (/ 1 betap))) + +(defmethod variance ::Exponential [e] + (let [betap (get-in e [:Parameters :beta])] + (/ 1 (pow betap 2)))) + +(defmethod mgf ::Exponential [e t] + {:pre [(< t (get-in e [:Parameters :beta]))]} + (let [betap (get-in e [:Parameters :beta])] + (/ betap (- betap t)))) + +(defmethod rawmoments ::Exponential [e k] + (let [betap (get-in e [:Parameters :beta])] + (* (pow (/ 1 betap) k) (factorial k)))) + +(defmethod stop-loss ::Exponential [e d] + (* (mean e) (survival e d))) + +(defmethod truncated-mean ::Exponential [e d] + (- (* (limited-mean e d)) + (* d (survival e d)))) + +(defmethod value-at-risk ::Exponential [e kappa] + {:pre [(<= 0 kappa 1)]} + (let [betap (get-in e [:Parameters :beta])] + (- (* (/ 1 betap) + (log (- 1 kappa)))))) + +;; ## Beta functions + +(defn beta-complete + "Complete beta function" + [alpha betap] + {:pre [(> alpha 0) (> betap 0)]} + (/ (* (gamma alpha) + (gamma betap)) + (gamma (+ alpha betap)))) + +(defn beta-incomplete + "Incomplete beta function" + [alpha betap x] + {:pre [(> alpha 0) (> betap 0)]} + (* (Statistics/incompleteBeta alpha betap x) + (beta-complete alpha betap))) + +;; ## Beta distribution + +(defmethod pdf ::Beta [b x] + {:pre (in-domain b x)} + (let [alpha (get-in b [:Parameters :alpha]) + betap (get-in b [:Parameters :beta])] + (* (/ 1 (beta-complete alpha betap)) + (pow x (- alpha 1)) (pow (- 1 x)(- betap 1))))) + +(defmethod cdf ::Beta [b x] + {:pre (in-domain b x)} + (let [alpha (get-in b [:Parameters :alpha]) + betap (get-in b [:Parameters :beta])] + (defn term-of-sum [alpha betap j] + (* (/ (factorial (+ alpha betap (- 1))) + (* (factorial j) + (factorial (- (+ alpha betap) j 1)))) + (pow x j) + (pow (- 1 x) + (-(+ alpha betap) 1 j)))) + (if (and (integer? alpha) (integer? betap)) + (reduce + (map #(term-of-sum alpha betap %) + (range alpha (+ alpha betap)))) + (/ (beta-incomplete alpha betap x) + (beta-complete alpha betap))))) + +(defn beta-mgf-term-of-product + [alpha betap j] + {:pre [(> alpha 0) (> betap 0)]} + (/ (+ alpha j) + (+ alpha betap j))) + +(defn beta-mgf-innerproduct + [alpha betap k] + {:pre [(> alpha 0) (> betap 0)]} + (reduce * (map #(beta-mgf-term-of-product alpha betap %) + (range k)))) + +(defn beta-mgf-term-of-sum + [alpha betap t k] + {:pre [(> alpha 0) (> betap 0)]} + (* (beta-mgf-innerproduct alpha betap k) + (/ (pow t k) + (factorial k)))) + +(defn beta-mgf-sums [alpha betap t] + {:pre [(> alpha 0) (> betap 0)]} + (map #(beta-mgf-term-of-sum alpha betap t %) + (rest (range)))) + +(defmethod mgf ::Beta [b t] + {:pre []} + (let [alpha (get-in b [:Parameters :alpha]) + betap (get-in b [:Parameters :beta])] + (reduce + (take 100 (beta-mgf-sums alpha betap t))))) + +(defmethod rawmoments ::Beta [b k] + (let [alpha (get-in b [:Parameters :alpha]) + betap (get-in b [:Parameters :beta])] + (/(*(gamma (+ alpha k)) + (gamma (+ alpha betap))) + (*(gamma alpha) + (gamma (+ alpha betap k)))))) + +(defmethod truncated-mean ::Beta [b d] + {:pre [(>= d 0) (<= d 1)]} + (let [alpha (get-in b [:Parameters :alpha]) + betap (get-in b [:Parameters :beta])] + (if (= betap 1) + (/ (* alpha (pow d (+ alpha 1))) + (+ alpha 1)) + (if (= alpha 1) + (+ (- 0 (* d (pow (- 1 d) betap))) + (/ (- 1 (pow (- 1 d) (+ betap 1))) + (+ betap 1))) + (* (/ alpha + (+ alpha betap)) + (cdf (beta-dist (+ alpha 1) betap) d)))))) + +(defmethod value-at-risk ::Beta [b kappa] + {:pre [(<= 0 kappa 1)]} + (let [alpha (get-in b [:Parameters :alpha]) + betap (get-in b [:Parameters :beta])] + (if (= betap 1) + (pow kappa (/ 1 alpha)) + (if (= alpha 1) + (- 1 (pow (- 1 kappa) (/ 1 betap))) + (rootNP #(- (cdf b %) kappa) (mean b)))))) + +;; ## Gamma and Erlang distributions + +(defmethod pdf ::Gamma [g x] + {:pre (in-domain g x)} + (let [alpha (get-in g [:Parameters :alpha]) + betap (get-in g [:Parameters :beta])] + (* (/ (pow betap alpha) + (gamma alpha)) + (pow x (- alpha 1)) + (exp (- (* betap x)))))) + +(defn incompletegamma [s x] + (* (pow x s) + (gamma s) + (exp (- x)) + (reduce + (map #(/ (pow x %) (gamma (+ s 1 %))) (range 100))))) + +(defn upperincompletegamma [s x] + (- (gamma s) + (incompletegamma s x))) + +(defmethod cdf ::Gamma [g x] + {:pre (in-domain g x)} + (let [alpha (get-in g [:Parameters :alpha]) + betap (get-in g [:Parameters :beta])] + (if (integer? alpha) + (- 1 (* (exp (- (* betap x))) + (reduce + (map #(/ (pow (* betap x) %) + (factorial %)) + (range alpha))))) + (/ (incompletegamma alpha (* betap x)) + (gamma alpha))))) + +(defmethod survival ::Gamma [g x] + {:pre (in-domain g x)} + (let [alpha (get-in g [:Parameters :alpha]) + betap (get-in g [:Parameters :beta])] + (if (integer? alpha) + (- 1 (cdf g x)) + (/ (upperincompletegamma + alpha + (* betap x)) + (gamma alpha))))) + + +(defmethod mgf ::Gamma [g t] + {:pre (< t (:beta g))} + (let [alpha (get-in g [:Parameters :alpha]) + betap (get-in g [:Parameters :beta])] + (pow (/ betap(- betap t)) alpha))) + +(defmethod rawmoments ::Gamma [g k] + {:pre (integer? k)} + (let [alpha (get-in g [:Parameters :alpha]) + betap (get-in g [:Parameters :beta])] + (/ (reduce * (map #(+ alpha %) (range k))) + (pow betap k)))) + +(defmethod truncated-mean ::Gamma [g d] + {:pre (> d 0)} + (let [alpha (get-in g [:Parameters :alpha]) + betap (get-in g [:Parameters :beta])] + (* (/ alpha betap) + (cdf (gamma-dist (+ alpha 1) betap) d)))) + +;; ## Auxiliairy distribution for the Weibull method definitions +(defn auxiliairy-gamma-weibull [w] + (let [betap (get-in w [:Parameters :beta]) + tau (get-in w [:Parameters :tau])] + (gamma-dist + (+ 1 (/ 1 tau)) + (pow betap tau)))) + +;; ## Weibull distribution + +(defmethod pdf ::Weibull [w x] + {:pre (in-domain w x)} + (let [betap (get-in w [:Parameters :beta]) + tau (get-in w [:Parameters :tau])] + (* betap + tau + (pow (* betap x) + (- tau 1)) + (exp (- (pow (* betap x) + tau)))))) + +(defmethod cdf ::Weibull [w x] + {:pre (in-domain w x)} + (let [betap (get-in w [:Parameters :beta]) + tau (get-in w [:Parameters :tau])] + (- 1 (exp (- (pow (* betap x) tau)))))) + +(defmethod rawmoments ::Weibull [w k] + (let [betap (get-in w [:Parameters :beta]) + tau (get-in w [:Parameters :tau])] + (/ (gamma (+ 1 (/ k tau))) + (pow betap k)))) + +(defmethod mgf ::Weibull [w k] + (let [betap (get-in w [:Parameters :beta]) + tau (get-in w [:Parameters :tau])] + (reduce + (map #(* (rawmoments w %) + (/ (pow tau %) + (* (pow betap %) + (factorial %)))) + (range 100))))) + +(defmethod truncated-mean ::Weibull [w d] + (let [tau (get-in w [:Parameters :tau])] + (* (mean w) + (cdf (auxiliairy-gamma-weibull w) + (pow d tau))))) + +(defmethod value-at-risk ::Weibull [w kappa] + {:pre [(<= 0 kappa 1)]} + (let [betap (get-in w [:Parameters :beta]) + tau (get-in w [:Parameters :tau])] + (/ (pow (- (log (- 1 kappa))) + (/ 1 tau)) betap))) + +;; ## Generalized Erlang Distribution +;; [reference: R package sdprisk] +;; (htaup://www.inside-r.org/packages/cran/sdprisk/docs/hypoexp) +;; *Helper functions* + +;; Inner big-product +(defn helper-product-generalized-erlang [i betas] + (let [nbeta (nth betas i)] + (reduce * (map #(/ % (- % nbeta)) + (filter #(not (= % nbeta)) betas))))) +;; Big-sum for functions of generalized erlang distribution +(defn helper-big-sum-generalized-erlang [i x betas f] + (* (helper-product-generalized-erlang i betas) + (f i x betas))) + +;; *methods* + +(defmethod pdf ::generalized-erlang [g x] + {:pre (in-domain g x)} + (let [betas (get-in g [:Parameters :betas]) + n (get-in g [:Parameters :n])] + (reduce + (map #(helper-big-sum-generalized-erlang + % x (vec betas) + (fn [i x betas] + (let [nbeta (nth betas i)] + (* nbeta + (exp (- (* x nbeta))))))) + (range n))))) + +(defmethod cdf ::generalized-erlang [g x] + {:pre (in-domain g x)} + (let [betas (get-in g [:Parameters :betas]) + n (get-in g [:Parameters :n])] + (reduce + (map #(helper-big-sum-generalized-erlang + % x (vec betas) + (fn [i x betas] + (- 1 (exp (- (* x (nth betas i))))))) + (range n))))) + +(defmethod mgf ::generalized-erlang [g t] + (let [betas (get-in g [:Parameters :betas])] + (reduce * (map #(/ % (- % t)) betas)))) + +(defmethod rawmoments ::generalized-erlang [g k] + (let [betas (get-in g [:Parameters :betas])] + (reduce * (map #(/ (gamma (+ k 1)) (pow % k)) betas)))) + +(defmethod truncated-mean ::generalized-erlang [g d] + (let [betas (get-in g [:Parameters :betas]) + n (get-in g [:Parameters :n])] + (reduce + (map #(helper-big-sum-generalized-erlang + % d (vec betas) + (fn [i d betas] + (let [nbeta (nth betas i)] + (+ (* (- d) (exp (-(* nbeta d)))) + (/ (- 1 (exp (-(* nbeta d)))) + nbeta))))) + (range n))))) + +;; ## Gaussian Distribution + +(defmethod pdf ::Gaussian [g x] + {:pre (in-domain g x)} + (let [sigma (get-in g [:Parameters :sigma])] + (/ (exp (- (* 0.5 (pow (standardize g x) 2)))) + (* sigma (sqrt (* 2 Math/PI)))))) + +(defmethod cdf ::Gaussian [g x] + {:pre (in-domain g x)} + (normalcdf (standardize g x))) + +(defmethod mean ::Gaussian [g] + (let [mu (get-in g [:Parameters :mu])] + mu)) + +(defmethod variance ::Gaussian [g] + (let [sigma (get-in g [:Parameters :sigma])] + (pow sigma 2))) + +(defmethod mgf ::Gaussian [g t] + (let [mu (get-in g [:Parameters :mu]) + sigma (get-in g [:Parameters :sigma])] + (exp (+ (* t mu) + (* 0.5 (pow t 2) (pow sigma 2)))))) + +(defn sigmasq-times-pdf-helper [g d] + (let [sigma (get-in g [:Parameters :sigma])] + (* (pow sigma 2) (pdf g d)))) + +(defmethod truncated-mean ::Gaussian [g d] + (let [mu (get-in g [:Parameters :mu])] + (- (* mu (cdf g d)) + (sigmasq-times-pdf-helper g d)))) + +(defmethod value-at-risk ::Gaussian [g kappa] + {:pre [(<= 0 kappa 1)]} + (let [mu (get-in g [:Parameters :mu]) + sigma (get-in g [:Parameters :sigma])] + (+ mu + (* sigma + (normalvar kappa))))) + +;; ## Lognormal distribution + +(defmethod standardize ::Lognormal [l x] + (let [mu (get-in l [:Parameters :mu]) + sigma (get-in l [:Parameters :sigma])] + (/ (- (log x) mu (pow sigma 2)) sigma))) + +(defmethod pdf ::Lognormal [l x] + {:pre (in-domain l x)} + (let [mu (get-in l [:Parameters :mu]) + sigma (get-in l [:Parameters :sigma])] + (/ (pdf (gaussian-dist mu sigma) (log x)) x))) + +(defmethod cdf ::Lognormal [l x] + {:pre (in-domain l x)} + (let [mu (get-in l [:Parameters :mu]) + sigma (get-in l [:Parameters :sigma])] + (cdf (gaussian-dist mu sigma) (log x)))) + +(defmethod rawmoments ::Lognormal [l k] + (let [mu (get-in l [:Parameters :mu]) + sigma (get-in l [:Parameters :sigma])] + (exp (+ (* k mu) + (/ (* (pow k 2) + (pow sigma 2)) 2))))) + +(defmethod truncated-mean ::Lognormal [l d] + (let [mu (get-in l [:Parameters :mu]) + sigma (get-in l [:Parameters :sigma])] + (* (exp (+ mu + (/ (pow sigma 2) 2))) + (normalcdf (standardize l d))))) + +(defmethod value-at-risk ::Lognormal [l kappa] + {:pre [(<= 0 kappa 1)]} + (let [mu (get-in l [:Parameters :mu]) + sigma (get-in l [:Parameters :sigma])] + (exp (+ mu + (* sigma + (value-at-risk + (gaussian-dist) + kappa)))))) + +;; ## Inverse gaussian distribution + +(defn e2mb [i] + (let [mu (get-in i [:Parameters :mu]) + betap (get-in i [:Parameters :beta])] + (exp (/ (* 2 mu) betap)))) + +(defn helper1 [i x s] + (let [mu (get-in i [:Parameters :mu]) + betap (get-in i [:Parameters :beta])] + (normalcdf (* (sqrt (/ 1 (* betap x))) (- (* x s) mu))))) + +(defmethod pdf ::InverseGaussian [i x] + {:pre (in-domain i x)} + (let [mu (get-in i [:Parameters :mu]) + betap (get-in i [:Parameters :beta])] + (/ (* mu (exp (- (/ (pow (- x mu) 2) + (* 2 betap x))))) + (sqrt (* 2 Math/PI betap (pow x 3)))))) + +(defmethod cdf ::InverseGaussian [i x] + {:pre (in-domain i x)} + (+ (helper1 i x 1) + (* (e2mb i) + (helper1 i x -1)))) + +(defmethod mean ::InverseGaussian [i] + (let [mu (get-in i [:Parameters :mu])] + mu)) + +(defmethod variance ::InverseGaussian [i] + (let [mu (get-in i [:Parameters :mu]) + betap (get-in i [:Parameters :beta])] + (* mu betap))) + +(defmethod mgf ::InverseGaussian [i t] + (let [mu (get-in i [:Parameters :mu]) + betap (get-in i [:Parameters :beta])] + (exp (* (/ mu beta) + (- 1 (sqrt (- 1 (* 2 beta t)))))))) + +(defmethod truncated-mean ::InverseGaussian [i d] + (let [mu (get-in i [:Parameters :mu])] + (- d + (* (- (* 2 d) mu) (helper1 i d 1)) + (* (+ (* 2 d) mu) + (e2mb i) + (helper1 i d -1))))) + +;; ## Pareto distribution + +(defmethod pdf ::Pareto [p x] + {:pre (in-domain p x)} + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda])] + (/ (* alpha (pow lambda alpha)) + (pow (+ lambda x) (inc alpha))))) + +(defmethod cdf ::Pareto [p x] + {:pre (in-domain p x)} + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda])] + (- 1 (pow (/ lambda (+ lambda x)) alpha)))) + +(defmethod rawmoments ::Pareto [p k] + {:pre (< k (get-in p [:Parameters :alpha]))} + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda])] + (if (and (< k alpha) (integer? k)) + (/ (* (pow lambda k) (factorial k)) + (reduce * (map #(- alpha %) (range 1 (inc k))))) + (if (and (< -1 k) (< k alpha)) + (/ (* (pow lambda k) + (gamma (inc k)) + (gamma (- alpha k))) + (gamma alpha)) + 0)))) + +(defmethod truncated-mean ::Pareto [p d] + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda])] + (- (* (mean p) (cdf (pareto-dist (dec alpha) lambda) d)) + (* d (survival p d))))) + +(defmethod value-at-risk ::Pareto [p kappa] + {:pre [(<= 0 kappa 1)]} + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda])] + (* lambda (dec (pow (- 1 kappa) (- (/ 1 alpha))))))) + +;; ## Generalized Pareto + +(defmethod pdf ::GeneralizedPareto [p x] + {:pre (in-domain p x)} + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda]) + tau (get-in p [:Parameters :tau])] + (/ (* (gamma (+ alpha tau)) + (pow lambda alpha) + (pow x (dec tau))) + (* (gamma alpha) + (gamma tau) + (pow (+ lambda x) (+ alpha tau)))))) + +(defmethod cdf ::GeneralizedPareto [p x] + {:pre (in-domain p x)} + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda]) + tau (get-in p [:Parameters :tau])] + (cdf (beta-dist tau alpha) (/ x (+ lambda x))))) + +(defmethod rawmoments ::GeneralizedPareto [p k] + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda]) + tau (get-in p [:Parameters :tau])] + (* (pow lambda k) + (/ (reduce * (map #(+ tau %) (range k))) + (reduce * (map #(- alpha %) (range 1 (inc k)))))))) + +(defmethod truncated-mean ::GeneralizedPareto [p d] + (let [alpha (get-in p [:Parameters :alpha]) + lambda (get-in p [:Parameters :lambda]) + tau (get-in p [:Parameters :tau])] + (* (mean p) + (cdf (beta-dist (inc tau) (dec alpha)) + (/ d (+ lambda d)))))) + +;; ## Burr + +(defmethod pdf ::Burr [b x] + {:pre (in-domain b x)} + (let [alpha (get-in b [:Parameters :alpha]) + lambda (get-in b [:Parameters :lambda]) + tau (get-in b [:Parameters :tau])] + (/ (* alpha tau (pow lambda alpha) (pow x (dec tau))) + (pow (+ lambda (pow x tau)) + (inc alpha))))) + +(defmethod cdf ::Burr [b x] + {:pre (in-domain b x)} + (let [alpha (get-in b [:Parameters :alpha]) + lambda (get-in b [:Parameters :lambda]) + tau (get-in b [:Parameters :tau])] + (- 1 (pow (/ lambda + (+ lambda + (pow x tau))) + alpha)))) + +(defmethod rawmoments ::Burr [b k] + {:pre [(< (-(get-in b [:Parameters :tau])) k) + (< k (* (get-in b [:Parameters :alpha]) + (get-in b [:Parameters :tau])))]} + (let [alpha (get-in b [:Parameters :alpha]) + lambda (get-in b [:Parameters :lambda]) + tau (get-in b [:Parameters :tau])] + (/ (* (pow lambda (/ k tau)) + (gamma (inc (/ k tau))) + (gamma (- alpha (/ k tau)))) + (gamma alpha)))) + +(defmethod truncated-mean ::Burr [b d] + (let [alpha (get-in b [:Parameters :alpha]) + lambda (get-in b [:Parameters :lambda]) + tau (get-in b [:Parameters :tau])] + (* (mean b) + (cdf (beta-dist (inc (/ 1 tau))(- alpha (/ 1 tau))) + (/ (pow d tau)(+ lambda (pow d tau))))))) + +(defmethod value-at-risk ::Burr [b kappa] + {:pre [(<= 0 kappa 1)]} + (let [alpha (get-in b [:Parameters :alpha]) + lambda (get-in b [:Parameters :lambda]) + tau (get-in b [:Parameters :tau])] + (pow (* lambda + (dec (pow (- 1 kappa) + (- (/ 1 alpha))))) + (/ 1 tau)))) + +;; ## Log-logistique +(defmethod pdf ::LogLogistic [l x] + {:pre (in-domain l x)} + (let [lambda (get-in l [:Parameters :lambda]) + tau (get-in l [:Parameters :tau])] + (/ (* tau + (pow x (dec tau))) + (pow (+ (pow lambda tau) + (pow x tau)) + 2)))) + +(defmethod cdf ::LogLogistic [l x] + {:pre (in-domain l x)} + (let [px (pow x (get-in l [:Parameters :tau])) + pl (pow (get-in l [:Parameters :lambda]) + (get-in l [:Parameters :tau]))] + (/ px (+ pl px)))) + +(defmethod rawmoments ::LogLogistic [l k] + {:pre (< (abs k) (get-in l [:Parameters :tau]))} + (let [lambda (get-in l [:Parameters :lambda]) + tau (get-in l [:Parameters :tau])] + (* (pow lambda k) + (gamma (inc (/ k tau))) + (gamma (- 1 (/ k tau)))))) + +(defmethod truncated-mean ::LogLogistic [l d] + {:pre (> (get-in l [:Parameters :tau]) 1)} + (let [lambda (get-in l [:Parameters :lambda]) + tau (get-in l [:Parameters :tau])] + (* (mean l) + (cdf (beta-dist (inc (/ 1 tau))(- 1 (/ 1 tau))) + (/ (pow d tau)(+ lambda (pow d tau))))))) + +(defmethod value-at-risk ::LogLogistic [l kappa] + {:pre [(<= 0 kappa 1)]} + (let [lambda (get-in l [:Parameters :lambda]) + tau (get-in l [:Parameters :tau])] + (* lambda (pow (dec (/ 1 kappa)) + (- (/ 1 tau)))))) + +;; ## Student + +(defmethod pdf ::Student [s x] + {:pre (in-domain s x)} + (let [nu (get-in s [:Parameters :nu])] + (case nu + 1 (/ 1 (* Math/PI (+ 1 (pow x 2)))) + 2 (/ 1 (pow (+ 2 (pow x 2)) (/ 3 2))) + (/ (* (gamma (/ (+ 1 nu) 2)) + (pow (+ 1 (/ (pow x 2) nu)) (- (/ (+ nu 1) 2)))) + (* (sqrt (* nu Math/PI)) + (gamma (/ nu 2))))))) + +(defmethod cdf ::Student [s x] + {:pre (in-domain s x)} + (let [nu (get-in s [:Parameters :nu]) + corepart (/ (cdf (beta-dist (/ nu 2) 0.5) + (/ nu (+ nu (pow x 2)))) 2)] + (case nu + 1 (+ (/ 1 2) (/ (Math/atan x) Math/PI)) + 2 (/ (+ 1 (/ x (sqrt (+ 2 (pow x 2))))) 2) + (cond + (< x 0.0) corepart + (= x 0.0) 0.5 + (> x 0.0) (- 1 corepart) + )))) + +(defmethod rawmoments ::Student [s k] + {:pre [(< 0 k (get-in s [:Parameters :nu]))]} + (let [nu (get-in s [:Parameters :nu])] + (cond + (odd? k) 0 + (even? k) (/ (* (gamma (/ (inc k) 2)) + (gamma (/ (- nu k) 2)) + (pow nu (/ k 2))) + (* (sqrt Math/PI) (gamma (/ nu 2)))) + ))) + +(defmethod truncated-mean ::Student [s d] + (let [nu (get-in s [:Parameters :nu]) + main-part (* (sqrt (/ nu Math/PI)) + (/ (gamma (/ (inc nu) 2)) (gamma (/ nu 2))) + (pow (inc (/ (pow d 2) nu)) + (/ (- 1 nu) 2))) + ] + (cond + (< d 0) main-part + (> d 0) (- main-part))) + ) + + diff --git a/src/actuarial/http.clj b/src/actuarial/http.clj new file mode 100644 index 0000000..78e6e8d --- /dev/null +++ b/src/actuarial/http.clj @@ -0,0 +1,50 @@ +(ns actuarial.http + (:use ring.adapter.jetty) + (:use ring.middleware.params) + (:use actuarial.distributions)) + +(require 'ring.middleware.stacktrace) + +;(describe (build-distribution (resolve (symbol (str "gaussian" "-dist"))) [1 3]) 0.5 0.25 3) + +(defn wrap-spy [handler] + (fn [request] + (println "-------------------------------") + (println "Incoming Request:") + (clojure.pprint/pprint request) + (let [response (handler request)] + (println "Outgoing Response Map:") + (clojure.pprint/pprint response) + (println "-------------------------------") + response))) + +(defmacro get-map-dist [map name param] + `(apply (get ~map ~name ) ~param)) + +(defn handler [req] + (let [distr-name (get-in req [:query-params "dist"]) + distr-params (read-string (get-in req [:query-params "params"])) + distr-x (read-string (get-in req [:query-params "x"] "0")) + distr-d (read-string (get-in req [:query-params "d"] "0")) + distr-t (read-string (get-in req [:query-params "t"] "1")) + distr-map (get-map-dist distributions-map distr-name distr-params) + sortie (describe distr-map distr-x distr-d distr-t) + ] + (println distr-name distr-params) + {:status 200 + :headers {"Content-Type" "text/html"} + :body (str "

" sortie "

")} + )) + +(def app + (-> #'handler + (ring.middleware.stacktrace/wrap-stacktrace) + (wrap-spy) + (wrap-params) + )) + +(defonce server (ring.adapter.jetty/run-jetty #'app {:port 8080 :join? false})) + +(.stop server) + +(.start server) \ No newline at end of file diff --git a/src/actuarial/measures.clj b/src/actuarial/measures.clj new file mode 100644 index 0000000..60dd623 --- /dev/null +++ b/src/actuarial/measures.clj @@ -0,0 +1,6 @@ +(ns actuarial.measures (:use actuarial.basic)) +(use '(incanter core)) + +(defn node-entropy [& values] + (- (reduce + (map #(* % (/ (log %) (log 2))) + (map #(/ % (reduce + values)) values))))) diff --git a/src/actuarial/mixed_distributions.clj b/src/actuarial/mixed_distributions.clj new file mode 100644 index 0000000..2244e80 --- /dev/null +++ b/src/actuarial/mixed_distributions.clj @@ -0,0 +1,32 @@ +(ns actuarial.mixed-distributions (:use actuarial.distributions)) + +(defn mixedDist + [prob-dist cont-dist] + {:Distribution ::Mixed + :SubDistributions {:prob-dist prob-dist + :cont-dist cont-dist} + :type ::Mixed}) + +(defmacro getmixed [f m & params] + `(let [probs# (get-in ~m [:SubDistributions :prob-dist :Parameters :prob]) + dists# (get-in ~m [:SubDistributions :cont-dist])] + (reduce + (map #(* %1 (~f %2 ~@params)) probs# dists#)))) + +(defmethod pdf ::Mixed [m x] + (getmixed pdf m x)) + +(defmethod cdf ::Mixed [m x] + (getmixed cdf m x)) + +(defmethod mgf ::Mixed [m t] + (getmixed mgf m t)) + +(defmethod rawmoments ::Mixed [m k] + (getmixed rawmoments m k)) + +(defmethod truncated-mean ::Mixed [m d] + (getmixed truncated-mean m d)) + + + + diff --git a/test/actuarial/test/core.clj b/test/actuarial/test/core.clj new file mode 100644 index 0000000..cee153a --- /dev/null +++ b/test/actuarial/test/core.clj @@ -0,0 +1,123 @@ +(ns actuarial.core-test + (:require [clojure.test :refer :all] + [actuarial.core :refer :all] + [actuarial.distributions :as dist] + [actuarial.mixed-distributions :as mix])) +;; # Discrete distributions + +;; ## Arithmetic distribution + +(defmacro dbg [& body] + `(let [x# ~@body] + (println (str "dbg: " (quote ~@body) "=" x#)) + x#)) + +(defn close? [tolerance x y] + (< (- x y) tolerance)) + +(defn =? [x y] (partial close? 0.000000001)) + +(deftest arith-dist + (let [arith1 (dist/arithmetic-dist 1 0.2 0.3 0.4 0.1)] + (testing "Arithmetic Distribution" + (is (=? (dist/pmf arith1 3) 0.1)) + (is (=? (dist/cdf arith1 2) 0.9)) + (is (=? (dist/survival arith1 2) 0.1)) + (is (=? (dist/mean arith1) 1.4)) + (is (=? (dist/variance arith1) 1.4)) + (is (=? (dist/value-at-risk arith1 0.6) 2)) + (is (=? (dist/tail-value-at-risk arith1 0.6) 2.25)) + ))) + +;; ## Poisson distribution + +(deftest poisson-dist + (let [poisson1 (dist/poisson-dist 3) + =? (partial close? 0.00001)] + (testing "Poisson distribution" + (is (=? (dist/pmf poisson1 0) (Math/exp -3))) + (is (=? (dist/mean poisson1) 3)) + (is (=? (dist/variance poisson1) 3))+ + (is (=? (dist/pgf poisson1 2) (Math/exp 3))) + (is (=? (dist/mgf poisson1 2) (Math/exp(* 3 (- (Math/exp 2) 1))))) + ))) + +;; ## Mixed Distributions + +;; tests + +(def dist1 (dist/arithmetic-dist 1 0.4 0.6)) +(def dist2 (map #(apply dist/exponential-dist %) '([0.4] [0.5]))) +(def mix1 (mix/mixedDist dist1 dist2)) +(println + (dist/pdf mix1 0.25) + (dist/cdf mix1 0.25) + (dist/mgf mix1 0.1) + (dist/rawmoments mix1 3)) + +(dist/describe mix1 0.25 0.5 2 0.99) + +;; # Continuous distributions + +;; Uniform distribution +(def uniform1 (dist/uniform-dist 1.0 3.0)) +(println (dist/describe uniform1 2.0 2.5 2.0 [0.95])) + +;; Exponential distribution +(def exp1 (dist/exponential-dist 2.0)) +(println (dist/describe exp1 2.0 2.5 1.0 [0.95])) + +;; Beta distribution +(def beta1 (dist/beta-dist 1 3)) +(println (dist/describe beta1 0.75 0.25 2 [0.95])) + +;; Gamma distribution +(def gamma1 (dist/gamma-dist 2.001 4.999)) +(println (dist/describe gamma1 0.75 0.25 2 [0.95])) + +;; Erlang distribution +(def erlang1 (dist/erlang-dist 2 5)) +(println (dist/describe erlang1 0.75 0.25 2 [0.95])) + +;; Weibull distribution +(def weibull1 (dist/weibull-dist 2.2 3.1)) +(println (dist/describe weibull1 0.75 0.25 2 [0.95])) + +;; Generalized Erlang distribution +(def generlang1 (dist/generalized-erlang-dist 2.1 5 0.25)) +(println (dist/describe generlang1 0.75 0.25 2 [0.95])) + +;; Gaussian distribution +(def gaussian1 (dist/gaussian-dist 1.5 0.5)) +(println (dist/describe gaussian1 0.75 0.25 2 [0.95])) + +;;Lognormal distribution +(def lognormal1 (dist/lognormal-dist 1.5 0.5)) +(println (dist/describe lognormal1 0.75 0.25 2 [0.95])) + +;;Inverse Gaussian distribution +(def inversegaussian1 (dist/inverse-gaussian-dist 1.5 0.5)) +(println (dist/describe inversegaussian1 0.75 0.25 2 [0.95])) + +;;Pareto distribution +(def pareto1 (dist/pareto-dist 4 3)) +(println (dist/describe pareto1 0.75 0.25 2 [0.95])) + +;;Generalized Pareto distribution +(def generalizedpareto1 (dist/generalized-pareto-dist 4 3 1)) +(println (dist/describe generalizedpareto1 0.75 0.25 2 [0.95])) + +;;Burr distribution +(def burr1 (dist/burr-dist 4 3 1)) +(println (dist/describe burr1 0.75 0.25 2 [0.95])) + +;;Log-Logistic distribution +(def loglog1 (dist/log-logistic-dist 3 2.1)) +(println (dist/describe loglog1 0.75 0.25 2 [0.95])) + +;;Student distribution +(def student1 (dist/student-dist 3)) +(println (dist/describe student1 0.75 0.25 0.5 [0 1])) + +(run-tests) +