0

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
Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131

0 Answers0