I have a problem with the following code. In the evidence part, the value is very small so, in the end, the probabilities are not calculated. I need to normalize it, but in which part should it be?
The code in MATLAB is:
clear all; close all; clc;
randn('seed', 1234);
resistivities = [50 200 2000 1500];
thicknesses = [500 100 200];
Par_real = [resistivities, thicknesses];
dataFreq = logspace(log10(0.001), log10(1000), 100);
[Ydata, phase] = modelMT2(Par_real, dataFreq);
sigma = 0.1;
Yexp = Ydata + sigma*randn(size(Ydata));
plot(dataFreq, Yexp, '.'); hold on; plot(dataFreq, Ydata, '-')
nsamples = 20000;
R1 = 5;
R2 = 2050;
P1 = 25;
P2 = 500;
Resis = R1 + (R2-R1)*rand(nsamples, 7);
Profs = P1 + (P2-P1)*rand(nsamples, 6);
for ii=1:nsamples
par3C = [Resis(ii, 1:3), Profs(ii, 1:2)];
par4C = [Resis(ii, 1:4), Profs(ii, 1:3)];
par5C = [Resis(ii, 1:5), Profs(ii, 1:4)];
par7C = [Resis(ii, 1:7), Profs(ii, 1:6)];
Like_M3C(ii) = log_likelihood(@modelMT2, dataFreq, Yexp, sigma, par3C);
Like_M4C(ii) = log_likelihood(@modelMT2, dataFreq, Yexp, sigma, par4C);
Like_M5C(ii) = log_likelihood(@modelMT2, dataFreq, Yexp, sigma, par5C);
Like_M7C(ii) = log_likelihood(@modelMT2, dataFreq, Yexp, sigma, par7C);
end
figure()
subplot(1, 2, 1)
plot(exp(Like_M5C))
subplot(1, 2, 2)
hist(exp(Like_M5C))
Evidencia(1) = mean(exp(Like_M3C));
Evidencia(2) = mean(exp(Like_M4C));
Evidencia(3) = mean(exp(Like_M5C));
Evidencia(4) = mean(exp(Like_M7C));
Denominador = sum(Evidencia);
PPMM = Evidencia/Denominador;
fprintf('La probabilidad de los modelos : \n');
fprintf('--------------------------------\n');
fprintf('Modelo M3C: %.4f \n', PPMM(1));
fprintf('Modelo M4C: %.4f \n', PPMM(2));
fprintf('Modelo M5C: %.4f \n', PPMM(3));
fprintf('Modelo M7C: %.4f \n', PPMM(4));
figure()
model = [1, 2, 3, 4];
bar(model, PPMM), grid on
ylim([0, 1])
xlabel('Modelo')
ylabel('Probabilidad del modelo')
function [LogPDF_post] = log_likelihood(Mod, xx, data, sigma, oldpar)
erro = (Mod(oldpar, xx) - data)';
LogPDF_post = -0.5 * erro' * 1/sigma^2 * erro;
end
I have tried to normalize the likelihood as follows, but it doesn't work. It gives equal probability in all cases.
function [LogPDF_norma] = log_likelihood(Mod, xx, data, sigma, oldpar)
erro = (Mod(oldpar, xx) - data)';
LogPDF_post = -0.5 * erro' * 1/sigma^2 * erro;
LogPDF_norma = (1/max(LogPDF_post))*LogPDF_post;
end