(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)))
  )