function mixexp1()
  
% Fit a piece-wise linear regression model.
% Here is the model
%
%  X \
%  | |
%  Q |
%  | /
%  Y
%
% where we condition everything on X, Q is a softmax, and Y is a linear Gaussian.
% Q is hidden, X and Y are observed.

rand('seed', 0);
randn('seed', 0);

sigma = 1;
dag = zeros(3,3);
dag(1,2) = 1; dag(1,3) = 1;
dag(2,3) = 1;
node_sizes = [1 2 1];
discrete_nodes = [2];
bnet = mk_bnet(dag, node_sizes, discrete_nodes);
bnet.bnodes{1} = mk_input_node(bnet, 1);
bnet.bnodes{2} = mk_softmax_node(bnet, 2);
bnet.bnodes{2}.max_iter = 50;
bnet.bnodes{2}.verbose = 0;
bnet.bnodes{3} = mk_gaussian_node(bnet, 3, 'full', 0);
bnet.bnodes{3}.cov = sigma*ones(node_sizes(3), node_sizes(3), node_sizes(2));
bnet.bnodes{3}.adjustable_cov = 0;
bnet = compile_params(bnet);

% Check inference

x = 0.1;
ystar = 1;

%x = data(1,1);
%ystar = data(1,2);

vals = [x ystar];
obs_nodes = [1 3];
[bnet, ll] = enter_evidence(bnet, vals, obs_nodes);
Qpost = marginal_nodes(bnet, 2);

% theta(i,:) = regression vector for expert i
% eta(i,:) = softmax (gating) params for expert i
x1 = [1 x];
ns = node_sizes;
eta = [bnet.bnodes{2}.glim.b1; bnet.bnodes{2}.glim.w1]';
% Weight matrix is size Y * X * Q = ns(3) * ns(1) * ns(2), where Y=1
theta = reshape(bnet.bnodes{3}.weights, ns(1), ns(2))';
theta = [bnet.bnodes{3}.mean' theta];
% yhat(i) = E[y | Q=i, x] = prediction of i'th expert
yhat = theta * x1'; 
% gate_prior(i,:) = Pr(Q=i | x)
gate_prior = exp(eta * x1');
gate_prior = gate_prior / sum(gate_prior);
% lik(i) = Pr(y | Q=i, x)
lik = (1/(sqrt(2*pi)*sigma)) * exp(-(0.5/sigma^2) * ((ystar - yhat) .* (ystar - yhat)));
% gate_posterior(i,:) = Pr(Q=i | x, y)
gate_posterior = gate_prior .* lik;
gate_posterior = gate_posterior / sum(gate_posterior);

assert(approxeq(gate_posterior(:), Qpost.T(:)));



% Check learning
global BNTroot
load([BNTroot '/examples/mixexp_data.txt'], '-ascii');
% Just use 1/10th of the data, to speed things up
data = data_mixexp(1:10:1000, :);
%data = data_mixexp;

plot(data(:,1), data(:,2), '.')

if 0
  % stochastic gradient descent method
  numexp = node_sizes(2); % number of experts
  theta = [0 -1; 0 1];
  eta = [0 0; -5 5];
  
  max_iter = 2;
  [theta1, eta1] = mixture_of_experts(numexp, data, max_iter, theta, eta);
  plot_mixexp(theta1, eta1, data);
end


bnet.bnodes{2}.glim.b1 = [0 0];
bnet.bnodes{2}.glim.w1 = [-5 5];
%bnet.bnodes{2}.glim.b1 = [0 -5];
%bnet.bnodes{2}.glim.w1 = [0 5];
bnet.bnodes{3}.weights = reshape([-1 1], [1 1 2]);
bnet.bnodes{3}.mean = [0 0];

bnet = randomize_params(bnet, 1);
bnet = compile_params(bnet);

eta0 = [bnet.bnodes{2}.glim.b1; bnet.bnodes{2}.glim.w1]';
theta0 = reshape(bnet.bnodes{3}.weights, ns(1), ns(2))';
theta0 = [bnet.bnodes{3}.mean' theta0];
figure(1)
plot_mixexp(theta0, eta0, data);

max_iter = 20;
bnet_hat = learn_params(bnet, data, obs_nodes, max_iter);

eta2 = [bnet_hat.bnodes{2}.glim.b1; bnet_hat.bnodes{2}.glim.w1]';
theta2 = reshape(bnet_hat.bnodes{3}.weights, ns(1), ns(2))';
theta2 = [bnet_hat.bnodes{3}.mean' theta2];

figure(2)
plot_mixexp(theta2, eta2, data);


