changeset 25:80f02ff5d37a

I think that previous normalisation was wrong!

--- a/yeti/em_onecolumn.yeti	Mon Mar 31 12:46:24 2014 +0100
+++ b/yeti/em_onecolumn.yeti	Mon Mar 31 15:41:15 2014 +0100
@@ -16,20 +16,26 @@
     else v
     fi);
 
+normaliseSources s =
+   (denoms = fold do acc inst: vec.add [acc, (mat.getColumn inst s)] done
+       (vec.zeros (mat.height s)) [0..(mat.width s)-1];
+    mat.fromColumns
+       (map do inst: vec.divide (mat.getColumn inst s) denoms done
+            [0..(mat.width s)-1]));
+
 initialise ranges templates notes =
    (instruments = sort (keys ranges);
     {
         pitches = // z in the original. 1 per note
             normalise (vec.randoms notes),
-        sources = // u in the original. 1 per note-instrument
-            mat.fromRows
-               (map do note:
-                   normalise
-                      (vec.fromList
-                          (map do instrument:
-                               if inRange ranges instrument note then 1 else 0 fi
-                           done instruments))
-                done [0..notes-1]),
+        sources = normaliseSources // u in the original. 1 per note-instrument
+           (mat.fromColumns
+               (map do instrument:
+                  (vec.fromList
+                      (map do note:
+                           if inRange ranges instrument note then 1 else 0 fi
+                       done [0..notes-1]))
+                done instruments)),
         instruments = array instruments,
         instCount = length instruments,
         noteCount = notes,
@@ -69,7 +75,7 @@
             done acc [data.ranges[inst].lowest .. 
                       data.ranges[inst].highest]
         done (vec.consts epsilon (vec.length column)) [0..data.instCount-1];
-    vec.divide column estimate);
+    { estimate, q = vec.divide column estimate });
 
 performMaximisation data column q =
    (column = normalise column;
@@ -88,32 +94,19 @@
         done [0..data.noteCount-1]);
     pitches = vec.divideBy (vec.sum pitches) pitches;
 
-    sources = mat.fromRows
-       (map do note: vec.fromList
-           (map do inst:
-                if not inRange data.ranges inst note then epsilon else
+    sources = mat.fromColumns
+       (map do inst: vec.fromList
+           (map do note:
+               (if not inRange data.ranges inst note then epsilon else
                     { w, p, s } = distributionsFor data inst note;
                     fold do acc bin:
                         acc + s * p * (vec.at w bin) * (vec.at q bin);
                     done epsilon [0..vec.length column - 1]
-                fi;
-            done [0..data.instCount-1])
-        done [0..data.noteCount-1]);
-   
-    sourceDenoms = fold do acc inst:
-        vec.add [acc, (mat.getColumn inst sources)]
-    done (vec.zeros data.noteCount) [0..data.instCount-1];
+                fi);
+            done [0..data.noteCount-1])
+        done [0..data.instCount-1]);
 
-    //!!! shouldn't this normalisation be the other way? (rows, not columns)    
-/*    sources = mat.fromColumns
-       (map do inst: 
-            vec.divide (mat.getColumn inst sources) sourceDenoms;
-        done [0..data.instCount-1]);
-*/
-    sources = mat.fromRows
-       (map do note: 
-            vec.divide (mat.getRow note sources) sourceDenoms;
-        done [0..data.noteCount-1]);
+    sources = normaliseSources sources;
 
     data with { 
         pitches,

--- a/yeti/silvet.yeti	Mon Mar 31 12:46:24 2014 +0100
+++ b/yeti/silvet.yeti	Mon Mar 31 15:41:15 2014 +0100
@@ -70,13 +70,21 @@
 
 \() (plot.plot [ Caption "Source distribution beforehand", Grid em1data.sources]);
 
+normalise v =
+   (s = vec.sum v;
+    if s > 0 then vec.divideBy s v 
+    else v
+    fi);
+
 oneIteration em1data col n =
-   (q = em1.performExpectation em1data col;
+   ({ estimate, q } = em1.performExpectation em1data col;
 //    \() (plot.plot [ Caption "Frequency distribution and output of E-step for iteration \(n)", Vector col, Vector q ]);
     newdata = em1.performMaximisation em1data col q;
     if (n % 6 == 0) then
         \() (plot.plot [ Caption "Pitch distribution before and after M-step update for iteration \(n)", Vector (em1data.pitches), Vector (newdata.pitches) ]);
         \() (plot.plot [ Caption "Source distribution after M-step update for iteration \(n)", Grid newdata.sources ]);
+        \() (plot.plot [ Caption "Q function for E-step iteration \(n)", Vector q ]);
+        \() (plot.plot [ Caption "Estimate from E-step iteration \(n) against original source distribution, and difference between them", Vector estimate, Vector (normalise col), Vector (vec.subtract estimate (normalise col)) ]);
     fi;
     newdata);
 
@@ -123,12 +131,6 @@
 
 println "Instruments: \(map do i: (d.instruments[i]) done instruments)";
 
-normalise v =
-   (s = vec.sum v;
-    if s > 0 then vec.divideBy s v 
-    else v
-    fi);
-
 if not (empty? sounding) then
    p = head sounding;
    i = head instruments;