/[sapaclisp]/sapaclisp/random.lisp
ViewVC logotype

Contents of /sapaclisp/random.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (vendor branch)
Mon May 9 20:25:06 2005 UTC (8 years, 11 months ago) by mmommer
Branch: MAIN, T
CVS Tags: initial_checkin, HEAD
Changes since 1.1: +0 -0 lines
Initial Checkin
1 ;;;-*- Mode: LISP; Package: :SAPA; Syntax: COMMON-LISP -*-
2 ;-----------------------------------------------------------------------------
3 ; (c) 1993, Donald B. Percival <dbp@apl.washington.edu>
4 ;
5 ; This code is licensed under the terms of the modified BSD license
6 ; ("sans advertising clause"). See the file COPYING for details.
7 ;
8 ; Comments about this software can be addressed to dbp@apl.washington.edu
9 ;-----------------------------------------------------------------------------
10 ;
11 ; random.lisp
12 ;
13 ; a collection of Lisp functions to simulate stationary random processes ...
14 ; Note: before compiling and loading random.lisp,
15 ; you should compile and load (in the order listed)
16 ; sapa-package.lisp, utilities.lisp and dft-and-fft.lisp
17 ;
18 ;-----------------------------------------------------------------------------
19 (if (not (find-package :SAPA))
20 (defpackage :SAPA (:USE :COMMON-LISP)))
21
22 (in-package :SAPA)
23
24 (export '(;;; functions to generate normally distributed random deviates ...
25 ranorm
26 ranorms
27
28 ;;; white noise with different distributions ...
29 generate-white-noise
30
31 ;;; normally distributed moving average and autoregressive processes
32 generate-ma-time-series
33 generate-ar-time-series
34 step-down-Levinson-Durbin-recursions
35
36 ;;; simulate a stationary process using frequency domain techniques
37 simulate-time-series-from-sdf
38 acvs-for-time-series-simulated-from-sdf
39 ))
40
41 ;-------------------------------------------------------------------------------
42 (defconstant +sapa-sqrt-8-over-e+ (sqrt (/ 8.0d0 (exp 1.0d0))))
43 (defconstant +sapa-4-time-exp-of-1-over-4+ (* 4.0d0 (exp 0.25d0)))
44 (defconstant +sapa-4-time-exp-of-minus-1-point-35+ (* 4.0d0 (exp (- 1.35d0))))
45
46 ;-------------------------------------------------------------------------------
47 ;-------------------------------------------------------------------------------
48 ;;; The functions ranorm
49 ;;; ranorms
50 ;;; generate one or more uncorrelated normally distributed deviates.
51 ;-------------------------------------------------------------------------------
52 ;-------------------------------------------------------------------------------
53 (defun ranorm
54 ()
55 "returns a random deviate from a normal distribution
56 with zero mean and unit variance"
57 (let ((u (random 1.0)))
58 (cond
59 ((= u 0.0) (ranorm)) ;bad choice!
60 (t
61 (let* ((x (/ (* (- (random 1.0) 0.5)
62 +sapa-sqrt-8-over-e+) u))
63 (xs (* x x)))
64 (cond
65 ((<= xs (- 5.0 (* u +sapa-4-time-exp-of-1-over-4+)))
66 x) ;done
67 ((>= xs (+ 1.4 (/ +sapa-4-time-exp-of-minus-1-point-35+ u)))
68 (ranorm)) ;do it again
69 ((<= xs (- (* 4.0 (log u))))
70 x)
71 (t
72 (ranorm))))))))
73
74 #|
75 ;;; each evaluation of ranorm produces a different random deviate,
76 ;;; so don't expect to get the same result!
77 (ranorm) ;==> -0.2719851788703296
78 |#
79
80 ;-------------------------------------------------------------------------------
81 (defun ranorms
82 (n
83 &key
84 (result (make-array n)))
85 "given:
86 [1] n (required)
87 ==> number of normal random deviates
88 to be generated
89 [2] result (keyword; vector of length n)
90 <== vector to hold random deviates
91 returns:
92 [1] result, a vector with n normal random deviates
93 (i.e., a realization of length n of a white
94 noise process from a normal distribution)"
95 (dotimes (i n result)
96 (setf (svref result i) (ranorm))))
97
98 #|
99 (ranorms 3)
100 ;==> #(-1.1128679745343053 1.4875994154684462 -0.06612547414267966)
101 |#
102
103 ;-------------------------------------------------------------------------------
104 ;-------------------------------------------------------------------------------
105 ;;; The function generate-white-noise
106 ;;; generate a sequence of white noise with a specified distribution.
107 ;-------------------------------------------------------------------------------
108 ;-------------------------------------------------------------------------------
109 (defun generate-white-noise
110 (n
111 &key
112 (distribution :normal)
113 (result (make-array n)))
114 "given
115 [1] n (required)
116 ==> a sample size
117 [2] distribution (keyword; :normal)
118 ==> either a keyword or a function with no arguments
119 that returns a random deviate from a particular
120 distribution. Choices are
121 :binary :cauchy :chi-square-2
122 :double-exponential :exponential
123 :extreme-value :Gaussian (same as :normal)
124 :logistic :lognormal :normal :uniform
125 [3] result (keyword; vector of length n)
126 <== a sequence in which results are to be stored
127 return
128 [1] result, containing n samples from a white noise process
129 with a distribution specified by the keyword distribution"
130 (let ((generating-function
131 (if (keywordp distribution)
132 (case distribution
133 (:binary #'(lambda ()
134 (if (<= (random 1.0) 0.5)
135 0 1)))
136 (:cauchy #'(lambda ()
137 (tan (* pi (- (random 1.0) 0.5)))))
138 (:chi-square-2 #'(lambda ()
139 (* 2 (- (log (- 1.0 (random 1.0)))))))
140 (:double-exponential #'(lambda ()
141 (let ((U (random 1.0)))
142 (if (<= U 0.5)
143 (log (* 2 U))
144 (- (log (* 2 (- 1.0 U))))))))
145 (:exponential #'(lambda ()
146 (- (log (- 1.0 (random 1.0))))))
147 (:extreme-value #'(lambda ()
148 (log (- (log (- 1.0 (random 1.0)))))))
149 (:Gaussian #'ranorm)
150 (:logistic #'(lambda ()
151 (let ((U (random 1.0)))
152 (log (/ U (- 1.0 U))))))
153 (:lognormal #'(lambda ()
154 (exp (ranorm))))
155 (:normal #'ranorm)
156 (:uniform #'(lambda () (random 1.0)))
157 (otherwise
158 (error "~&generate-white-noise: unknown distribution keyword (~A)"
159 distribution)))
160 distribution)))
161 (dotimes (i n result)
162 (setf (elt result i) (funcall generating-function)))))
163
164 #|
165 (generate-white-noise 10 :distribution :binary)
166 ;==> #(1 0 1 0 1 0 0 0 0 0)
167 (generate-white-noise 3 :distribution :cauchy)
168 ;==> #(-27.51099440214296 -0.6075812207412569 0.36550059169755994)
169 (generate-white-noise 3 :distribution :lognormal)
170 ;==> #(0.5822708648379962 3.833558848750849 0.3592735989502123)
171 |#
172
173 ;-------------------------------------------------------------------------------
174 ;-------------------------------------------------------------------------------
175 ;;; The functions generate-ma-time-series
176 ;;; generate-ar-time-series
177 ;;; step-down-Levinson-Durbin-recursions
178 ;;; can be used to generate realizations of specified lengths from
179 ;;; normally distributed moving average (MA) or autoregressive processes (AR)
180 ;;; (these processes are discussed in Chapter 2 and 9 of the SAPA book).
181 ;-------------------------------------------------------------------------------
182 ;-------------------------------------------------------------------------------
183 (defun generate-ma-time-series
184 (coeffs
185 variance
186 sample-size
187 &key
188 (process-variance-p t)
189 (result (make-array sample-size)))
190 "given
191 [1] coeffs (required)
192 ==> sequence of length q with MA coefficients:
193 x_t = e_t - coeffs_0*e_{t-1}
194 - coeffs_1*e_{t-2}
195 - ... - coeffs_{q-1}*e_{t-q}
196 (see Equation (43a) in the SAPA book)
197 [2] variance (required)
198 ==> process variance or innovations variance
199 (see keyword process-variance-p)
200 [3] sample-size (required)
201 ==> length of generated time series
202 [4] process-variance-p (keyword; t)
203 ==> if t, variance is taken to be process variance
204 if nil, variance is taken to be innovations variance
205 [5] result (keyword; vector of length sample-size)
206 <== vector with simulated series
207 generates realization of zero mean normally distributed MA(q) process and
208 returns
209 [1] vector of length sample-size with realization"
210 (let ((q+1 (1+ (length coeffs)))
211 (sd (sqrt (if process-variance-p
212 (/ variance (1+ (sum-of-squares coeffs)))
213 variance)))
214 (scratch-1 (make-array (1+ (length coeffs))))
215 (scratch-2 (make-array (1+ (length coeffs)))))
216 #+mcl(declare (dynamic-extent scratch-1 scratch-2))
217 (generate-white-noise q+1 :result scratch-1)
218 (a*x! sd scratch-1)
219 (dotimes (i (1- q+1))
220 (setf (aref scratch-2 (1+ i)) (- (elt coeffs i))))
221 (setf (aref scratch-2 0) 1.0)
222 (dotimes (i sample-size result)
223 (setf (aref result i)
224 (dot-product scratch-2 scratch-1))
225 (circular-shift-sequence scratch-1 :result scratch-1)
226 (setf (aref scratch-1 0) (* sd (ranorm))))))
227
228 #|
229 ;;; See top plot of Figure 44 of the SAPA book.
230 (generate-ma-time-series #(1.0) 1.0 4 :process-variance-p nil)
231 ;==> #(-0.555869655695869 -0.6755180985339049 0.3536345814801179 0.7421943411372343)
232 ;;; See bottom plot of Figure 44 of the SAPA book.
233 (generate-ma-time-series #(-1.0) 1.0 4 :process-variance-p nil)
234 #(0.2105828208289749 1.1611109339530934 1.7492356283253514 0.07032022771512858)
235 |#
236
237 ;-------------------------------------------------------------------------------
238 (defun generate-ar-time-series
239 (coeffs
240 variance
241 sample-size
242 &key
243 (process-variance-p t)
244 (list-of-lower-order-phi nil)
245 (list-of-lower-order-pev nil)
246 (result (make-array
247 sample-size
248 :initial-element 0.0)
249 result-supplied-p))
250 "given
251 [1] coeffs (required)
252 ==> sequence of length p with AR coefficients:
253 x_t = coeffs_0*x_{t-1} + coeffs_1*x_{t-2}
254 + ... + coeffs_{p-1}*x_{t-p} + e_t
255 (see Equation (392a) in the SAPA book;
256 the coefficients can be real or complex-valued)
257 [2] variance (required)
258 ==> process variance or innovations variance
259 (see keyword process-variance-p)
260 [3] sample-size (required)
261 ==> length of generated time series
262 [4] process-variance-p (keyword; t)
263 ==> if t, variance is taken to be process variance
264 if nil, variance is taken to be innovations variance
265 [5] list-of-lower-order-phi (keyword; nil)
266 ==> to bypass call to step-down-Levinson-Durbin-recursions
267 (useful for multiple realizations)
268 [6] list-of-lower-order-pev (keyword; nil)
269 ==> to bypass call to step-down-Levinson-Durbin-recursions
270 [7] result (keyword; vector of length sample-size)
271 <== vector with simulated series
272 generates realization of zero mean normally distributed AR(p) process and
273 returns
274 [1] result, vector of length sample-size with realization
275 [2] list-of-lower-order-phi (can be used on subsequent calls
276 to this function to bypass call to
277 step-down-Levinson-Durbin-recursions)
278 [3] list-of-lower-order-pev (can be also used on subsequent calls)
279 ---
280 Note: this function generates the proper stationary initial conditions"
281 (if result-supplied-p (fill result 0.0 :end sample-size))
282 (if (not list-of-lower-order-pev)
283 (multiple-value-setq (list-of-lower-order-phi list-of-lower-order-pev)
284 (step-down-Levinson-Durbin-recursions
285 coeffs variance
286 :process-variance-p process-variance-p)))
287 (let ((p (length coeffs))
288 (sd (sqrt (nth 0 list-of-lower-order-pev))))
289 (dotimes (i sample-size (values result
290 list-of-lower-order-phi
291 list-of-lower-order-pev))
292 (cond
293 ((< i p)
294 (dotimes (j i)
295 (incf (aref result i)
296 (* (aref result (- i j 1))
297 (aref (nth (1- i) list-of-lower-order-phi) j))))
298 (incf (aref result i) (* sd (ranorm)))
299 (setf sd (sqrt (nth (1+ i) list-of-lower-order-pev))))
300 (t
301 (dotimes (j p)
302 (incf (aref result i) (* (aref result (- i j 1))
303 (elt coeffs j))))
304 (incf (aref result i) (* sd (ranorm))))))))
305
306 #|
307 ;;; AR(2) process described by Equation (45) of the SAPA book:
308 (generate-ar-time-series #(0.75 -0.5) 1.0 4 :process-variance-p nil)
309 ;==> #(-0.5263377203675295 -1.5334622586392335 -1.628199801947855 1.0935751830627396)
310 ; (#(0.5) #(0.75 -0.5))
311 ; (1.7777777777777777 1.3333333333333333 1.0)
312
313 ;;; AR(4) process described by Equation (46a) of the SAPA book:
314 (generate-ar-time-series #(2.7607 -3.8106 2.6535 -0.9238)
315 1.0 4 :process-variance-p nil)
316 ;==> #(-6.887266139675448 -31.09831681568327 -34.70711181816937 -17.248058702547237)
317 (#(0.7164772070165697)
318 #(1.4197722065109775 -0.9816013581547822)
319 #(2.1105749802378733 -1.9807672315209488 0.7037508332562511)
320 #(2.7607 -3.8106 2.6535 -0.9238))
321 ; (761.7172900314962 370.6976500615111 13.515181723106767 6.821582066770185 1.0)
322 |#
323
324 ;-------------------------------------------------------------------------------
325 (defun step-down-Levinson-Durbin-recursions
326 (coeffs
327 variance
328 &key
329 (process-variance-p t))
330 "given
331 [1] coeffs (required)
332 ==> sequence of length p with AR coefficients:
333 x_t = coeffs_0*x_{t-1} + coeffs_1*x_{t-2}
334 + ... + coeffs_{p-1}*x_{t-p} + e_t
335 (see Equation (392a) in the SAPA book;
336 the coefficients can be real or complex-valued)
337 [2] variance (required)
338 ==> process variance or innovations variance
339 (see keyword process-variance-p)
340 [3] process-variance-p (keyword; t)
341 ==> if t, variance is taken to be process variance
342 if nil, variance is taken to be innovations variance
343 computes best linear prediction coefficients of orders 1, 2, ..., p-1
344 and prediction error variances of orders 0 (process variance), 1, ..., p and
345 returns
346 [1] list of vectors with best linear prediction coefficients
347 going from order 1 to order p;
348 [2] list of prediction error variances going from order 0 to order p
349 ---
350 Note: see item [4] of the Comments and Extensions
351 to Section 9.4 of the SAPA book. The values returned by
352 this function can be used to set the keyword parameters
353 list-of-lower-order-phi and list-of-lower-order-pev in
354 the function generate-ar-time-series"
355 (let* ((p (length coeffs))
356 (p-1 (1- p))
357 (list-of-lower-order-phi `(,coeffs))
358 (list-of-lower-order-pev `(,variance))
359 k-1 phi-k phi-k-1 den)
360 ;;; generate lower order coefficients via step-down Levinson-Durbin
361 (dotimes (i p-1)
362 (setf k-1 (- p-1 i)
363 phi-k (car list-of-lower-order-phi)
364 phi-k-1 (car (push (make-array k-1) list-of-lower-order-phi))
365 den (- 1.0 (realpart
366 (* (aref phi-k k-1) (conjugate (aref phi-k k-1))))))
367 (dotimes (j k-1)
368 (setf (aref phi-k-1 j) (/ (+ (aref phi-k j)
369 (* (aref phi-k k-1)
370 (conjugate (aref phi-k (- k-1 j 1)))))
371 den))))
372 (values (if (plusp p) list-of-lower-order-phi)
373 (cond
374 (process-variance-p
375 (dotimes (i p (reverse list-of-lower-order-pev))
376 (push (* (car list-of-lower-order-pev)
377 (- 1.0
378 (realpart
379 (* (aref (nth i list-of-lower-order-phi) i)
380 (conjugate
381 (aref (nth i list-of-lower-order-phi) i))))))
382 list-of-lower-order-pev)))
383 (t
384 (dotimes (i p list-of-lower-order-pev)
385 (push (/ (car list-of-lower-order-pev)
386 (- 1.0
387 (realpart
388 (* (aref (nth (- p-1 i) list-of-lower-order-phi)
389 (- p-1 i))
390 (conjugate
391 (aref (nth (- p-1 i) list-of-lower-order-phi)
392 (- p-1 i)))))))
393 list-of-lower-order-pev)))))))
394
395 #|
396 ;;; AR(2) process described by Equation (45) of the SAPA book:
397 (step-down-Levinson-Durbin-recursions #(0.75 -0.5)
398 1.0 :process-variance-p nil)
399 ;==> (#(0.5) #(0.75 -0.5))
400 ; (1.7777777777777777 1.3333333333333333 1.0)
401
402 ;;; AR(4) process described by Equation (46a) of the SAPA book:
403 (step-down-Levinson-Durbin-recursions #(2.7607 -3.8106 2.6535 -0.9238)
404 1.0 :process-variance-p nil)
405 ;==> (#(0.7164772070165697)
406 #(1.4197722065109775 -0.9816013581547822)
407 #(2.1105749802378733 -1.9807672315209488 0.7037508332562511)
408 #(2.7607 -3.8106 2.6535 -0.9238))
409 ; (761.7172900314962 370.6976500615111 13.515181723106767 6.821582066770185 1.0)
410 |#
411
412 ;-------------------------------------------------------------------------------
413 ;-------------------------------------------------------------------------------
414 ;;; The functions simulate-time-series-from-sdf
415 ;;; acvs-for-time-series-simulated-from-sdf
416 ;;; can be used to generate realizations of specified lengths
417 ;;; from normally distributed stationary processes with a specified
418 ;;; spectral density function or autocovariance sequence. These functions
419 ;;; use frequency domain techniques. The techniques are discussed
420 ;;; in the paper ``Simulating Gaussian Random Processes
421 ;;; with Specified Spectra'' by Percival in "Computing Science and Statistics",
422 ;;; vol 24, 534--8 (contains Proc. of 24th Symp. on the Interface of Computer
423 ;;; Science and Statistics).
424 ;-------------------------------------------------------------------------------
425 ;-------------------------------------------------------------------------------
426 (defun simulate-time-series-from-sdf
427 (n-series
428 sdf
429 &key
430 (sampling-time 1.0)
431 (n-total
432 (* 4 (next-power-of-2 n-series)))
433 (result (make-array n-series)))
434 "given
435 [1] n-series (required)
436 ==> length of time series to be simulated
437 [2] sdf (required)
438 ==> spectral density function (sdf) defined for
439 0 <= f <= 1/(2.0 sampling-time) -- the
440 sdf is assumed to be two-sided and symmetric
441 about f = 0
442 [3] sampling-time (keyword; 1.0)
443 ==> the assumed sampling time (delta t)
444 [4] n-total (keyword; 4 * next power of 2 for n-series)
445 ==> a power of 2 controlling degree of accuracy
446 of the approximation (the larger, the better --
447 see Percival, 1992, for details)
448 [5] result (keyword; vector of length n-series)
449 <== a vector to contain simulated time series
450 returns
451 [1] a vector of length n-series generated from sdf
452 ---
453 Note: the method used here is an approximate frequency domain
454 technique; in particular, it will NOT simulate the DC component
455 correctly, so beware of using this function in studies
456 where the process mean is of important"
457 ;;; This assertion fails if N is not a power of two
458 (assert n-total
459 ()
460 "n-total = ~D not a power of 2"
461 n-total)
462 (let* ((n-sampling-time (float (* n-total sampling-time)))
463 (vector-for-fft (make-array n-total))
464 (j 0)
465 (n-j n-total)
466 (n-over-2 (/ n-total 2))
467 (1-over-sqrt-n-sampling-time (/ (sqrt (* n-total sampling-time)))))
468 (setf (svref vector-for-fft 0)
469 (* (sqrt (funcall sdf 0.0))
470 (ranorm))
471 (svref vector-for-fft n-over-2)
472 (* (sqrt (funcall sdf (/ 0.5 sampling-time)))
473 (ranorm)))
474 (dotimes (i (1- n-over-2))
475 (incf j) (decf n-j)
476 (setf (svref vector-for-fft n-j)
477 (conjugate
478 (setf (svref vector-for-fft j)
479 (* (sqrt (* 0.5 (funcall sdf (/ j n-sampling-time))))
480 (complex (ranorm) (ranorm)))))))
481 (fft! vector-for-fft)
482 (dotimes (i n-series result)
483 (setf (svref result i)
484 (* (realpart (svref vector-for-fft i))
485 1-over-sqrt-n-sampling-time)))))
486
487 #|
488 (simulate-time-series-from-sdf
489 4
490 #'(lambda (f) (1+ f)))
491 ;==> #(-0.775118374075842 -0.13513773684982933 -1.4345372479995293 0.31204493073149375)
492 |#
493
494 ;-------------------------------------------------------------------------------
495 (defun acvs-for-time-series-simulated-from-sdf
496 (n-series
497 sdf
498 &key
499 (sampling-time 1.0)
500 (n-total (* 4 (next-power-of-2 n-series)))
501 (result (make-array n-series)))
502 "given
503 [1] n-series (required)
504 ==> length of time series simulated using
505 simulate-time-series-from-sdf
506 (must be a power of 2)
507 [2] sdf (required)
508 ==> spectral density function (sdf) defined for
509 0 <= f <= 1/(2.0 sampling-time) -- the
510 sdf is assumed to be two-sided and symmetric
511 about f = 0
512 [3] sampling-time (keyword; 1.0)
513 ==> the assumed sampling time (delta t)
514 [4] n-total (keyword; 4 * next power of 2 for n-series)
515 ==> a power of 2 controlling degree of accuracy
516 of the approximation (the larger, the better --
517 see Percival, 1992, for details)
518 [5] result (keyword; vector of length n-series)
519 <== a vector to contain simulated time series
520 returns
521 [1] a vector with the theoretical acvs from lag 0
522 to n-series - 1 for a time series generated via a call
523 to simulate-time-series-from-sdf with n-series,
524 sdf, sampling-time and n-total set to the same values"
525 (assert (power-of-2 n-total))
526 (let* ((sigma^2 (make-array n-total))
527 (n-total-sampling-time (* n-total sampling-time))
528 (j 1)
529 (N-j (1- n-total)))
530 (setf (aref sigma^2 0)
531 (/ (funcall sdf 0.0) n-total-sampling-time)
532 (aref sigma^2 (/ n-total 2))
533 (/ (funcall sdf (/ 0.5 sampling-time)) n-total-sampling-time))
534 (dotimes (i (- (/ n-total 2) 1))
535 (setf (aref sigma^2 j)
536 (/ (funcall sdf (/ j n-total-sampling-time))
537 n-total-sampling-time)
538 (aref sigma^2 N-j) (aref sigma^2 j))
539 (incf j)
540 (decf N-j))
541 (fft! sigma^2)
542 (dotimes (i n-series result)
543 (setf (aref result i) (realpart (aref sigma^2 i))))))
544
545 #|
546 (acvs-for-time-series-simulated-from-sdf
547 4
548 #'(lambda (f) (1+ f)))
549 ;==> #(1.25 -0.10263336862925071 0.0 -0.012655581284545123)
550 |#

  ViewVC Help
Powered by ViewVC 1.1.5