Changeset 14646 in project


Ignore:
Timestamp:
05/15/09 04:50:47 (11 years ago)
Author:
Ivan Raikov
Message:

Ported probdist to Chicken 4

Location:
release/4/probdist
Files:
8 edited
1 copied

Legend:

Unmodified
Added
Removed
  • release/4/probdist/trunk/normal-pdf.scm

    r6771 r14646  
    99;;
    1010
    11 (require-extension syntax-case)
    12 (require-extension datatype)
    13 (require-extension matrix-utils)
    14 (require-extension blas)
    15 (require-extension atlas-lapack)
    16 (require-extension srfi-4)
    17 (require-extension srfi-4-utils)
     11(module normal-pdf
    1812
    19 (define-extension normal-pdf)
     13 ( normal-pdf?
     14   make-normal-pdf
     15   normal-pdf:expectation
     16   normal-pdf:covariance
     17   normal-pdf:density
     18   normal-pdf:sample)
     19                 
     20   (import scheme chicken extras srfi-4)
    2021
    21 (define-utility-matrices f64)
    22 
    23 (declare (export normal-pdf?
    24                  make-normal-pdf
    25                  normal-pdf:expectation
    26                  normal-pdf:covariance
    27                  normal-pdf:density
    28                  normal-pdf:sample))
    29 
    30                  
     22   
     23   (require-extension datatype srfi-4 srfi-4-utils blas matrix-utils atlas-lapack)
    3124
    3225(define normal-pdf:order blas:RowMajor)
     
    184177
    185178
     179)
  • release/4/probdist/trunk/probdist-eggdoc.scm

    r7358 r14646  
    2424     (name "probdist")
    2525     (description "Probability density functions.")
    26      (author (url "http://chicken.wiki.br/ivan raikov" "Ivan Raikov"))
     26     (author (url "http://chicken.wiki.br/users/ivan-raikov" "Ivan Raikov"))
    2727
    2828     (history
     29      (version "1.5" "Ported to Chicken 4.")
    2930      (version "1.4" "Added some error checking in make-normal-pdf and a record printer for the normal-pdf datatype.")
    3031      (version "1.3" "Build script updated for better cross-platform compatibility")
     
    265266
    266267     (license
    267       "Copyright 2007 Ivan Raikov and the Okinawa Institute of Science and Technology
     268      "Copyright 2007-2009 Ivan Raikov and the Okinawa Institute of Science and Technology.
    268269
    269270This program is free software: you can redistribute it and/or modify
  • release/4/probdist/trunk/probdist.meta

    r9305 r14646  
     1;;;; -*- Hen -*-
     2
    13((egg "probdist.egg") ; This should never change
    24
     
    1618 ; A list of eggs probdist depends on.
    1719
    18  (needs eggdoc datatype srfi-4 srfi-4-utils blas matrix-utils atlas-lapack)
     20 (needs eggdoc datatype matchable srfi-4-utils blas matrix-utils atlas-lapack)
    1921
    2022 (eggdoc "probdist-eggdoc.scm")
     
    2224 (author "Ivan Raikov")
    2325
    24  (synopsis "Probability distributions"))
     26 (synopsis "Probability distributions."))
  • release/4/probdist/trunk/probdist.setup

    r6772 r14646  
    1 
    2 (define has-exports? (string>=? (chicken-version) "2.310"))
    3 
     1;;;; -*- Hen -*-
    42
    53(define (dynld-name fn)         
    64  (make-pathname #f fn ##sys#load-dynamic-extension))   
    75
    8 (compile -O2 -s
    9          ,@(if has-exports? '(-check-imports -emit-exports normal-pdf.exports) '())
    10          normal-pdf.scm -lchicken -ldl -lm)
     6(compile -O2 -s normal-pdf.scm -j normal-pdf)
     7(compile -O2 -s normal-pdf.import.scm)
    118
    12 (compile -O2 -s
    13          ,@(if has-exports? '(-check-imports -emit-exports sampled-pdf.exports) '())
    14          sampled-pdf.scm -lchicken -ldl -lm)
     9(compile -O2 -s sampled-pdf.scm -j sampled-pdf)
     10(compile -O2 -s sampled-pdf.import.scm)
    1511
    1612(run (csi -qbs probdist-eggdoc.scm > probdist.html))
     
    2218
    2319  ; Files to install for your extension:
    24   `(,(dynld-name "sampled-pdf") ,(dynld-name "normal-pdf") "probdist.html"
    25     ,@(if has-exports? '("sampled-pdf.exports" "normal-pdf.exports") (list)) )
     20  `(,(dynld-name "sampled-pdf") ,(dynld-name "sampled-pdf.import")
     21    ,(dynld-name "normal-pdf") ,(dynld-name "normal-pdf.import") )
    2622
    2723  ; Assoc list with properties for your extension:
    28   '((version 1.4)
    29     (documentation "probdist.html")
    30     ,@(if has-exports? `((exports "sampled-pdf.exports" "normal-pdf.exports")) (list)) ))
     24  '((version 1.5)
     25    (documentation "probdist.html")))
    3126
     27   
  • release/4/probdist/trunk/sampled-pdf.scm

    r6152 r14646  
    2424;;
    2525
    26 (require-extension syntax-case)
    27 (require-extension datatype)
    28 (require-extension srfi-4)
    29 (require-extension srfi-4-utils)
    30 (require-extension blas)
    31 (require-extension matrix-utils)
    32 
    33 (define-extension sampled-pdf)
    34 
    35 (declare (export sampled-pdf?
    36                  weighted-samples?
    37                  SampledPdf WeightedSamples
    38                  make-sampled-pdf
    39                  sampled-pdf:expectation
    40                  sampled-pdf:covariance
    41                  sampled-pdf:find-sample
    42                  sampled-pdf:normalize
    43                  sampled-pdf:resample))
     26(module sampled-pdf
     27
     28         ( sampled-pdf?
     29           weighted-samples?
     30           SampledPDF WeightedSamples
     31           make-sampled-pdf
     32           sampled-pdf:expectation
     33           sampled-pdf:covariance
     34           sampled-pdf:find-sample
     35           sampled-pdf:normalize
     36           sampled-pdf:resample)
    4437                 
     38   (import scheme chicken srfi-1 srfi-4)
     39
     40   
     41   (require-extension datatype matchable srfi-4 srfi-4-utils blas matrix-utils)
    4542
    4643(define-utility-matrices f64)
     
    8178;
    8279(define-datatype sampled-pdf sampled-pdf?
    83   (SampledPdf (N      (lambda (x) (and (integer? x) (> x 0))))
     80  (SampledPDF (N      (lambda (x) (and (integer? x) (> x 0))))
    8481              (W      number?)
    8582              (wxs    weighted-samples?)
     
    154151               (mu     (expectation W wxs))
    155152               (sigma  (covariance W wxs mu)))
    156            (SampledPdf N W wxs mu sigma))))))
     153           (SampledPDF N W wxs mu sigma))))))
    157154
    158155
     
    166163                           "parameter u is not in [0,1]: u = " u))
    167164  (cases sampled-pdf x
    168      (SampledPdf (N W wxs mu sigma)
     165     (SampledPDF (N W wxs mu sigma)
    169166       (cases weighted-samples wxs
    170167         (WeightedSamples (S senv)
    171168            (let ((N  (senv 'size)))
    172               (let ((Wu+x  ((senv 'fold-limit) (lambda (ax) (<= (car ax) 0.0))
    173                             (lambda (wx ax) (cons (- (car Wu) (car wx)) wx))
    174                             (* W u))))
     169              (let ((Wu+x  ((senv 'fold-limit)
     170                            (lambda (ax) (<= (car ax) 0.0))
     171                            (lambda (wx ax) (match ax ((Wu x) (cons (- Wu (car wx)) (cdr wx)))))
     172                            (list (* W u) #f))))
    175173                (cdr Wu+x))))))))
    176174
     
    179177(define (sampled-pdf:normalize x)
    180178  (cases sampled-pdf x
    181      (SampledPdf (N W wxs mu sigma)
     179     (SampledPDF (N W wxs mu sigma)
    182180       (if (= W 1.0) x
    183181           (let ((WI  (/ 1 W)))
     
    205203(define (sampled-pdf:resample make-senv x . rest)
    206204  (cases sampled-pdf x
    207      (SampledPdf (N W wxs mu sigma)
     205     (SampledPDF (N W wxs mu sigma)
    208206       (let-optionals rest ((m N) (alpha 0.5))
    209207         (cases weighted-samples wxs
     
    223221                                 (mu     (expectation rW wxs))
    224222                                 (sigma  (covariance rW wxs mu)))
    225                             (SampledPdf m rW wxs mu sigma)))))))))))
     223                            (SampledPDF m rW wxs mu sigma)))))))))))
    226224
    227225
    228226(define (sampled-pdf:expectation x)                       
    229227  (cases sampled-pdf x
    230      (SampledPdf (N W wxs mu sigma) mu)))
     228     (SampledPDF (N W wxs mu sigma) mu)))
    231229                         
    232230
    233231(define (sampled-pdf:covariance x)                       
    234232  (cases sampled-pdf x
    235      (SampledPdf (N W wxs mu sigma) sigma)))
     233     (SampledPDF (N W wxs mu sigma) sigma)))
    236234
    237235                         
     236)
  • release/4/probdist/trunk/tests/test1.scm

    r6151 r14646  
    88;;
    99
    10 (require-extension srfi-4)
    11 (require-extension normal-pdf)
    12 (require-extension random-mtzig)
    13 (require-extension testbase)
    14 (require-extension testbase-output-compact)
     10(use srfi-4 normal-pdf random-mtzig test)
    1511
    1612;; Number of samples
     
    3632            (loop (fx+ 1 i) (+ sqsum (* delta delta))))
    3733          (* (/ 1 (- n 1)) sqsum)))))
     34
     35
     36;;  distribution parameters
     37(define mu     1.58798)
     38(define sigma  4.55175)
     39
     40(define gpdf (make-normal-pdf 1 mu sigma))
     41
     42(define rng      (random-mtzig:init))
     43(define input    (random-mtzig:f64vector-randn! N rng))
     44(define samples  (make-f64vector N 0))
     45(define density  (make-f64vector N 0))
    3846     
    39 (define-test normal-pdf-test "univariate normal PDF test"
    40   (initial
    41 
    42    ;;  distribution parameters
    43    (define mu     1.58798)
    44    (define sigma  4.55175)
    45 
    46    (define gpdf (make-normal-pdf 1 mu sigma))
    47 
    48    (define rng      (random-mtzig:init))
    49    (define input    (random-mtzig:f64vector-randn! N rng))
    50    (define samples  (make-f64vector N 0))
    51    (define density  (make-f64vector N 0))
     47(test-group "univariate normal PDF test"
    5248
    5349   (let loop ((i 0))
     
    5551         (begin (f64vector-set! samples i (normal-pdf:sample gpdf (f64vector-ref input i)))
    5652                (f64vector-set! density i (normal-pdf:density gpdf (f64vector-ref input i)))
    57                 (loop (fx+ 1 i))))) )
     53                (loop (fx+ 1 i)))))
    5854
    59   ;; test for mean and variance equality to one decimal place
    60   (test/equal 'sample-mean
    61                (truncate (* 10 mu))
    62                (truncate (* 10 (mean samples))))
     55   ;; test for mean and variance equality to one decimal place
     56   (test "sample-mean"
     57         (truncate (* 10 mu))
     58         (truncate (* 10 (mean samples))))
    6359
    64   (test/equal 'sample-variance
    65               (truncate (* 10 sigma))
    66               (truncate (* 10 (variance samples (mean samples))))) )
     60   (test "sample-variance"
     61         (truncate (* 10 sigma))
     62         (truncate (* 10 (variance samples (mean samples)))))
    6763
     64)
    6865
    69 (test::styler-set! normal-pdf-test test::output-style-compact)
    70 
    71 (run-test "probability distribution test")
    7266
    7367
  • release/4/probdist/trunk/tests/test2.scm

    r6151 r14646  
    88;;
    99
    10 (require-extension srfi-1)
    11 (require-extension srfi-4)
    12 (require-extension srfi-4-utils)
    13 (require-extension matrix-utils)
    14 (require-extension blas)
    15 (require-extension normal-pdf)
    16 (require-extension random-mtzig)
    17 (require-extension testbase)
    18 (require-extension testbase-output-compact)
     10(use srfi-1 srfi-4 srfi-4-utils matrix-utils blas random-mtzig
     11     normal-pdf test)
    1912
    2013(define-utility-matrices f64)
     
    6861         (f64vector-scale (/ 1 (- n 1)) cov)))))
    6962                         
    70      
    71 (define-test multi-normal-pdf-test "multivariate normal PDF test"
    72   (initial
    7363
    7464   ;; Given a vector of values in the range [0..1], scale all values in
     
    10393   (if (<= N 100) (print "samples = " samples))
    10494;   (print "density = " (cadr samples+density))
    105 
    106    )
     95     
     96(test-group "multivariate normal PDF test"
    10797
    10898  ;; test for mean and covariance equality to one decimal place
    109   (test/equal 'sample-mean
    110                (f64vector-map (lambda (x) (truncate (* 10 x))) mu)
    111                (f64vector-map (lambda (x) (truncate (* 10 x)))
    112                               (mean S samples)))
     99  (test "sample-mean"
     100        (f64vector-map (lambda (x) (truncate (* 10 x))) mu)
     101        (f64vector-map (lambda (x) (truncate (* 10 x)))
     102                       (mean S samples)))
    113103
    114   (test/equal 'sample-covariance
    115               (f64vector-map (lambda (x) (truncate (* 10 x))) (upper S S sigma))
    116               (f64vector-map (lambda (x) (truncate (* 10 x)))
    117                              (covariance S samples (mean S samples)))) )
     104  (test "sample-covariance"
     105        (f64vector-map (lambda (x) (truncate (* 10 x))) (upper S S sigma))
     106        (f64vector-map (lambda (x) (truncate (* 10 x)))
     107                       (covariance S samples (mean S samples))))
     108  )
    118109
    119 
    120 (test::styler-set! multi-normal-pdf-test test::output-style-compact)
    121 
    122 (run-test "probability distribution test")
    123110
    124111
  • release/4/probdist/trunk/tests/test3.scm

    r6151 r14646  
    1414;;
    1515
    16 (require-extension srfi-1)
    17 (require-extension srfi-4)
    18 (require-extension srfi-4-utils)
    19 (require-extension matrix-utils)
    20 (require-extension blas)
    21 (require-extension normal-pdf)
    22 (require-extension sampled-pdf)
    23 (require-extension random-mtzig)
    24 (require-extension rb-tree)
    25 (require-extension testbase)
    26 (require-extension testbase-output-compact)
     16
     17(use srfi-1 srfi-4 srfi-4-utils matrix-utils blas random-mtzig rb-tree
     18     normal-pdf sampled-pdf  test)
    2719
    2820;; Number of samples
     
    6658             (lambda (i j) (fx= i j)) (list)))))
    6759
     60;; Given a vector of values in the range [0..1], scale all values in
     61   ;; the vector to the specified range.
     62(define (f64vector-urange x lo hi)
     63  (if (f64vector< lo hi)
     64      (let ((d (f64vector-sub hi lo)))
     65        (f64vector-sum lo (f64vector-mul d x)))))
     66
     67(define rng      (random-mtzig:init))
     68
     69;;  distribution parameters
     70(define mu     (f64vector-urange (random-mtzig:f64vector-randu! S rng)
     71                                 (make-f64vector S -5) (make-f64vector S 5)))
     72
     73(define sigma  (let* ((r  (random-mtzig:f64vector-randn! (* S S) rng))
     74                      (x  (f64vector-urange r (make-f64vector (* S S) 0) (make-f64vector (* S S) 5))))
     75                 (blas:dgemm! blas:RowMajor blas:NoTrans blas:Trans S S S
     76                              1.0 x x 0.0 (make-f64vector (* S S)))))
     77
     78(define gpdf  (begin
     79                (make-normal-pdf S mu sigma)))
     80
     81(define var    (diag S S sigma))
     82(define sdev   (f64vector-map sqrt var))
     83
     84(define samples+weights
     85  (let ((hi (f64vector-scale 3 sdev))
     86        (lo (f64vector-scale -3 sdev)))
     87    (let loop ((i 0) (samples (list)))
     88      (if (< i N)
     89          (let ((x (f64vector-urange (random-mtzig:f64vector-randu! S rng) lo hi)))
     90            (let* ((sample   (normal-pdf:sample gpdf  x))
     91                   (density  (normal-pdf:density gpdf sample)))
     92              (loop (fx+ 1 i) (cons (cons density sample) samples))))
     93          samples))))
     94
     95(define spdf (make-sampled-pdf S make-rb-tree samples+weights)) 
     96
     97(print "given mean: " mu)
     98(print "given covariance: " sigma)
     99(print "sampled mean: " (sampled-pdf:expectation spdf))
     100(print "sampled covariance: " (sampled-pdf:covariance spdf))
    68101
    69102     
    70 ;(define-test multi-normal-pdf-test "multivariate sampled PDF test"
    71 ;  (initial
     103(test-group  "multivariate sampled PDF test"
    72104
    73    ;; Given a vector of values in the range [0..1], scale all values in
    74    ;; the vector to the specified range.
    75    (define (f64vector-urange x lo hi)
    76      (if (f64vector< lo hi)
    77          (let ((d (f64vector-sub hi lo)))
    78            (f64vector-sum lo (f64vector-mul d x)))))
    79    
    80    (define rng      (random-mtzig:init))
    81 
    82    ;;  distribution parameters
    83    (define mu     (f64vector-urange (random-mtzig:f64vector-randu! S rng)
    84                                     (make-f64vector S -5) (make-f64vector S 5)))
    85 
    86    (define sigma  (let* ((r  (random-mtzig:f64vector-randn! (* S S) rng))
    87                          (x  (f64vector-urange r (make-f64vector (* S S) 0) (make-f64vector (* S S) 5))))
    88                     (blas:dgemm! blas:RowMajor blas:NoTrans blas:Trans S S S
    89                                  1.0 x x 0.0 (make-f64vector (* S S)))))
    90 
    91    (define gpdf  (begin
    92                    (print "mu = " mu)
    93                    (print "sigma = " sigma)
    94                    (make-normal-pdf S mu sigma)))
    95        
    96    (define var    (diag S S sigma))
    97    (define sdev   (f64vector-map sqrt var))
    98 
    99    (define samples+weights
    100      (let ((hi (f64vector-scale 3 sdev))
    101            (lo (f64vector-scale -3 sdev)))
    102        (let loop ((i 0) (samples (list)))
    103          (if (< i N)
    104              (let ((x (f64vector-urange (random-mtzig:f64vector-randu! S rng) lo hi)))
    105                (let* ((sample   (normal-pdf:sample gpdf  x))
    106                       (density  (normal-pdf:density gpdf sample)))
    107                  (print "sample = " sample)
    108                  (print "density = " density)
    109                  (loop (fx+ 1 i) (cons (cons density sample) samples))))
    110              samples))))
    111      
    112    (print "samples = " samples+weights)
    113 
    114    (define spdf (make-sampled-pdf S make-rb-tree samples+weights)) 
    115 
    116    (print "given mean: " mu)
    117    (print "given covariance: " sigma)
    118    (print "sampled mean: " (sampled-pdf:expectation spdf))
    119    (print "sampled covariance: " (sampled-pdf:covariance spdf))
    120 
    121 ;;   ;; test for mean and covariance equality to one decimal place
    122 ;;   (test/equal 'sample-mean
    123 ;;             (f64vector-map (lambda (x) (truncate (* 10 x))) (normal-pdf:expectation gpdf))
    124 ;;             (f64vector-map (lambda (x) (truncate (* 10 x))) (sampled-pdf:expectation spdf)))
     105             ;; test for mean and covariance equality to one decimal place
     106             (test "sample-mean"
     107                   (f64vector-map (lambda (x) (truncate (* 10 x))) mu)
     108                   (f64vector-map (lambda (x) (truncate (* 10 x))) (sampled-pdf:expectation spdf)))
    125109
    126110
    127 ;;   (test/equal 'sample-covariance
    128 ;;            (f64vector-map (lambda (x) (truncate (* 10 x))) (normal-pdf:covariance gpdf))
    129 ;;            (f64vector-map (lambda (x) (truncate (* 10 x))) (sampled-pdf:expectation spdf))) )
     111             (test "sample-covariance"
     112                   (f64vector-map (lambda (x) (truncate (* 10 x))) sigma)
     113                   (f64vector-map (lambda (x) (truncate (* 10 x))) (sampled-pdf:covariance spdf))) )
    130114
    131 ;; (test::styler-set! multi-normal-pdf-test test::output-style-compact)
    132115
    133 ;; (run-test "probability distribution test")
     116
    134117
    135118
Note: See TracChangeset for help on using the changeset viewer.