function [X,Y,Z]=epdf_2dplot(pdf,xlims,ylims, points)
% function evaluates given pdf on support
%  xlims = [x_begin, x_end] --or-- empty
%  ylims = [x_begin, x_end] --or-- empty
%  points = if not given, 100

if nargin<4
    points = [100,100];
end

if length(points) ==1
    points = points*[1 1];
end

if nargin<3 | isempty(ylims)
    p_mean = epdf_mean(pdf);
    p_var  = epdf_variance(pdf);
    ylims = p_mean(2)*[1 1] + [-3 3]*sqrt(p_var(2));
end
if nargin<2 | isempty(xlims)
    p_mean = epdf_mean(pdf);
    p_var  = epdf_variance(pdf);
    xlims = p_mean(1)*[1 1] + [-3 3]*sqrt(p_var(1));
end

steps = [xlims(2)-xlims(1), ylims(2)-ylims(1)];
[X,Y]=meshgrid(xlims(1):steps(1)/(points(1)-1):xlims(2), ...
                ylims(1):steps(2)/(points(2)-1):ylims(2));
data = [reshape(X, 1, prod(points)); reshape(Y,1,prod(points))];
v = epdf_evallog_mat(pdf,data);
Z = reshape(exp(v), points(1), points(2));

if 0
    vol = prod(steps./(points-1));
    nc = sum(sum(Z))*vol
    % renormalize
    p = exp(v)/nc;
    num_mean = vol*data*p
    num_cov = vol*data*(p*ones(1,2).*data') - num_mean * num_mean'
end
contour(X,Y,Z);