function [E,samps] = BR_MC2(sigma,m,n,fname)
% BR_MC2 - produce samples for function, as well as mean of samples
% On input:
%     sigma (float): variance for q function (normal)
%     m (int): burnin amount
%     n (int): desired number of samples
%     fname (string): name of function (distribution) to be considered
% On output:
%     E (float): mean of function at samples
%     samps (nx1 vector): samples produced by method
% Call:
%     [E,samps] = BR_MC2(1,1000,10000,'BR_U');
% Author: 
%     T. Henderson
%     UU
%     Winter 2019
%

x = 0.5;
samps = zeros(m,1);
t = 0;
while t<m
    t = t + 1;
    y = normrnd(x,sigma);
    U = rand;
    v1 = feval(fname,x);
    v2 = feval(fname,y);
    v = v2/v1;
    if isnan(v)|isinf(v)
        v = 1;
    end
    alpha = min(1,v);
    if U<=alpha
        samps(t) = y;
        x = y;
    else
        samps(t) = x;
    end
end
samps = zeros(n,1);
t = 0;
while t<n
    t = t + 1;
    y = normrnd(x,sigma);
    U = rand;
    v1 = feval(fname,x);
    v2 = feval(fname,y);
    v = v2/v1;
    if isnan(v)|isinf(v)
        v = 1;
    end
    alpha = min(1,v);
    if U<=alpha
        samps(t) = y;
        x = y;
    else
        samps(t) = x;
    end
end

E = mean(feval(fname,samps));
tch = 0;
