            D (diagonal-matrix (map (fn [x] (sqrt (max x 0))) S))]
        (add mean-vector (mmul U D z))))
     ([d n] (take n (repeatedly #(sample d)))))
+  GeneralizedDoubleParetoDistribution
+  (sample
+    ([d]
+     (let [shape (:shape d)
+           scale (:scale d)
+           eta (* shape scale)
+           lambda (sample (gamma shape eta))
+           tau (sample (exponential (/ (square lambda) 2)))]
+       (sample (normal 0 tau))))
+    ([d n] (take n (repeatedly #(sample d)))))
   )
 
 (defn expectation

        p-double (double-array p)]
    (new EnumeratedIntegerDistribution x-int p-double)))
(defn exponential
  [mean]
  (new ExponentialDistribution mean))
(defn f-distribution
  [df1 df2]
  (new FDistribution df1 df2))
(defn gamma
  [shape rate]
  (new GammaDistribution shape (/ 1 rate)))
(defn geometric
  [p]
  (new GeometricDistribution p))

  RealDistribution
  (variance [d] (.getNumericalVariance d))
  )

(defn inv [x] (/ 1 x))

(defprotocol proximal
  (prox [d] [d h x]))
(extend-protocol proximal
  NormalDistribution
  (prox
    ([d] (fn [h x] (* x (/ (inv h) (+ (inv h) (inv (variance d)))))))
    ([d h x] ((prox d) h x)))
  LaplaceDistribution
  (prox
    ([d] (fn [h x] (* (signum x) (max 0.0 (- (abs x) (* h (rate d)))))))
    ([d h x] ((prox d) h x)))
  GeneralizedDoubleParetoDistribution
  (prox
    ([d]
     (let [shape (:shape d)
           scale (:scale d)]
       (fn [h x]
         (let [g (* h (inc shape))
               a (* shape scale)
               d (max 0 (- (* a (abs x)) g))]
           (* 0.5 (signum x) (+ (abs x) (negate a) (sqrt (+ (* 4 d) (square (- a (abs x)))))))))))
    ([d h x]
     ((prox d) h x))))

(defprotocol random
  (sample [d] [d n]))
