MCMC

Hello,

I am currently writing a MCMC code. Here is what I have so far:-

clear all

s = load('domain name on my computer.txt');

x = s(:,2);

y = s(:,3);

for q = 1:1000

    a = 9.9840;
    b = 9.9939;
    c(q) = a + (b-a) * rand;
    d = 5.0040;
    e = 5.0139;
    f(q) = d + (d-e) * rand;
    chi2P1(q) = sum((y - (c(q)+(f(q).*x))).^2);
    g = 0;
    h = 1;
    i(q) = abs(g + (h-g)*randn);
    j(q) = c(q)+i(q);
    k(q) = f(q)+i(q);
    chi2PN(q) = sum((y - (j(q)+(k(q).*x))).^2);
    if (chi2PN(q)<=chi2P1(q));
    elseif (chi2PN(q)>=chi2P1(q));
    end
end

At the point of the if loop I am stuck. If j and k values give a smaller chi2PN value than chi2P1 (dependent on c and f) how would I make j and k become the new c and f?

Also if this isn't the case how would I disregard the j and k values?

Hope this makes sense and thanks for taking a look at my problem.

Answer 1

[Oleg Komarov]

Hi, you don't any loop there:

x = s(:,2);
y = s(:,3);

q = 1000;

a = 9.9840;
b = 9.9939;
c = a + (b-a) * rand(1,q);
d = 5.0040;
e = 5.0139;
f = d + (d-e) * rand(1,q);
chi2P1 = sum((y - c + f .*x ).^2);

g = 0;
h = 1;
i = abs(g + (h-g)*randn(1,q));
j = c + i;
k = f + i;
chi2PN = sum((y - j + k.*x).^2);

Now you can locate values of chi2PN smaller equal than chi2P1. Do you want to replace the values of the latter array with those in the former or the other way around?

idx = chi2PN < chi2P1;
chi2P1(idx) = chi2PN(idx);

[Dominic Lawson]

thanks very much for the help.

I would like to replace values of chi2P1 with values of chi2PN that are smaller than the values of chi2P1.

Does that make sense?

[Oleg Komarov]

I edited the post. Look at the last line.

[Dominic Lawson]

thank you