diff --git a/projects/sceneContext/obsModelTest.m b/projects/sceneContext/obsModelTest.m
new file mode 100644
index 000000000..3329a15a1
--- /dev/null
+++ b/projects/sceneContext/obsModelTest.m
@@ -0,0 +1,82 @@
+
+
+%% Check the reasonableness of the local observation model 
+
+% Data
+loadData('sceneContextSUN09', 'ismatfile', false)
+load('SUN09data')
+
+
+train = data.train;
+test = data.test;
+objectnames = data.names;
+
+[Ntrain, Nobjects] = size(train.presence);
+[Ntest, Nobjects2] = size(test.presence);
+
+obstypes = {'gauss', 'quantize'};
+
+for oo=1:numel(obstypes)
+  obstype = obstypes{oo};
+
+labels = train.presence;
+scores = train.detect_maxprob;
+%[quantizedScores, discretizeParams] = discretizePMTK(scores, 10);
+[obsmodel] = obsModelFit(labels, scores, obstype);
+
+
+% we plot the distribution of scores for 2 classes
+
+for c=[1 110]
+  
+  % Empirical distributon
+  scores = train.detect_maxprob;
+  ndx=(train.presence(:,c)==1);
+  figure;
+  subplot(2,2,1)
+  [counts, bins]=hist(scores(ndx,c));
+  binstr =cellfun(@(b) sprintf('%2.1f', b), num2cell(bins), 'uniformoutput', false);
+  bar(counts); set(gca, 'xticklabel', binstr)
+  title(sprintf('%s present, m %5.3f, v %5.3f', ...
+    objectnames{c}, mean(scores(ndx,c)),var(scores(ndx,c))));
+  
+  subplot(2,2,2)
+  [counts, bins] = hist(scores(~ndx,c));
+  binstr =cellfun(@(b) sprintf('%2.1f', b), num2cell(bins), 'uniformoutput', false);
+  bar(counts); set(gca, 'xticklabel', binstr)
+  title(sprintf('%s absent, m %5.3f, v %5.3f', ...
+    objectnames{c}, mean(scores(~ndx,c)), var(scores(~ndx,c))));
+  
+  % Model distribution
+  
+  switch obsmodel.obsType
+    case 'gauss'
+      xmin = min(scores(:,c));
+      xmax = max(scores(:,c));
+      xvals = linspace(xmin, xmax, 100);
+      mu = squeeze(obsmodel.mu(1,:,c));
+      Sigma = permute(obsmodel.Sigma(:,:,:,c), [3 4 1 2]);
+      p = gaussProb(xvals, mu(2), Sigma(2));
+      subplot(2,2,3)
+      plot(xvals, p, 'b-');
+      title(sprintf('model for %s presence', objectnames{c}))
+      subplot(2,2,4)
+      p = gaussProb(xvals, mu(1), Sigma(1));
+      plot(xvals, p, 'r:');
+      title(sprintf('model for %s absence', objectnames{c}))
+    case 'quantize'
+      % CPT(label, feature, node)
+      subplot(2,2,3)
+      bar(squeeze(obsmodel.CPT(2,:,c)))
+      title(sprintf('model for %s presence', objectnames{c}))
+      bins = obsmodel.discretizeParams.bins{c};
+      binstr =cellfun(@(b) sprintf('%2.1f', b), num2cell(bins), 'uniformoutput', false);
+      set(gca,'xticklabel', binstr)
+      subplot(2,2,4)
+      bar(squeeze(obsmodel.CPT(1,:,c)))
+      title(sprintf('model for %s absence', objectnames{c}))
+      set(gca,'xticklabel',binstr)
+  end
+end
+
+end
diff --git a/projects/sceneContext/tagContextDemo3.m b/projects/sceneContext/tagContextDemo3.m
index 155f19d97..66814cb66 100644
--- a/projects/sceneContext/tagContextDemo3.m
+++ b/projects/sceneContext/tagContextDemo3.m
@@ -11,7 +11,7 @@
 
 
 %% Models/ methods
-methodNames  = { 'dag' };
+methodNames  = { 'tree' };
 %methodNames  = { 'dag', 'mix10', 'mix15', 'tree' };
 
 
@@ -19,15 +19,16 @@
 % model = fn(truth(N, D), features(N, D, :))
 % where truth(n,d) in {0,1}
 
+%{
 fitMethods = {
   @(labels, features) dgmFit(labels)
   };
+%}
+
 
-%{
 fitMethods = {
-  @(labels, features) noisyMixModelFit(labels, [], 1)
+  @(labels, features) treegmFit(labels)
   };
-%}
 
 
 %{
@@ -49,7 +50,6 @@
 infMethods = {
   @(model, features, softev) argout(2, @noisyMixModelInferNodes, model, [], softev)
   };
-%}
 
 
 infMethods = {
@@ -58,13 +58,18 @@
   @(model, features, softev) argout(2, @noisyMixModelInferNodes, model, [], softev), ...
   @(model, features, softev) argout(2, @treegmInferNodes, model, [], softev)
   };
+%}
+
+
+infMethods = {
+  @(model, features, softev) argout(2, @treegmInferNodes, model, [], softev)
+  };
 
 
 %{
 logprobMethods = {
   @(model, X) mixModelLogprob(model.mixmodel, X)
   };
-%}
   
 
 logprobMethods = {
@@ -74,6 +79,11 @@
   @(model, X) treegmLogprob(model, X)
   };
 
+  %}
+  
+logprobMethods = {
+  @(model, X) treegmLogprob(model, X)
+  };
 
 
 %% CV
@@ -88,7 +98,7 @@
 detect_maxprob = [data.train.detect_maxprob; data.test.detect_maxprob];
 N = size(presence, 1);
 assert(N==Ntrain+Ntest);
-Nfolds = 1;
+Nfolds = 3;
 if Nfolds == 1
   % use original train/ test split
   trainfolds{1} = 1:Ntrain;
@@ -102,83 +112,20 @@
   train.detect_maxprob = detect_maxprob(trainfolds{fold}, :);
   test.presence = presence(testfolds{fold}, :);
   test.detect_maxprob = detect_maxprob(testfolds{fold}, :);
-  
- 
   [Ntrain, Nobjects] = size(train.presence);
   [Ntest, Nobjects2] = size(test.presence);
   
     %% Train  p(scores | labels)
-  %obstype = 'localev';
-  obstype = 'gauss';
-  %obstype = 'quantize';
   
   labels = train.presence;
   scores = train.detect_maxprob;
   %[quantizedScores, discretizeParams] = discretizePMTK(scores, 10);
-  [obsmodel] = obsModelFit(labels, scores, obstype);
-  
- 
-  %% Check the reasonableness of the local observation model for class c
-  % note that p(score|label) is same for all models
-
-  %{
-for c=[1 110]
-
-% Empirical distributon
-scores = train.detect_maxprob;
-ndx=(train.presence(:,c)==1);
-figure;
-subplot(2,2,1)
-[counts, bins]=hist(scores(ndx,c));
-binstr =cellfun(@(b) sprintf('%2.1f', b), num2cell(bins), 'uniformoutput', false);
-bar(counts); set(gca, 'xticklabel', binstr)
-title(sprintf('%s present, m %5.3f, v %5.3f', ...
-  objectnames{c}, mean(scores(ndx,c)),var(scores(ndx,c))));
-
-subplot(2,2,2)
-[counts, bins] = hist(scores(~ndx,c));
-binstr =cellfun(@(b) sprintf('%2.1f', b), num2cell(bins), 'uniformoutput', false);
-bar(counts); set(gca, 'xticklabel', binstr)
-title(sprintf('%s absent, m %5.3f, v %5.3f', ...
-  objectnames{c}, mean(scores(~ndx,c)), var(scores(~ndx,c))));
-
-% Model distribution
-
-switch obsmodel.obstype
-  case 'gauss'
-    xmin = min(scores(:,c));
-    xmax = max(scores(:,c));
-    xvals = linspace(xmin, xmax, 100);
-    mu = squeeze(obsmodel.mu(1,:,c));
-    Sigma = permute(obsmodel.Sigma(:,:,:,c), [3 4 1 2]);
-    p = gaussProb(xvals, mu(2), Sigma(2));
-    subplot(2,2,3)
-    plot(xvals, p, 'b-');
-    title(sprintf('model for %s presence', objectnames{c}))
-    subplot(2,2,4)
-    p = gaussProb(xvals, mu(1), Sigma(1));
-    plot(xvals, p, 'r:');
-    title(sprintf('model for %s absence', objectnames{c}))
-  case 'quantize'
-    % CPT(label, feature, node)
-    subplot(2,2,3)
-    bar(squeeze(obsmodel.CPT(2,:,c)))
-    title(sprintf('model for %s presence', objectnames{c}))
-    bins = obsmodel.discretizeParams.bins{c};
-    binstr =cellfun(@(b) sprintf('%2.1f', b), num2cell(bins), 'uniformoutput', false);
-    set(gca,'xticklabel', binstr)
-    subplot(2,2,4)
-    bar(squeeze(obsmodel.CPT(1,:,c)))
-    title(sprintf('model for %s absence', objectnames{c}))
-    set(gca,'xticklabel',binstr)
-end
-end
-  
+  obstype = 'gauss';
+  [obsmode] = obsModelFit(labels, scores, obstype);
   
-  %}
   
   
-  %% Training p(labels, scores)
+  %% Training p(labels)
   
   Nmethods = numel(methodNames);
   models = cell(1, Nmethods);
@@ -194,50 +141,7 @@
     models{m} = fitMethods{m}(labels, scores);
   end
   
-  keyboard
-  
-  %{
-  if isfield(models{1}, 'mixmodel') && models{1}.mixmodel.nmix==1
-    assert(approxeq(priorProb, models{1}.mixmodel.cpd.T(1,2,:)))
-    priorProb(1:5)
-    squeeze(models{1}.mixmodel.cpd.T(1,2,1:5))
-  end
-  %}
-  
-  
-  % Fit a depnet to the labels
-model = depnetFit(data.train.presence, 'nodeNames', data.names)
-% folder =  fileparts(which(mfilename())
-folder = '/home/kpmurphy/Dropbox/figures';
-% for some reason, the directed graph is much more readable
-graphviz(model.G, 'labels', model.nodeNames, 'directed', 1, ...
-  'filename', fullfile(folder, 'SUN09depnet'));
-  
 
-  %{
-% Visualize tree
-folder =  fileparts(which(mfilename())
-folder = '/home/kpmurphy/Dropbox/figures';
-% for some reason, the directed graph is much more readable
-graphviz(model.edge_weights, 'labels', train.names, 'directed', 1, ...
-  'filename', fullfile(folder, 'SUN09treeNeg'));
-  %}
-  
-  %{
-% Visualize mix model
-model = models{2};
-K = model.mixmodel.nmix;
-[nr,nc] = nsubplots(K);
-figure;
-for k=1:K
-  T = squeeze(model.mixmodel.cpd.T(k,2,:));
-  subplot(nr, nc, k)
-  bar(T);
-  [probs, perm] = sort(T, 'descend');
-  memberNames = sprintf('%s,', train.names{perm(1:5)})
-  title(sprintf('%5.3f, %s', model.mixmodel.mixWeight(k), memberNames))
-end
-  %}
   
   %% Probability of labels
   % See if the models help with p(y(1:T))
@@ -279,20 +183,21 @@
     for n=1:Ntest
       frame = n;
       if (n==1) || (mod(n,500)==0), fprintf('testing image %d of %d\n', n, Ntest); end
-      
-    %{
-      % needs database of images
-    img = imread(fullfile(HOMEIMAGES, test.folder{frame}, test.filename{frame}));
-    figure(1); clf; image(img)
-    trueObjects = sprintf('%s,', test.names{find(test.presence(frame,:))});
-    title(trueObjects)
-     %}
-      
       softev = softevBatch(:,:,n); % Nstates * Nnodes * 1
       [presence_indep(n,:)] = features(n, :);
       %[presence_indep(n,:)] = softev(2, :);
       bel = infMethods{m}(models{m}, features(n,:), softev);
       presence_model(n,:,m) = bel(2,:);
+      
+         
+      % visualize predictions - % needs database of images
+    if fold==1 && (n <= 3)
+      img = imread(fullfile(HOMEIMAGES, test.folder{frame}, test.filename{frame}));
+      figure(1); clf; image(img)
+      trueObjects = sprintf('%s,', test.names{find(test.presence(frame,:))});
+      title(trueObjects)
+    end
+      
     end
   end
   
diff --git a/projects/sceneContext/tagContextDemo4.m b/projects/sceneContext/tagContextDemo4.m
new file mode 100644
index 000000000..38601a634
--- /dev/null
+++ b/projects/sceneContext/tagContextDemo4.m
@@ -0,0 +1,268 @@
+
+% This is like tagContextDemo except we separate out the observation
+% model, p(features|labels), which is shared by different priors p(labels)
+% We also do cross validation 
+
+%% Data
+loadData('sceneContextSUN09', 'ismatfile', false)
+load('SUN09data')
+objectnames  = data.names;
+
+
+
+%% Models/ methods
+
+methods = [];
+m = 0;
+
+m = m + 1;
+methods(m).modelname = 'tree';
+methods(m).obstype = 'gauss';
+methods(m).fitFn = @(labels, features) dgmFit(labels);
+methods(m).infFn = @(model, features, softev) argout(2, @treegmInferNodes, model, [], softev);
+methods(m).logprobFn = @(model, labels) treegmLogprob(model, labels);
+
+m = m + 1;
+methods(m).modelname = 'mix1';
+methods(m).obstype = 'detector';
+methods(m).fitFn = @(labels, features) noisyMixModelFit(labels, [], 1);
+methods(m).infFn = @(model, features, softev) argout(2, @noisyMixModelInferNodes, model, [], softev);
+methods(m).logprobFn = @(model, labels) mixModelLogprob(model, labels);
+
+
+Nmethods = numel(methods);
+
+
+%[logZ, nodeBel] = treegmInferNodes(treeModel, localFeatures, softev);
+%[pZ, pX] = noisyMixModelInferNodes(mixModel{ki}, localFeatures, softev);
+
+
+%% CV
+
+setSeed(0);
+[Ntrain Nobjects] = size(data.train.presence);
+objectnames = data.names;
+Ntest = size(data.test.presence,1);
+presence = [data.train.presence; data.test.presence];
+detect_maxprob = [data.train.detect_maxprob; data.test.detect_maxprob];
+N = size(presence, 1);
+assert(N==Ntrain+Ntest);
+Nfolds = 3;
+if Nfolds == 1
+  % use original train/ test split
+  trainfolds{1} = 1:Ntrain;
+  testfolds{1} = (Ntrain+1):(Ntrain+Ntest);
+else
+  [trainfolds, testfolds] = Kfold(N, Nfolds);
+end
+for fold=1:Nfolds
+  fprintf('fold %d of %d\n', fold, Nfolds);
+  train.presence = presence(trainfolds{fold}, :);
+  train.detect_maxprob = detect_maxprob(trainfolds{fold}, :);
+  test.presence = presence(testfolds{fold}, :);
+  test.detect_maxprob = detect_maxprob(testfolds{fold}, :);
+  [Ntrain, Nobjects] = size(train.presence);
+  [Ntest, Nobjects2] = size(test.presence);
+  
+    %% Train  p(scores | labels)
+  
+  labels = train.presence;
+  scores = train.detect_maxprob;
+  %[quantizedScores, discretizeParams] = discretizePMTK(scores, 10);
+  [obsmodelGauss] = obsModelFit(labels, scores, 'gauss');
+  
+  
+  
+  %% Training p(labels)
+
+  models = cell(1, Nmethods);
+  % indep model
+  Npresent = sum(train.presence, 1);
+  priorProb = Npresent/Ntrain;
+
+  
+  % Train up models
+  for m=1:Nmethods
+    methodNames{m} = sprintf('%s-%s', methods(m).modelname, methods(m).obstype);
+    fprintf('fitting %s\n', methodNames{m});
+    models{m} = methods(m).fitFn(labels, scores);
+  end
+  
+
+  
+  %% Probability of labels
+  % See if the models help with p(y(1:T))
+  ll_indep = zeros(1, Ntest); %#ok
+  ll_model = zeros(Ntest, Nmethods);
+  labels = test.presence+1; % 1,2
+  
+  %logPrior = [log(1-priorProb+eps); log(priorProb+eps)];
+  logPrior = [log(1-priorProb); log(priorProb)];
+  ll = zeros(Ntest, Nobjects);
+  for j=1:Nobjects
+    ll(:,j) = logPrior(labels(:, j), j);
+  end
+  ll_indep = sum(ll,2);
+  
+  for m=1:Nmethods
+    ll_model(:, m) = methods{m}.logprobFn(models{m}, labels);
+  end
+  
+  ll = [sum(ll_indep) sum(ll_model,1)];
+  loglik(fold, :) = ll;
+  
+ 
+
+
+  %% Inference
+  presence_indep = zeros(Ntest, Nobjects);
+  presence_model = zeros(Ntest, Nobjects, Nmethods);
+  
+  features = test.detect_maxprob; %Ncases*Nnodes*Ndims
+  % as a speedup, we compute soft evidence from features in batch form
+  softevBatchGauss = obsModelEval(obsmodelGauss, features);
+  %[quantizedScores, discretizeParams] = discretizePMTK(scores, 10);
+   
+  for m=1:Nmethods
+    fprintf('running method %s\n', methodNames{m});
+    
+    for n=1:Ntest
+      frame = n;
+      if (n==1) || (mod(n,500)==0), fprintf('testing image %d of %d\n', n, Ntest); end
+      [presence_indep(n,:)] = test.detect_maxprob(n, :);
+      
+      switch method(m).obstype
+        case 'detector'
+          bel = infMethods{m}(models{m}, [], test.detect_maxprob(n, :));
+        case 'gauss'
+          bel = infMethods{m}(models{m}, [], softevBatchGauss(:, :, n));
+      end
+      
+      probPresent = bel(2,:);
+      presence_model(n,:,m) = probPresent;
+      
+         
+      % visualize some predictions -  needs database of images
+    if 0 % fold==1 && (n <= 3)
+      srcfolder = '';
+      destfolder = '';
+      img = imread(fullfile(HOMEIMAGES, srcfolder{frame}, test.filename{frame}));
+      figure(1); clf; image(img)
+      trueObjects = sprintf('%s,', objectnames{find(test.presence(frame,:))});
+      title(trueObjects)
+      fname = fullfile(destfolder, sprintf('testimg%d-truth.png', n));
+      print(gcf, '-dpng', fname);
+      
+      figure(2); clf; image(img)
+      predObjects = sprintf('%s,', objectsnames{find(probPresent > 0.5));
+      title(predObjects)
+      fname = fullfile(destfolder, sprintf('testimg%d-%s.png', n, methodNames{m}));
+      print(gcf, '-dpng', fname);
+    end
+      
+    end % for n
+  end % for m
+  
+  %% Performance evaluation
+  evalPerf =   @(confidence, truth) rocPMTK(confidence, truth);
+  perfStr = 'aROC';
+  
+  % If the object is absent in a given fold, we may get NaN for
+  % the performance. We want to exclude these from the evaluation.
+  score_indep = nan(1, Nobjects);
+  score_models = nan(Nobjects, Nmethods);
+  absent = zeros(1, Nobjects);
+  for c=1:Nobjects
+    absent(c) = all(test.presence(:,c)==0);
+    if absent(c), continue; end
+    [score_indep(c)] = evalPerf(presence_indep(:,c), test.presence(:,c));
+    for m=1:Nmethods
+      score_models(c,m) = evalPerf(presence_model(:,c,m), test.presence(:,c));
+    end
+  end
+  if ~isempty(absent)
+    fprintf('warning: in fold %d of %d, %s are absent from test\n', ...
+      fold, Nfolds, sprintf('%s,', objectnames{find(absent)}));
+  
+  mean_perf_indep(fold) = mean(score_indep(~absent));
+  for m=1:Nmethods
+    mean_perf_models(fold,m) = mean(score_models(~absent,m));
+  end
+  
+  
+end % fold
+
+
+
+%% Plotting
+
+[styles, colors, symbols, plotstr] =  plotColors();
+
+
+%{
+    % plot improvement over baseline for each method as separate figs
+    for m=1:Nmethods
+      figure;
+      [delta, perm] = sort(score_models(:,m) - score_indep(:), 'descend');
+      bar(delta)
+      str = sprintf('mean of %s using indep %5.3f, using %s  %5.3f', ...
+        perfStr, mean_perf_indep(fold), methodNames{m},  mean_perf_models(fold,m));
+      disp(str)
+      title(str)
+      xlabel('category')
+      ylabel(sprintf('improvement in %s over baseline', perfStr))
+    end
+  end
+%}
+
+% Plot some ROC curves for some classes on a single fold
+if 1
+  for c=1:3
+    absent(c) = all(test.presence(:,c)==0);
+    if absent(c), 
+      error(sprintf('%s is absent in this test fold', objectnames{c}))
+    end
+    figure;
+    [aROC, falseAlarmRate, detectionRate] = figROC(presence_indep(:,c), test.presence(:,c), colors(1));
+    for m=1:Nmethods
+      [aROC, falseAlarmRate, detectionRate] = figROC(presence_model(:,c,m), test.presence(:,c), colors(m+1));
+    end
+    legendstr = {'indep', methodNames{:}};
+    legend(legendstr)
+    title(objectnames{c})
+  end
+end
+
+
+% Mean AUC
+figure;
+if Nfolds==1
+  %bar([mean_perf_indep mean_perf_models])
+  plot([mean_perf_indep mean_perf_models], 'x', 'markersize', 12, 'linewidth', 2)
+  axis_pct
+else
+  perf = [mean_perf_indep(:)   mean_perf_models]; % folds * methods
+  boxplot(perf)
+end
+legendstr = {'indep', methodNames{:}};
+set(gca, 'xtick', 1:numel(legendstr))
+set(gca, 'xticklabel', legendstr)
+title(sprintf('mean %s', perfStr))
+
+
+% NLL on labels
+figure;
+if Nfolds==1
+  plot(-loglik, 'x', 'markersize', 12, 'linewidth', 2)
+  legendstr = {'indep', methodNames{:}};
+  set(gca, 'xtick', 1:numel(legendstr))
+  set(gca, 'xticklabel', legendstr)
+  axis_pct
+else
+  boxplot(-loglik)
+end
+legendstr = {'indep', methodNames{:}};
+set(gca, 'xtick', 1:numel(legendstr))
+set(gca, 'xticklabel', legendstr)
+title(sprintf('negloglik'))
+
diff --git a/projects/sceneContext/tagModels.m b/projects/sceneContext/tagModels.m
index 8638cda7f..3a3584915 100644
--- a/projects/sceneContext/tagModels.m
+++ b/projects/sceneContext/tagModels.m
@@ -87,6 +87,15 @@
   
   %{
 % Visualize mix model
+  
+  
+  if isfield(models{1}, 'mixmodel') && models{1}.mixmodel.nmix==1
+    assert(approxeq(priorProb, models{1}.mixmodel.cpd.T(1,2,:)))
+    priorProb(1:5)
+    squeeze(models{1}.mixmodel.cpd.T(1,2,1:5))
+  end
+ 
+  
 model = models{2};
 K = model.mixmodel.nmix;
 [nr,nc] = nsubplots(K);