You can edit this page using wiki syntax for markup.

Article contents:== Outdated egg! This is an egg for CHICKEN 4, the unsupported old release. You're almost certainly looking for [[/eggref/5/probdist|the CHICKEN 5 version of this egg]], if it exists. If it does not exist, there may be equivalent functionality provided by another egg; have a look at the [[https://wiki.call-cc.org/chicken-projects/egg-index-5.html|egg index]]. Otherwise, please consider porting this egg to the current version of CHICKEN. [[tags:egg]] == probdist Probability density functions. [[toc:]] == Usage * (require-extension sampled-pdf) * (require-extension normal-pdf) == Documentation {{probdist}} is a library of routines to compute univariate or multivariate normal probability function for a given mean and covariance; and to approximate an arbitrary probability distribution as a weighted set of samples (sampled PDF). This code is based on code in the {{dysii}} project by Lawrence Murray. The sampled PDF algorithm is based on the paper by Kitagawa: Monte Carlo Filter and Smoother for Non-Gaussian Nonlinear State Space Models_, Journal of Computational and Graphical Statistics, 1996, 5, 1-25. === Normal PDF (define-datatype normal-pdf normal-pdf? (NormalPDF S mu sigma R sigmaInv sigmaDet ZI mu-zero? )) A representation of a normal PDF: ; {{S}} : number of dimensions ; {{mu}} : expectation ; {{sigma}} : covariance ; {{R}} : upper triangular Cholesky decomposition of covariance matrix ; {{sigmaInv}} : inverse of covariance matrix ; {{sigmaDet}} : determinant of covariance matrix ; {{ZI}} : the constant factor {{(2*pi)^(-S/2) / prod(diag(R))}} ; {{mu-zero?}} : true if expectation is zero <procedure>make-normal-pdf:: S * MU * SIGMA -> NORMAL-PDF</procedure> Creates a new normal PDF object with the specified dimension, expectation, and covariance. If the dimension {{S}} is 1, then the expectation {{MU}} and the covariance {{SIGMA}} must be scalars; otherwise, they must be SRFI-4 {{f64}} vectors. In the latter case, {{MU}} must be a vector of size {{S}}, and {{SIGMA}} must be a matrix of size {{S x S}}. <procedure>normal-pdf:density:: NORMAL-PDF * X -> DENSITY</procedure> Computes the density of the distribution at the given point {{X}}. If the dimension {{S}} of the distribution is 1, then {{X}} must be a scalar, otherwise it must be an SRFI-4 {{f64}} vector of length {{S}}. <procedure>normal-pdf:sample:: NORMAL-PDF * X -> SAMPLE</procedure> Sample from the distribution. Let {{X}} be a sample from a standard normal distribution. Then the sample is {{SAMPLE = MU + R*X}}, where {{R}} is the upper triangular Cholesky decomposition of the covariance matrix. <procedure>normal-pdf:expectation:: NORMAL-PDF -> MU</procedure> Return the expectation value of the distribution (scalar or f64vector). <procedure>normal-pdf:covariance:: NORMAL-PDF -> SIGMA</procedure> Return the covariance value of the distribution (scalar or f64vector). === Sampled PDF (define-datatype sampled-pdf sampled-pdf? (SampledPDF N W WXS MU SIGMA)) A representation of a weighted sampled PDF: ; {{N}} : sample set size ; {{W}} : total weight of the sample set ; {{wxs}} : weighted samples (see below) ; {{mu}} : expectation ; {{sigma}} : covariance (define-datatype weighted-samples weighted-samples? (WeightedSamples S senv make-senv)) A representation of weighted samples: ; {{S}} : number of dimensions ; {{senv}} : an ordered dictionary data structure that follows the API of [[rb-tree|rb-tree]] ; make-senv : a procedure to create an empty samples dictionary <procedure>make-sampled-pdf:: S * MAKE-SENV * XS [* WS * XCAR * XCDR * XNULL?] -> SAMPLED-PDF</procedure> Creates a new sampled PDF object with the specified dimension, samples environment creation procedure that follows the API of [[rb-tree]], and weighted samples. The samples can be specified in one of two ways. If both {{XS}} and {{WS}} are provided, then {{XS}} is a sequence that contains the samples, and {{WS}} is a sequence that contains the corresponding weights. If only {{XS}} is provided, or {{WS}} is false, then {{XS}} must consist of cons cells, where the car is the weight, and the cdr is the sample. Optional arguments {{XCAR XCDR XNULL?}} could be used for sequences other than lists (e.g. streams). <procedure>sampled-pdf:expectation:: SAMPLED-PDF -> MU</procedure> Return the expectation value of the distribution (scalar or f64vector). <procedure>sampled-pdf:covariance:: SAMPLED-PDF -> SIGMA</procedure> Return the covariance value of the distribution (scalar or f64vector). <procedure>sampled-pdf:find-sample:: U * SAMPLED-PDF -> X</procedure> Given a number {{u in [0,1]}}, returns the sample {{x(i)}} such that {{\sum_{j=1}^{i-1} w_j < Wu <= \sum_{j=1}^{i} w_j}}, where {{W}} is the total weight of the sample set. <procedure>sampled-pdf:normalize:: SAMPLED-PDF -> SAMPLED-PDF</procedure> Normalizes the sample weights to sum to 1. <procedure>sampled-pdf:resample:: MAKE-SENV * SAMPLED-PDF [ * M] -> SAMPLED-PDF</procedure> Resamples the distribution. This produces a new approximation of the same distribution using a set of equally weighted sample points. Sample points are selected using the deterministic resampling method given in the appendix of Kitagawa. Optional argument {{M}} is the number of samples to take for the new distribution. If not given, defaults to the number of samples in the existing distribution. == Examples ;; ;; Multivariate test of normal-pdf. ;; ;; This test creates a 3-dimensional normal distribution with a ;; random mean and variance. It then samples from this distribution ;; and calculates the sample mean and variance for comparison. ;; (require-extension srfi-1) (require-extension srfi-4) (require-extension srfi-4-utils) (require-extension matrix-utils) (require-extension blas) (require-extension normal-pdf) (require-extension random-mtzig) (define-utility-matrices f64) ;; Number of samples (define N 100000) ;; Number of dimensions (sample size) (define S 3) (define (f64vector-scale a x) (let ((n (f64vector-length x))) (dscal n a x))) (define (f64vector-sum x y) (let ((n (f64vector-length x))) (daxpy n 1.0 x y))) (define (f64vector-sub x y) (let ((n (f64vector-length x))) (daxpy n -1.0 y x))) (define (f64vector-mul a b) (f64vector-map (lambda (v1 v2) (* v1 v2)) a b)) (define f64matrix-map (make-matrix-map f64vector-ref f64vector-set!)) (define (upper M N A) (f64matrix-map RowMajor M N A (lambda (i j v) (if (>= j i) v 0.0)))) ;; ;; Computes the arithmetic mean of a list of f64 vectors using the ;; recurrence relation mean_(n) = mean(n-1) + (v[n] - ;; mean(n-1))/(n+1) ;; (define (mean s vs) (let loop ((i 0) (vs vs) (mean (make-f64vector s 0.0))) (if (null? vs) mean (loop (fx+ 1 i) (cdr vs) (let ((d (f64vector-sub (car vs) mean))) (f64vector-sum mean (f64vector-scale (/ 1 (+ 1 i)) d))))))) (define (covariance s vs mean) (let* ((cov (matrix-zeros S S))) (let ((ds (map (lambda (x) (f64vector-sub x mean)) vs))) (let ((n (fold (lambda (d n) (dsyr! RowMajor Upper s 1.0 d cov) (fx+ 1 n)) 0 ds))) (f64vector-scale (/ 1 (- n 1)) cov))))) ;; Given a vector of values in the range [0..1], scale all values in ;; the vector to the specified range. (define (f64vector-urange x lo hi) (if (< lo hi) (let ((d (- hi lo))) (f64vector-map (lambda (x) (+ lo (* d x))) x)))) (define rng (random-mtzig:init)) ;; distribution parameters (define mu (f64vector-urange (random-mtzig:f64vector-randu! S rng) -5 5)) (define sigma (let ((x (f64vector-urange (random-mtzig:f64vector-randu! (* S S) rng) 0 5))) (dgemm! RowMajor NoTrans Trans S S S 1.0 x x 0.0 (make-f64vector (* S S))))) (define gpdf (begin (print "mu = " mu) (print "sigma = " sigma) (make-normal-pdf S mu sigma))) (define input (list-tabulate N (lambda (i) (random-mtzig:f64vector-randn! S rng)))) (define samples (let loop ((i 0) (input input) (samples (list))) (if (null? input) samples (let ((samples (cons (normal-pdf:sample gpdf (car input)) samples))) (loop (fx+ 1 i) (cdr input) samples))))) (if (<= N 100) (print "input = " input)) (if (<= N 100) (print "samples = " samples)) ;; test for mean and covariance equality to one decimal place (equal? (f64vector-map (lambda (x) (truncate (* 10 x))) mu) (f64vector-map (lambda (x) (truncate (* 10 x))) (mean S samples))) (equal? (f64vector-map (lambda (x) (truncate (* 10 x))) (upper S S sigma)) (f64vector-map (lambda (x) (truncate (* 10 x))) (covariance S samples (mean S samples)))) == About this egg === Author [[/users/ivan-raikov|Ivan Raikov]] === Version history ; 1.7 : Updated use of blas ; 1.6 : Documentation converted to wiki format ; 1.5 : Ported to Chicken 4. ; 1.4 : Added some error checking in make-normal-pdf and a record printer for the normal-pdf datatype. ; 1.3 : Build script updated for better cross-platform compatibility ; 1.2 : More metafile fixes ; 1.1 : Metafile fixes ; 1.0 : Initial release === License Copyright 2007-2011 Ivan Raikov and the Okinawa Institute of Science and Technology. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program 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 General Public License for more details. A full copy of the GPL license can be found at <http://www.gnu.org/licenses/>. Description of your changes:

I would like to authenticate

Spam control

What do you get when you multiply 5 by 6?