changeset 20:982aa1197a7e

Getting there, slowly, sort of, with EM

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Thu Mar 27 11:52:07 2014 +0000
@@ -0,0 +1,4 @@
+syntax: glob
+*~
+*.orig
+

--- a/yeti/em_onecolumn.yeti	Wed Mar 26 18:49:12 2014 +0000
+++ b/yeti/em_onecolumn.yeti	Thu Mar 27 11:52:07 2014 +0000
@@ -5,6 +5,8 @@
 vec = load may.vector;
 mat = load may.matrix;
 
+plot = load may.plot;
+
 inRange ranges instrument note =
     note >= ranges[instrument].lowest and note <= ranges[instrument].highest;
 
@@ -53,31 +55,45 @@
     vec.divide column estimate);
 
 performMaximisation data column q =
-   (pitchTop = vec.fromList
+   (pitches = vec.fromList
        (map do note:
             fold do acc instrument:
                 { w, p, s } = distributionsFor data instrument note;
                 fold do acc bin:
-                    acc + s * (vec.at w bin) * (vec.at column bin);
+                    acc + s * (vec.at w bin) * (vec.at q bin);
                 done acc [0..vec.length column - 1]
             done epsilon data.instruments;
         done [data.lowest .. data.highest]);
-    pitchBottom = vec.sum pitchTop;
+    pitches = vec.divideBy (vec.sum pitches) pitches;
 
-/*    sourceTop = vec.fromList
-       (map do instrument:
-            fold do acc note:
-                { w, p, s } = distributionsFor data instrument note;
-                fold do acc bin:
-                    acc + p * (vec.at w bin) * (vec.at column bin);
-                done acc [0..vec.length column - 1]
-            done epsilon [data.lowest .. data.highest]
-        done data.instruments);
-    sourceBottom = vec.sum sourceTop;
-*/
+    sources = mapIntoHash id
+        do instrument: vec.fromList
+           (map do note:
+                if not inRange data.ranges instrument note then epsilon else
+                    { w, p, s } = distributionsFor data instrument note;
+                    fold do acc bin:
+                        acc + (vec.at w bin) * (vec.at q bin);
+                    done epsilon [0..vec.length column - 1]
+                fi;
+            done [data.lowest .. data.highest])
+        done data.instruments;
+   
+    sourceDenoms = fold do acc instrument:
+        vec.add [acc, sources[instrument]]
+    done (vec.zeros (data.highest - data.lowest + 1)) data.instruments;
+
+        \() (plot.plot [ Caption "Source numerators for piano", Vector (sources["piano-maps-SptkBGCl"]) ]);
+
+        \() (plot.plot [ Caption "Source denominators", Vector sourceDenoms ]);
+    
+    sources = mapIntoHash id
+        do instrument: 
+            vec.divide sources[instrument] sourceDenoms;
+        done (keys sources);
+ 
     data with { 
-        pitches = vec.divideBy pitchBottom pitchTop,
-//        sources = vec.divideBy sourceBottom sourceTop
+        pitches,
+        sources,
     });
 
 {

--- a/yeti/silvet.yeti	Wed Mar 26 18:49:12 2014 +0000
+++ b/yeti/silvet.yeti	Thu Mar 27 11:52:07 2014 +0000
@@ -66,19 +66,30 @@
 
 col = head (drop 35 columns);
 
-\() (plot.plot [ Vector col ]);
+\() (plot.plot [ Caption "Source frequency distribution", Vector col ]);
 
-oneIteration em1data col =
+sourceGrid d =
+    mat.fromColumns (map do k: d.sources[k] done (sort (keys d.sources)));
+
+\() (plot.plot [ Caption "Source distribution beforehand", Grid (sourceGrid em1data)]);
+
+oneIteration em1data col n =
    (q = em1.performExpectation em1data col;
-    \() (plot.plot [ Vector col, Vector q ]);
+//    \() (plot.plot [ Caption "Frequency distribution and output of E-step for iteration \(n)", Vector col, Vector q ]);
     newdata = em1.performMaximisation em1data col q;
-    \() (plot.plot [ Vector (em1data.pitches), Vector (newdata.pitches) ]);
+    if (n % 2 == 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 (sourceGrid newdata) ]);
+    fi;
     newdata);
 
-em1data = oneIteration em1data col;
-em1data = oneIteration em1data col;
-em1data = oneIteration em1data col;
+iterations = 6;
 
+var d = em1data;
+
+for [1..iterations] do i:
+    d := oneIteration d col i;
+done;    
 
 ();