#lang scheme/base ;;; PLT Scheme Science Collection ;;; random-distributions/lognormal.ss ;;; Copyright (c) 2004-2008 M. Douglas Williams ;;; ;;; This library is free software; you can redistribute it and/or ;;; modify it under the terms of the GNU Lesser General Public ;;; License as published by the Free Software Foundation; either ;;; version 2.1 of the License, or (at your option) any later version. ;;; ;;; This library is distributed in the hope that it will be useful, ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ;;; Lesser General Public License for more details. ;;; ;;; You should have received a copy of the GNU Lesser General Public ;;; License along with this library; if not, write to the Free ;;; Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA ;;; 02111-1307 USA. ;;; ;;; ------------------------------------------------------------------- ;;; ;;; Yhis module implements the lognormal distribution. It is based on ;;; the Random Number Distributions in the GNU Scientific Library. ;;; ;;; Version Date Description ;;; 1.0.0 09/23/04 Marked as ready for Release 1.0 Added ;;; contracts for functions. (Doug Williams) ;;; 2.0.0 11/19/07 Added unchecked versions of functions and ;;; getting ready for PLT Scheme 4.0 (Doug Williams) ;;; 3.0.0 06/09/08 Changes required for V4.0. (Doug Williams) (require (lib "contract.ss")) (provide (rename-out (random-lognormal unchecked-random-logmormal) (lognormal-pdf unchecked-lognormal-pdf) (lognormal-cdf unchecked-lognormal-cdf))) (provide/contract (random-lognormal (case-> (-> random-source? real? (>=/c 0.0) real?) (-> real? (>=/c 0.0) real?))) (lognormal-pdf (-> real? real? (>=/c 0.0) (>=/c 0.0))) (lognormal-cdf (-> real? real? (>=/c 0.0) (real-in 0.0 1.0)))) (require "../math.ss") (require "../random-source.ss") (require "../special-functions/error.ss") ;;; random-lognormal: random-source x real x real -> real ;;; random-lognormal: real x real -> real ;;; This function returns a random variate from the lognormal ;;; distribution with parameters mu and sigma. (define random-lognormal (case-lambda ((r mu sigma) (let ((u 0.0) (v 0.0) (r2 0.0) (normal 0.0)) (let loop () (set! u (+ -1.0 (* 2.0 (unchecked-random-uniform r)))) (set! v (+ -1.0 (* 2.0 (unchecked-random-uniform r)))) (set! r2 (+ (* u u) (* v v))) (when (> r2 1.0) (loop))) (set! normal (* u (sqrt (/ (* -2.0 (log r2)) r2)))) (exp (+ (* sigma normal) mu)))) ((mu sigma) (random-lognormal (current-random-source) mu sigma)))) ;;; lognormal-pdf: real x real x real -> real ;;; This function computes the probability density p(x) at x for a ;;; lognormal distribution with parameters mu and sigma. (define (lognormal-pdf x mu sigma) (if (<= x 0.0) 0.0 (let* ((u (/ (- (log x) mu) sigma)) (p (* (/ 1.0 (* x (abs sigma) (sqrt (* 2.0 pi)))) (exp (/ (- (* u u)) 2.0))))) p))) ;;; lognormal-cdf: real x real x real -> real ;;; This function computes the cummulative density d(x) at x for a ;;; lognormal distribution with parameters mu and sigma. (define (lognormal-cdf x mu sigma) (if (<= x 0.0) 0.0 (* 0.5 (+ 1.0 (unchecked-erf (/ (- (log x) mu) (* sigma sqrt2)))))))