Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also compare across forks.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also compare across forks.
base: 598e68e9cd
...
compare: 995232be9e
  • 3 commits
  • 10 files changed
  • 0 commit comments
  • 1 contributor
View
4 src/continuous-time.lisp
@@ -19,7 +19,7 @@ two underlying random variables.")
(keys (when keys (coerce keys 'vector)))
(n (length rates))
(total-rate (sum rates))
- ((&values transition-rate probabilities)
+ ((&values transition-rate probabilities)
(if transition-rate
(let ((no-change-rate (- transition-rate total-rate)))
(assert (<= 0 no-change-rate) ()
@@ -39,7 +39,7 @@ two underlying random variables.")
(draw (&key)
(values (draw duration)
(let ((j (draw jump)))
- (cond
+ (cond
((= j n) no-change)
(keys (aref keys j))
(t j))))))
View
10 src/design-matrix.lisp
@@ -5,11 +5,11 @@
;;; building blocks for a DSL for design matrices
;;; (list term1 ...)
-;;;
+;;;
;;; term := constant | covariate | interaction
-;;;
+;;;
;;; covariate := symbol | (^ symbol exponent)
-;;;
+;;;
;;; interaction := (* covariate1 covariate2 ...)
(defun interaction-matrix (&rest matrices)
@@ -99,7 +99,7 @@ symbol size)."
(collecting '#:* :into interaction-names))
(if (atom name)
(collecting name :into interaction-names)
- (progn
+ (progn
(collecting (first name) :into interaction-names)
(multiplying (second name) :into size)))
(finally
@@ -118,7 +118,7 @@ PROCESS-FACTOR. When CONSTANT is non-nil, a constant column with this name will
be added. Return IX specification, the matrix, and a list of factors and levels
as values.
-Example:
+Example:
(design-matrix (clo :integer
1 2 3 :/
4 5 6)
View
2  src/log-infinity.lisp
@@ -14,7 +14,7 @@
(defun ~- (a &rest args)
"Subtraction, treating NIL as -Infinity."
(if args
- (let ((args-sum
+ (let ((args-sum
(iter
(for arg :in args)
(unless arg
View
20 src/multivariate.lisp
@@ -7,7 +7,7 @@
(defun check-mean-variance-compatibility (mean variance)
"Assert that the mean is a vector, and its dimensions are compatible with
variance."
- (assert (and (vectorp mean)
+ (assert (and (vectorp mean)
(= (length mean) (nrow variance) (ncol variance)))))
(defun normal-quadratic-form (x mean variance-left-sqrt)
@@ -20,7 +20,7 @@ square root of variance."
:instance instance)
((n :type fixnum :documentation "Number of dimensions.")
(mean :type vector :reader t)
- (variance-left-sqrt :reader t :documentation
+ (variance-left-sqrt :reader t :documentation
"Left square root of the variance matrix."))
(let ((variance-left-sqrt (left-square-root variance)))
(check-mean-variance-compatibility mean variance-left-sqrt)
@@ -53,8 +53,8 @@ square root of variance."
;;; inverse-chi-square, nu degrees of freedom) is returned as the second
;;; value.
-(define-rv r-multivariate-t (mean sigma nu
- &key multivariate-normal scaling-factor
+(define-rv r-multivariate-t (mean sigma nu
+ &key multivariate-normal scaling-factor
(s^2 1d0 s^2?))
(:documentation "Multivariate T distribution with given MEAN, scale
SIGMA*S^2 and NU degrees of freedom.")
@@ -63,7 +63,7 @@ square root of variance."
(scaling-factor :type r-inverse-gamma :documentation
"distribution that scales the variance of draws."
:reader t))
- (make :multivariate-normal
+ (make :multivariate-normal
(aif multivariate-normal
(prog1 it
(check-type it r-multivariate-normal)
@@ -80,7 +80,7 @@ square root of variance."
"Can't initialize both S^2 and SCALING-FACTOR.")
it)
(r-inverse-chi-square nu s^2)))
- (mean ()
+ (mean ()
(assert (< 1 (nu scaling-factor)))
(mean multivariate-normal))
(variance ()
@@ -112,9 +112,9 @@ square root of variance."
(normal (sub multivariate-normal index-specification))
((&accessors-r/o nu s^2) scaling-factor))
(etypecase normal
- (r-normal
+ (r-normal
(r-t (mean normal) (* s^2 (variance normal)) nu))
- (r-multivariate-normal
+ (r-multivariate-normal
(r-multivariate-t nil nil nil
:multivariate-normal normal
:scaling-factor scaling-factor)))))
@@ -154,12 +154,12 @@ distribution (dimension k x k)."
(check-type nu (and fixnum (satisfies plusp)))
(make :nu nu :scale-left-sqrt scale-left-sqrt :k (nrow scale-left-sqrt)))
(mean () (e* nu (mm scale-left-sqrt t)))
- (draw (&key)
+ (draw (&key)
(mm (mm scale-left-sqrt (draw-standard-wishart-left-sqrt nu k)) t)))
;;; Inverse-Wishart
;;;
-;;; If A ~ Inverse-Wishart[nu,inverse-scale], then
+;;; If A ~ Inverse-Wishart[nu,inverse-scale], then
;;; (invert A) ~ Wishart(nu,inverse-scale).
(define-rv r-inverse-wishart (nu inverse-scale)
View
51 src/optimization.lisp
@@ -67,7 +67,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
(setf (aref y j) (aref x j))))
(lla::make-matrix% n n ddfx :kind :hermitian)))
-(defmacro loop-max-iter ((max-iter
+(defmacro loop-max-iter ((max-iter
&optional (end-form '(error 'reached-max-iter)))
&body body)
"Execute loop count down using max-iter, call end-form when it reaches zero."
@@ -101,7 +101,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
(tau2 :initarg :tau2 :initform 0.1d0 :documentation "0 < tau2 < tau3 <= 0.5
keep the bracketing algorithm away from the endpoints of the interval. tau2
<= sigma is suggested, typical values are tau2=0.1 and tau3=0.5.")
- (tau3 :initarg :tau3 :initform 0.5d0
+ (tau3 :initarg :tau3 :initform 0.5d0
:documentation "see tau2 of the same class.")
;; parameters for Armijo line search algorithm
(step-reduction :initarg :step-reduction :initform 0.2d0
@@ -145,7 +145,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; ;;; Algorithm for strong Wolfe conditions.
;; (defun sw-find-acceptable-alpha (f df a fa dfa b fb dfb
-;; f0 slope curvature max-iter epsilon
+;; f0 slope curvature max-iter epsilon
;; bfgs-parameters)
;; ;; notes: dfb may be NIL, in which case quadratic interpolation is used
;; (bind (((:slots-r/o tau2 tau3) bfgs-parameters))
@@ -186,7 +186,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; (>= f-alpha f-alpha-prev))
;; (return
;; (sw-find-acceptable-alpha
-;; f df
+;; f df
;; alpha-prev f-alpha-prev df-alpha-prev
;; alpha f-alpha nil
;; f0 slope curvature linesearch-max-iter epsilon bfgs-parameters)))
@@ -194,21 +194,21 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; (when (<= (abs df-alpha) curvature)
;; (return (values alpha f-alpha)))
;; (when (<= 0 df-alpha)
-;; (return
+;; (return
;; (sw-find-acceptable-alpha
-;; f df
+;; f df
;; alpha f-alpha df-alpha
-;; alpha-prev f-alpha-prev df-alpha-prev
+;; alpha-prev f-alpha-prev df-alpha-prev
;; f0 slope curvature linesearch-max-iter epsilon bfgs-parameters)))
;; (let* ((left (- (* 2 alpha) alpha-prev))
-;; (alpha-next
+;; (alpha-next
;; (if (<= alpha-max left)
;; alpha-max
;; (cubic-minimum alpha f-alpha df-alpha
;; alpha-prev f-alpha-prev df-alpha-prev
;; left
;; (min alpha-max
-;; (+ alpha
+;; (+ alpha
;; (* tau1 (- alpha alpha-prev))))))))
;; (setf alpha-prev alpha
;; f-alpha-prev f-alpha
@@ -217,7 +217,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
-;; (defun bfgs-minimize (f x &key (df nil df?)
+;; (defun bfgs-minimize (f x &key (df nil df?)
;; (delta (expt double-float-epsilon 1/2))
;; (epsilon (expt double-float-epsilon 1/2))
;; gamma
@@ -228,7 +228,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; count-f-eval?
;; (bfgs-parameters *default-bfgs-parameters*))
;; "Parameters:
-
+
;; f : An R^n=>R function to be optimized.
;; x : Initial point.
@@ -242,7 +242,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; epsilon : Convergence criterion in f(x):
;; |change in f(x)| <= epsilon*(1+|f(x)|). Ignored when nil.
-;; gamma : Convergence criterion in df(x).
+;; gamma : Convergence criterion in df(x).
;; |change in df(x)|_sup <= gamma. Ignored when nil.
;; Convergence happens when all three are met (except for those which are
@@ -258,7 +258,7 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; difference will be used. Note that the latter is quite cheap, since it
;; requires only two evaluations.
-
+
;; Return
;; (values
@@ -301,14 +301,14 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; (df-uni (if (and df? use-df-for-linesearch?)
;; (lambda (alpha) (dot (funcall df (e+ x (e* p alpha))) p))
;; (numdiff1* f-uni numdiff-epsilon)))
-;; ((:values alpha f-alpha)
-;; (funcall linesearch f-uni df-uni fx (dot dfx p) 1d0
+;; ((:values alpha f-alpha)
+;; (funcall linesearch f-uni df-uni fx (dot dfx p) 1d0
;; epsilon bfgs-parameters))
;; (s (e* p alpha)))
;; ;; check convergence (using relative changes)
;; (when (and (below? (/ (abs (- fx f-alpha)) (1+ (abs fx))) epsilon)
;; (below? (/ (norm2 s) (1+ (norm2 x))) delta))
-;; (d:d "rel conv: ~a, ~a"
+;; (d:d "rel conv: ~a, ~a"
;; (/ (abs (- fx f-alpha)) (1+ (abs fx)))
;; (/ (norm2 s) (1+ (norm2 x))))
;; (done iter-count))
@@ -317,8 +317,8 @@ returned. Uses Richardson extrapolation w/ forward differencing."
;; (setf x (e+ x s)
;; fx f-alpha
;; dfx (funcall df x)
-;; H (bfgs-update-inverse-hessian H (e- dfx dfx-prev) s
-;; (and (not H?)
+;; H (bfgs-update-inverse-hessian H (e- dfx dfx-prev) s
+;; (and (not H?)
;; (first-iteration-p)))))))))
(defun negative-quadratic-form (x A negative-Ax)
@@ -350,7 +350,7 @@ All vectors are assumed to have element type DOUBLE-FLOAT."
(setf (aref negative-Ax i) negative-sum)))
negative-xAX))
-(defun x+direction (x direction alpha &optional
+(defun x+direction (x direction alpha &optional
(z (lla-array (length x) :double)) (relative 16d0))
"Calculate X+ALPHA*DIRECTION, place it in Z. All vectors (X, DIRECTION,
Z) are assumed to have element type double-float (not checked), and the same
@@ -380,11 +380,11 @@ iff all are =. Z is returned as the second value."
;;; alpha - starting alpha
;;; z - resulting x ends up here
;;; bfgs-parameters - all other parameters picked from here
-;;;
+;;;
;;; Return values
;;; alpha - the suggested alpha
;;; f-alpha - the value at this alpha
-;;; convergence indicator: NIL if OK, :ZERO if converged to X,
+;;; convergence indicator: NIL if OK, :ZERO if converged to X,
(defun linesearch-armijo (f df-uni x direction fx df0 alpha z bfgs-parameters)
"Armijo line search. Last x evaluated is available in z."
@@ -428,7 +428,7 @@ and ACCEPTED?"
((< (funcall df-uni alpha f-alpha) curvature) (setf left alpha))
;; both hold
(t (done nil)))
- (setf alpha
+ (setf alpha
(if right
(if (<= (incf bisections) max-bisections)
(/ (+ left right) 2d0)
@@ -474,7 +474,7 @@ too small), in this case, no argument is changed.."
t))
-(defun bfgs-minimize (f x &key (df nil df?)
+(defun bfgs-minimize (f x &key (df nil df?)
; (delta (expt double-float-epsilon 1/2))
(epsilon (expt double-float-epsilon 1/2))
;; gamma
@@ -547,7 +547,7 @@ too small), in this case, no argument is changed.."
(funcall f y)))
0d0)))
((&values alpha f-alpha status)
- (funcall linesearch #'f #'df-uni x s fx df0 1d0 z
+ (funcall linesearch #'f #'df-uni x s fx df0 1d0 z
bfgs-parameters)))
(when (<= (/ (abs (- fx f-alpha)) (+ (abs fx) epsilon)) epsilon)
(done))
@@ -556,7 +556,7 @@ too small), in this case, no argument is changed.."
;; s is scaled to be the step
(dotimes (i n)
(multf (aref s i) alpha))
- ;; calculate new gradient and difference,
+ ;; calculate new gradient and difference,
;; overwrite old gradient
(setf y (funcall df z))
(rotatef y dfx) ; y: old, dfx: new
@@ -576,4 +576,3 @@ too small), in this case, no argument is changed.."
(reset "linesearch converged to current point")))
;; nonnegative gradient, try to reset
(reset "nonnegative gradient"))))))
-
View
4 src/random.lisp
@@ -93,7 +93,9 @@ Also, within BODY, slots are accessible by their names."
(declaim (inline check-probability))
(defun check-probability (p &optional limits)
- "Assert that P is a probability (ie a real number between 0 and 1)."
+ "Assert that P is a probability (ie a real number between 0 and 1). When
+LIMITS is given, it is checked that p is not 0 (:LEFT), 1 (:RIGHT), or
+0/1 (:BOTH)."
(assert (<= 0 p 1) () "~A is not a valid probability." p)
(when limits
(let ((msg "The given probability is only attained in the limit."))
View
10 src/special-functions.lisp
@@ -7,10 +7,10 @@
calculating the Gamma function using the Lanczos method."
;; coefficients are from http://lib.stat.cmu.edu/apstat/245
;; !!! compare with http://home.att.net/~numericana/answer/info/godfrey.htm
- (let ((coefficients #(0.9999999999995183d0 676.5203681218835d0
- -1259.139216722289d0 771.3234287757674d0
- -176.6150291498386d0 12.50734324009056d0
- -0.1385710331296526d0 0.9934937113930748d-05
+ (let ((coefficients #(0.9999999999995183d0 676.5203681218835d0
+ -1259.139216722289d0 771.3234287757674d0
+ -176.6150291498386d0 12.50734324009056d0
+ -0.1385710331296526d0 0.9934937113930748d-05
0.1659470187408462d-06)))
`(+ ,(aref coefficients 0) ,@(iter
(for c :in-vector coefficients :from 1)
@@ -53,7 +53,7 @@ conditions can be ensured."
(t
(error 'division-by-zero)))))
-(defun gamma% (z)
+(defun gamma% (z)
"Internal function implementing for calculating Gamma(z) for z > 0,
z double-float. Not safe to be called directly, unless these
conditions can be ensured."
View
35 src/univariate.lisp
@@ -66,7 +66,7 @@ which has density exp(-x)."
(with-doubles (p)
(check-probability p :right)
(/ (log (- 1 p)) (- beta))))
- (draw (&key)
+ (draw (&key)
(/ (draw-standard-exponential) beta)))
@@ -90,7 +90,7 @@ which has density exp(-x)."
(y (+ (abs v) 0.386595d0))
(q (+ (expt x 2) (* y (- (* 0.19600d0 y) (* 0.25472d0 x))))))
(if (and (> q 0.27597d0)
- (or (> q 0.27846d0)
+ (or (> q 0.27846d0)
(plusp (+ (expt v 2) (* 4 (expt u 2) (log u))))))
(go top)
(return-from draw-standard-normal (/ v u))))))
@@ -102,7 +102,7 @@ which has density exp(-x)."
(/ (- x mu) sigma))
(defun from-standard-normal (x mu sigma)
- "Scale x from standard normal."
+ "Scale x from standard normal."
(+ (* x sigma) mu))
(defun cdf-normal% (x mean sd)
@@ -198,7 +198,11 @@ respectively), on the interval [left, \infinity).")
(cdf (x) (if (<= left x)
(/ (1- (+ (cdf-normal% x mu sigma) m0)) m0)
0d0))
- (mean () (truncated-normal-moments% 1 mu sigma left nil))
+ (quantile (q)
+ (with-doubles (q)
+ (check-probability q :right)
+ (rmath:qnorm5 (+ (* q m0) (1c m0)) mu sigma 1 0)))
+ (mean () (truncated-normal-moments% 1 mu sigma left nil))
(variance () (truncated-normal-moments% 2 mu sigma left nil)))
(defun r-truncated-normal (left right &optional (mu 0d0) (sigma 1d0))
@@ -333,7 +337,7 @@ respectively), on the interval [left, \infinity).")
;; left alpha)))
;; (<= x right) x))
;; ;; optimal to use the uniform-based reject/accept
-;; (lambda* (draw-left-right-truncated-standard-normal
+;; (lambda* (draw-left-right-truncated-standard-normal
;; left width (* (expt left 2) 0.5d0))))))
;; ;; whole support below 0, will flip
;; (t
@@ -347,7 +351,7 @@ respectively), on the interval [left, \infinity).")
;; left alpha)))
;; (<= x right) x))
;; ;; optimal to use the uniform-based reject/accept
-;; (lambda*- (draw-left-right-truncated-standard-normal
+;; (lambda*- (draw-left-right-truncated-standard-normal
;; left width (* (expt left 2) 0.5d0))))))))))
;; ;; truncated on the left
;; (left
@@ -355,7 +359,7 @@ respectively), on the interval [left, \infinity).")
;; (if (<= left 0d0)
;; (lambda* (try ((x (draw-standard-normal)))
;; (<= left x) x))
-;; (lambda* (draw-left-truncated-standard-normal
+;; (lambda* (draw-left-truncated-standard-normal
;; left
;; (truncated-normal-optimal-alpha left))))))
;; ;; truncated on the right, flip
@@ -364,12 +368,12 @@ respectively), on the interval [left, \infinity).")
;; (if (<= left 0d0)
;; (lambda*- (try ((x (draw-standard-normal)))
;; (<= left x) x))
-;; (lambda*- (draw-left-truncated-standard-normal
+;; (lambda*- (draw-left-truncated-standard-normal
;; left
;; (truncated-normal-optimal-alpha left))))))
;; ;; this is a standard normal, no truncation
;; (t (lambda* (draw-standard-normal)))))))
-
+
;;; Lognormal distribution
@@ -380,13 +384,13 @@ respectively), on the interval [left, \infinity).")
(with-doubles (log-mean log-sd)
(assert (plusp log-sd))
(make :log-mean log-mean :log-sd log-sd))
-
+
(mean () (exp (+ log-mean (/ (expt log-sd 2) 2))))
(variance () (let ((sigma^2 (expt log-sd 2)))
(* (1- (exp sigma^2))
(exp (+ (* 2 log-mean) sigma^2)))))
(log-pdf (x &optional ignore-constant?)
- (maybe-ignore-constant ignore-constant?
+ (maybe-ignore-constant ignore-constant?
(with-doubles (x)
(let ((log-x (log x)))
(- (/ (expt (- log-x log-mean) 2)
@@ -459,7 +463,7 @@ and c using the utility function above. "
(declare (optimize (speed 3)))
(declare (double-float d c))
(check-type alpha (double-float 1d0))
- (tagbody
+ (tagbody
top
(let+ (((&values x v) (prog () ; loop was not optimized for some reason
top
@@ -609,7 +613,7 @@ x^(alpha-1)*(1-x)^(beta-1).")
;;; Discrete distribution.
-;;;
+;;;
;;; ?? The implementation may be improved speedwise with declarations and
;;; micro-optimizations. Not a high priority. However, converting arguments
;;; to double-float provided a great speedup, especially in cases when the
@@ -653,7 +657,7 @@ x^(alpha-1)*(1-x)^(beta-1).")
;; save what's needed
(make :probabilities probabilities :prob prob :alias alias :n-double n-double))
(mean ()
- (iter
+ (iter
(for p :in-vector probabilities :with-index i)
(summing (* p i))))
(variance ()
@@ -670,11 +674,10 @@ x^(alpha-1)*(1-x)^(beta-1).")
(iter
(for p :in-vector probabilities :to i)
(summing p))
- (cumulative-sum probabilities
+ (cumulative-sum probabilities
:result-type 'double-float-vector)))
(draw (&key)
(multiple-value-bind (j p) (floor (random n-double))
(if (<= p (aref prob j))
j
(aref alias j)))))
-
View
1  src/utilities.lisp
@@ -154,4 +154,3 @@ them as a VECTOR-DOUBLE-FLOAT. Vector is always copied."
(assert (<= 0 x) (x) "Element is not positive.")
(/ x sum))
vector)))
-
View
8 tests/univariate.lisp
@@ -60,12 +60,16 @@
(ensure-same (pdf rv -0.7) 0)
(ensure-same (pdf rv -0.3) 0.5515669)
(ensure-same (cdf rv -0.7) 0)
- (ensure-same (cdf rv -0.3) 0.1063703))
+ (ensure-same (quantile rv 0) -0.5)
+ (ensure-same (cdf rv -0.3) 0.1063703)
+ (ensure-same (quantile rv 0.1063703) -0.3))
(let ((rv (r-truncated-normal -0.5 nil -0.4 2.5)))
(ensure-same (pdf rv -0.7) 0)
(ensure-same (pdf rv -0.3) 0.3090382)
(ensure-same (cdf rv -0.7) 0)
- (ensure-same (cdf rv -0.3) 0.06184061)))
+ (ensure-same (cdf rv -0.3) 0.06184061)
+ (ensure-same (quantile rv 0) -0.5)
+ (ensure-same (quantile rv 0.06184061) -0.3)))
(addtest (cl-random-tests)
;; LEFT

No commit comments for this range

Something went wrong with that request. Please try again.