	clear; clc; close all;
	

	rng(1);

	% set the distributions 
	d = 2;
	mu1 = -ones(d, 1);
	Sigma1 = eye(d);
	mu2 = ones(d, 1);
	Sigma2 = eye(d);
	n1 = 20;
	n2 = 20;
	Nte = 100;
	truep = 0.3;

	% generate the training datasets
	x1 = mvnrnd(mu1, Sigma1, n1)';
	x2 = mvnrnd(mu2, Sigma2, n2)';
	
	% generate the test dataset
	ntea = binornd(Nte, truep);
	nteb = Nte - ntea;
	
	xa = mvnrnd(mu1, Sigma1, ntea)';
	xb = mvnrnd(mu2, Sigma2, nteb)';
	xte = [xa xb];
	xte = xte(:, randperm(size(xte, 2)));

	% plot the original functions
	hfig = figure; hold on;
	hs1 = subplot(1, 2, 1);
	scatter(x1(1, :), x1(2, :), 'r', 'filled'); hold on;
	xlabel('x_1');
	ylabel('x_2');
	scatter(x2(1, :), x2(2, :), 'g', 'filled');
	legend('Class 1', 'Class 2');
	title('Labeled samples');
	v1 = axis();
	hs2 = subplot(1, 2, 2);
	scatter(xte(1, :), xte(2, :), 'k', 'filled');
	v2 = axis();
	xlabel('x_1');
	ylabel('x_2');
	title('Samples from mixture x_{te}');

	v(1) = min(v1(1), v2(1));
	v(2) = max(v1(2), v2(2));
	v(3) = min(v1(3), v2(3));
	v(4) = max(v1(4), v2(4));

	axis(v);
	subplot(hs1); axis(v);

	% estimate the class prior
	[xi_best, xi_list, LSDE] =  LSDDPriorEstMedian(xte, {x1, x2});

	% plot the result
	hfig1 = figure; hold on;
	plot(xi_list, LSDE, 'k', 'LineWidth', 2);
	[~, idx] = min((xi_list - xi_best).^2);
	scatter(xi_best, LSDE(idx), 'filled');
	xlabel('Candidate class prior (\theta)');
	ylabel('L_2(\theta)');
	title(sprintf('Est. theta %g, true theta %g', xi_best, truep));
	
	% 
	