Tuesday, January 26, 2010

1987 Evett Fiber Model

(defn numerator [y x]
(let [m (nrow x)
fnc1 (fn [t u]
(- (/ (* (Math/exp (- u)) (- t 1) (Math/pow u (- t 2)))
(- t 2))
(* (- (Math/exp (- u))) (Math/pow u (- t 1)))))
fnc2 (fn []
(/ (* (- m 1) (+ m 1)) m))
smean (fn [x]
(div (reduce plus x) m))
Sx (fn [x]
(let [vals (div
(reduce plus
(matrix (map #(pow (minus % (smean x)) 2) x))) m)
comat (matrix [[(first vals) 0] [0 (second vals)]])
icomat (matrix [[(second vals) 0] [0 (first vals)]])]
[comat icomat]))]

(/ (/ (fnc1 (/ m 2) (first y)) (* (Math/PI) (fnc1 (/ (- m 2) 2) (first y))))
(mult (pow (det (mult (fnc2) (first (Sx x)))) 0.5)
(Math/pow (+ 1 (mmult (mmult (mult (minus y (smean x)) (fnc2))
(second (Sx x)))
(trans (minus y (smean x)))))
(/ m 2))))))

(defn denominator [y z lamb]
(let [Normal-zlamb (fn [zi]
(* (/ 1 (Math/sqrt (* 2 (Math/PI) (Math/pow lamb 2))))
(Math/exp (- (/ (mmult (minus y zi) (trans (minus y zi)))
(* 2 (Math/pow lamb 2)))))))
SumZ (/ (apply + (map #(Normal-zlamb %) z)) (count z))]
SumZ))

(defn likelihood [x y z lamb]
(/ (numerator y x) (denominator y z lamb)))

No comments:

Post a Comment