clc
clear all
close all

% MAIN: computes Bayes Factor of model M1 (Poisson) and model M2
% (geometric) for the n data generated under M1 or M2, using a uniform and
% improper prior for Poisson, and a uniform but proper prior for geometric


%%%%%%% SET PARAMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 100; % number of observations
model = 1; % 1= True model is Poisson; 2=True model is geometric
theta_true=2; % Poisson param
phi_true=0.5; % geometric param
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% log likelihood of model M1 (Poisson): enter data and param
logf1 = @(x,theta) sum(x)*log(theta) - length(x)*theta - sum(log(factorial(x)));
% log likelihood of model M2 (Geometric): enter data and param
logf2 = @(x,phi) length(x)*log(phi) + sum(x)*log(1-phi);


for simu = 1 : 100 % 100 independent simulations
    if mod(simu,10)==0
        simu
    end
    
% data generation
switch model
    case 1 % TRUE model = Poisson
    data = poissrnd(theta_true,1,n);
    
    case 2 % TRUE model = Geometric
    data = geornd(phi_true,1,n);
end

%--- marginal likelihood of model M2 (prior is uniform on (0,1)) ---%
step2=0.01; grid2 = 0:step2:1;
y2=logf2(data,grid2);
m2 = sum(exp(y2)*step2);

%--- marginal likelihood of model M1 (improper prior) ---%
step1=0.01; grid1 = 0.01:step1:1000;
y1 = logf1(data,grid1);
m1 = sum(exp(y1)*step1); % marglike of complete data

m1end = 0;
rep = n; % number of training sets (size of trainset = n/rep)
for i = 1 : rep
    data_train = data( (1 + (i-1)*n/rep): i*n/rep );
    y1train = logf1(data_train,grid1); % only 1 observation
    m1train = sum(exp(y1train)*step1); % marglike of train sample
    
    m1end = m1end + 1/rep* (m1/m1train); % averaging, like in intrinsic BF (Berger and Pericchi, 1996)
end

% Intrinsic Bayes Factor
BF(simu) = m1end / m2;
end

min(BF)
max(BF)
switch model
    case 1
        disp('Exceptions M1')
        sum(BF<1)
    case 2
        disp('Exceptions M2')
        sum(BF>1)
end