(defn truncated-variance
" Returns the variance of a normal distribution truncated at a and b.

  Options:
    :mean (default 0) mean of untruncated normal distribution
    :sd (default 1) standard deviation of untruncated normal distribution
    :a (default -infinity) lower truncation point
    :b (default +infinity) upper truncation point

  Examples:

    (use '(incanter core stats))
    (truncated-variance :a -1.96 :b 1.96)
    (truncated-variance :a 0)
    (truncated-variance :b 0)

    (use 'incanter.charts)
    (def x (range -3 3 0.1))
    (def plot (xy-plot x (map #(truncated-variance :a %) x)))
    (view plot)
    (add-lines plot x (map #(truncated-variance :b %) x))

    (def samp (sample-normal 10000))
    (add-lines plot x (map #(variance (filter (fn [s] (> s %)) samp)) x))
    (add-lines plot x (map #(variance (mult samp (indicator (fn [s] (> s %)) samp))) x))


  References:
    DeMaris, A. (2004) Regression with social data: modeling continuous and limited response variables.
      Wiley-IEEE.

    http://en.wikipedia.org/wiki/Truncated_normal_distribution
"
  ([& options]
    (let [opts (when options (apply assoc {} options))
          mu (or (:mean opts) 0)
          sd (or (:sd opts) 1)
          sigma-sq (* sd sd)
          a (or (:a opts) Double/NEGATIVE_INFINITY)
          b (or (:b opts) Double/POSITIVE_INFINITY)
          lambda (fn [alpha] (/ (pdf-normal alpha) (- 1 (cdf-normal alpha))))
          delta (fn [alpha] (* (lambda alpha) (- (lambda alpha) alpha)))
          ;one-tail-var (fn [alpha s-sq]
          ;               (* s-sq
          ;                  (- 1 (cdf-normal alpha))
          ;                  (+ (- 1 (delta alpha))
          ;                     (* (pow (- alpha (lambda alpha)) 2)
          ;                        (cdf-normal alpha)))))
          a-std (if (= a Double/NEGATIVE_INFINITY) Double/NEGATIVE_INFINITY (/ (- a mu) sd))
          b-std (if (= b Double/POSITIVE_INFINITY) Double/POSITIVE_INFINITY (/ (- b mu) sd))
          pdf-a (if (= a Double/NEGATIVE_INFINITY) 0 (pdf-normal a-std))
          pdf-b (if (= b Double/POSITIVE_INFINITY) 0 (pdf-normal b-std))
          cdf-a (if (= a Double/NEGATIVE_INFINITY) 0 (cdf-normal a-std))
          cdf-b (if (= b Double/POSITIVE_INFINITY) 1 (cdf-normal b-std))]
      (cond
        (and (= b Double/POSITIVE_INFINITY) (= a Double/NEGATIVE_INFINITY))
          sigma-sq
        (and (= b Double/POSITIVE_INFINITY) (> a Double/NEGATIVE_INFINITY))
          (* sigma-sq (- 1 (delta a-std)))
          ;(one-tail-var a-std sigma-sq)
        (and (= a Double/NEGATIVE_INFINITY) (< b Double/POSITIVE_INFINITY))
          (* sigma-sq (- 1 (delta (- 1 b-std))))
          ;(- sigma-sq (one-tail-var b-std sigma-sq))
        :else
          (* sigma-sq
            (+ 1 (/ (- (* a-std pdf-a) (* b-std pdf-b))
                    (- cdf-b cdf-a))
              (- (pow (/ (- pdf-a pdf-b) (- cdf-b cdf-a)) 2))))))))