(Arne Babenhauserheide)
2016-11-07: also plot the development of the ensemble also plot the development of the ensemble
diff --git a/examples/ensemble-estimation.w b/examples/ensemble-estimation.w
--- a/examples/ensemble-estimation.w
+++ b/examples/ensemble-estimation.w
@@ -1,5 +1,6 @@
#!/usr/bin/env sh
# -*- wisp -*-
+guile -L $(dirname $(dirname $(realpath "$0"))) -c '(import (wisp-scheme) (language wisp spec))'
exec guile -L $(dirname $(dirname $(realpath "$0"))) --language=wisp -e '(@@ (examples ensemble-estimation) main)' -s "$0" "$@"
; !#
@@ -78,19 +79,24 @@ define* : write-multiple . x
;; Start with the simple case: One variable and independent observations (R diagonal)
;; First define a truth
-define x^seed '(0.5 0.6 2 0.1) ; 0.7 0.9 0.8 0.4)
+define x^seed '(0.5 2 0.6 0.1) ; 0.7 0.9 0.8 0.4)
;; The size is the length of the seed, squared, each multiplied by each
define x^true : append-ec (: i (length x^seed)) : list-ec (: j x^seed) : * j : list-ref x^seed i
;; And add an initial guess of the parameters
define x^b : append-ec (: i (length x^seed)) '(1 1 1 1) ; 1 1 1 1) ; initial guess
-define P : make-covariance-matrix-from-standard-deviations : append-ec (: i (length x^seed)) '(0.5 0.1 0.3 0.1) ; 0.2 0.2 0.2 0.2)
+;; set x^b as x^true to test losing uncertainty
+define x^b x^true
+define x^seed-std : append-ec (: i (length x^seed)) '(0.5 0.1 0.3 0.1) ; 0.2 0.2 0.2 0.2)
+define P : make-covariance-matrix-from-standard-deviations x^seed-std
;; Then generate observations
-define y⁰-num 1000
+define y⁰-num 600 ; careful: N² in memory, 1000 is still OK
define y⁰-pos-max 100
;; At the positions where they are measured. Drawn randomly to avoid
;; giving an undue weight to later values.
-define y⁰-pos : list-ec (: i y⁰-num) : * (random:uniform) y⁰-pos-max
+define y⁰-pos-sorted : list-ec (: i y⁰-num) : exact->inexact : * y⁰-pos-max : / i y⁰-num
+define y⁰-pos-random : list-ec (: i y⁰-num) : * (random:uniform) y⁰-pos-max
+define y⁰-pos y⁰-pos-random
define : H-single-parameter xi xi-pos pos
. "Observation function for a single parameter."
@@ -124,7 +130,7 @@ x are parameters to be optimized, pos is
;; the equivalent of measured observations
define y^true : list-ec (: i y⁰-pos) : H x^true i
;; now we disturb the observations with a fixed standard deviation. This assumes uncorrelated observations.
-define y⁰-std 10
+define y⁰-std 1
define y⁰ : list-ec (: i y^true) : + i : * y⁰-std : random:normal
;; and define the covariance matrix. This assumes uncorrelated observations.
define R : make-covariance-matrix-from-standard-deviations : list-ec (: i y⁰-num) y⁰-std
@@ -135,6 +141,8 @@ define R : make-covariance-matrix-from-s
;; define y⁰ : list-ec (: i y⁰-num) : + y⁰-mean : * y⁰-std : random:normal
+define x^steps '()
+
define : EnSRT H x P y R y-pos N
. "Observation function H, parameters x,
parameter-covariance P, observations y, observation covariance R
@@ -147,11 +155,12 @@ Limitations: y is a single value. R and
observation-variances : list-ec (: i (length y)) : list-ref (list-ref R i) i
observation-positions y-pos
x^b x
- x-deviations
+ x-deviations
list-ec (: i N)
list-ec (: j (length x))
* : random:normal
sqrt : list-ref (list-ref P j) j ; only for diagonal P!
+ set! x^steps : cons x-deviations x^steps
cond
: null? observations-to-process
list x^b x-deviations
@@ -160,7 +169,7 @@ Limitations: y is a single value. R and
; newline
let*
: y_cur : car observations-to-process
- R_cur : car observation-variances
+ R_cur : * 1 : car observation-variances
y-pos_cur : car observation-positions
Hx^b_i
list-ec (: i x-deviations)
@@ -207,36 +216,67 @@ Limitations: y is a single value. R and
. x^a-deviations
+define : x-deviations->y-deviations H x-opt x-deviations y⁰-pos
+ . "Calculate y-deviations for each measurement from the x-deviations.
+
+ return ((y_0'0 y_0'1 ...) ...)
+ "
+ define (~ li i) (list-ref li i)
+ ; for each ensemble member calculate the y-deviations
+ let*
+ : y-opt : map (λ (x) (H x-opt x)) y⁰-pos
+ y-ensemble
+ list-ec (: x-dev x-deviations) ; x-dev has one number (x'_i) per x-parameter
+ let*
+ : ; calculate x-parameters for the ensemble
+ x-opt+dev
+ list-ec (: j (length x-opt))
+ + : ~ x-opt j
+ ~ x-dev j
+ ; calculate all y-values for the ensemble
+ y-opt+dev : map (λ (x) (H x-opt+dev x)) y⁰-pos
+ ; calculate all differences
+ map (λ (x y) (- x y)) y-opt+dev y-opt
+ ; reshape into one ensemble per value
+ list-ec (: ensidx (length (~ y-ensemble 0)))
+ list-ec (: obsidx (length y-ensemble))
+ ~ (~ y-ensemble obsidx) ensidx
+
+define : flatten li
+ append-ec (: i li) i
+
define : main args
let*
- : optimized : EnSRT H x^b P y⁰ R y⁰-pos 40
+ : ensemble-member-count 16
+ optimized : EnSRT H x^b P y⁰ R y⁰-pos ensemble-member-count
x-opt : list-ref optimized 0
x-deviations : list-ref optimized 1
- ; std : sqrt : * {1 / {(length x-deviations) - 1}} : sum-ec (: i x-deviations) : expt i 2
- format #t "x⁰: ~A ± ~A\nx: ~A ± ~A\nx^t: ~A\nx-t/σ:~A\ny̅: ~A ± ~A\ny̅⁰: ~A ± ~A\ny̅^t: ~A\nnoise:~A\n"
+ x-std
+ list-ec (: i (length x-opt))
+ apply standard-deviation-from-deviations : list-ec (: j x-deviations) : list-ref j i
+ y-deviations : x-deviations->y-deviations H x-opt x-deviations y⁰-pos
+ x^b-deviations-approx
+ list-ec (: i ensemble-member-count)
+ list-ec (: j (length x^b))
+ * : random:normal
+ sqrt : list-ref (list-ref P j) j ; only for diagonal P!
+ y^b-deviations : x-deviations->y-deviations H x^b x^b-deviations-approx y⁰-pos
+ y-std
+ apply standard-deviation-from-deviations
+ flatten y-deviations
+ y-stds : list-ec (: i y-deviations) : apply standard-deviation-from-deviations i
+ y^b-stds : list-ec (: i y^b-deviations) : apply standard-deviation-from-deviations i
+ format #t "x⁰: ~A\n ± ~A\nx: ~A\n ± ~A\nx^t: ~A\nx-t/σ:~A\ny̅: ~A ± ~A\ny̅⁰: ~A ± ~A\ny̅^t: ~A\nnoise:~A\n"
. x^b
list-ec (: i (length x^b)) : list-ref (list-ref P i) i
. x-opt
- list-ec (: i (length x-opt))
- apply standard-deviation-from-deviations : list-ec (: j x-deviations) : list-ref j i
+ . x-std
. x^true
list-ec (: i (length x-opt))
/ : - (list-ref x-opt i) (list-ref x^true i)
apply standard-deviation-from-deviations : list-ec (: j x-deviations) : list-ref j i
- * {1 / (length y⁰)} : apply + : map (lambda (x) (H x-opt x)) y⁰-pos
- apply standard-deviation-from-deviations
- append-ec (: i (length x-deviations))
- let*
- :
- x-opt+dev
- list-ec (: j (length x-opt))
- + : list-ref x-opt j
- list-ref
- list-ref x-deviations i
- . j
- y-opt+dev : map (lambda (x) (H x-opt+dev x)) y⁰-pos
- y-opt : map (lambda (x) (H x-opt x)) y⁰-pos
- map (lambda (x y) (- x y)) y-opt+dev y-opt
+ mean : map (lambda (x) (H x-opt x)) y⁰-pos
+ . y-std
; list-ec (: i (length y-opt))
; - (list-ref y-opt+dev i) (list-ref y-opt i)
; apply standard-deviation-from-deviations : map H x-deviations ; FIXME: This only works for trivial H.
@@ -246,17 +286,28 @@ define : main args
. y⁰-std
; now plot the result
let : : port : open-output-pipe "python"
- format port "import pylab as pl\n"
+ format port "import pylab as pl\nimport matplotlib as mpl\n"
format port "y0 = [float(i) for i in '~A'[1:-1].split(' ')]\n" y⁰
format port "ypos = [float(i) for i in '~A'[1:-1].split(' ')]\n" y⁰-pos
format port "yinit = [float(i) for i in '~A'[1:-1].split(' ')]\n" : list-ec (: i y⁰-pos) : H x^b i
+ format port "yinitstds = [float(i) for i in '~A'[1:-1].split(' ')]\n" y^b-stds
format port "ytrue = [float(i) for i in '~A'[1:-1].split(' ')]\n" : list-ec (: i y⁰-pos) : H x^true i
format port "yopt = [float(i) for i in '~A'[1:-1].split(' ')]\n" : list-ec (: i y⁰-pos) : H x-opt i
- format port "pl.plot(*zip(*sorted(zip(ypos, yinit))), label='prior')\n"
+ format port "yoptstds = [float(i) for i in '~A'[1:-1].split(' ')]\n" y-stds
+ format port "pl.errorbar(*zip(*sorted(zip(ypos, yinit))), yerr=zip(*sorted(zip(ypos, yinitstds)))[1], label='prior')\n"
format port "pl.plot(*zip(*sorted(zip(ypos, ytrue))), label='true')\n"
- format port "pl.plot(*zip(*sorted(zip(ypos, yopt))), label='optimized')\n"
+ format port "pl.errorbar(*zip(*sorted(zip(ypos, yopt))), yerr=zip(*sorted(zip(ypos, yoptstds)))[1], label='optimized')\n"
format port "pl.plot(*zip(*sorted(zip(ypos, y0))), marker='+', linewidth=0, label='measurements')\n"
- format port "pl.legend()\n"
+ list-ec (: step 0 (length x^steps) 10) ; stepsize 10: sample one in 10 steps
+ list-ec (: member (list-ref x^steps (- (length x^steps) step 1))) ; reversed
+ begin
+ format port "paired = pl.get_cmap('Paired')
+cNorm = mpl.colors.Normalize(vmin=~A, vmax=~A)
+scalarMap = mpl.cm.ScalarMappable(norm=cNorm, cmap=paired)\n" 0 (length member)
+ list-ec (: param-idx 0 (length member) 4) ; step = 4
+ ; plot parameter 0
+ format port "pl.plot(~A, ~A, marker='.', color=scalarMap.to_rgba(~A), linewidth=0, label='', alpha=0.3, zorder=-1)\n" (/ step 10) (+ 80 (* 60 (list-ref member param-idx))) param-idx
+ format port "pl.legend(loc='upper right')\n"
format port "pl.xlabel('position [arbitrary units]')\n"
format port "pl.ylabel('value [arbitrary units]')\n"
format port "pl.title('ensemble optimization results')\n"