7 added 928 characters in body
    %%%%%% MCMC - Metropolis Algorithm %%%%%%%
% Declare parameters of array
nb = 1e5;
% Metropolis loop
% Index of main loop
i = 1;
% Reject number
j = 0;
% Accept number
k = 0;
% First solution : generate% inStarting oncevalues
%pStart = amp*rand(6,1);
% Second solution
l1 = [15[10, 10];3];
l2 = [5[10, 4];3];
l3 = [10, 5];3];
l4 = [10, 5];3];
l5 = [2[10, 2];3];
l6 = [2[10, 2];3];

% Initial param array
paramsArray = zeros(nb,6);
paramsArray(1,1:6) = [l1(1) + l1(2)*rand(1), l2(1) + l2(2)*rand(1), ...
l3(1) + l3(2)*rand(1), l4(1) + l4(2)*rand(1), ...
l5(1) + l5(2)*rand(1), l6(1) + l6(2)*rand(1)];
% Init parameters vector
pStart = paramsArray(1,1:6);';
% varianceStandard deviation vectors for proposal distributon
gamFactor = 5;
gamMatrix = gamFactor*diagdiag([randn[l1(62),1 l2(2), l3(2), l4(2), l5(2), l6(2)]);
% Starting prior distribution f(p|d)
w_x = Crit_J(pStart,D)*exp(-((pStart(1)-l1(1))^2/(2*l1(2)^2)+(pStart(2)-l2(1))^2/(2*l2(2)^2)+ ...
(pStart(3)-l3(1))^2/(2*l3(2)^2)+(pStart(4)-l4(1))^2/(2*l4(2)^2)+ ...
(pStart(5)-l5(1))^2/(2*l5(2)^2)+(pStart(6)-l6(1))^2/(2*l6(2)^2)));
while (i <= nb)
if (i == 1)
% First random
%ptest(1:6) = pStart(1:6) + gamVec*rand(6,1);
ptest = pStart;
ptest = abs(pStart + gamMatrix*randn(6,1));else
else
% Other random
ptemp = paramsArray(i,1:6);
ptest = abs(ptemppStart + gamMatrix*randn(6,1));
end
% Upper term for acceptance ratio
w_y = Crit_J(ptest,D);*exp(-((ptest(1)-l1(1))^2/(2*l1(2)^2)+(ptest(2)-l2(1))^2/(2*l2(2)^2)+ ...
% Ratio acceptance                          (ptest(3)-l3(1))^2/(2*l3(2)^2)+(ptest(4)-l4(1))^2/(2*l4(2)^2)+ ...
log_prob = log                          (w_yptest(5)-l5(1))^2/w_x(2*l5(2)^2)+(ptest(6)-l6(1))^2/(2*l6(2)^2)));
% Random number forGenerate acceptanceu testuniformly
u = rand(1);
% Ratio acceptance
log_prob = log(w_y/w_x);
% Test acceptation
if (log(u) < log_prob)
% Assing new paramsArray
paramsArray(i,1:6) = ptest(1:6);
w_x = Crit_J(paramsArray(i,1:6),D);*exp(-((ptest(1)-l1(1))^2/(2*l1(2)^2)+(ptest(2)-l2(1))^2/(2*l2(2)^2)+ ...
i = i+1;
k = k+1;
else
if (i ~= 1)
% Random values
randNew = rand       (ptest(3)-l3(1);)^2/(2*l3(2)^2)+(ptest(4)-l4(1))^2/(2*l4(2)^2)+ ...
% Assing new paramsArray
if (randNew < exp(log_prob))
paramsArray(i,1:6) = ptest(1:6);
i = i+1;
(ptest(5)-l5(1))^2/(2*l5(2)^2)+(ptest(6)-l6(1))^2/(2*l6(2)^2)));
ki = k+1;i+1;
k = elsek+1;
else
% Assing to previous
if paramsArray(i, ~= 1:6)
=     paramsArray(i-1,1:6);
w_x = Crit_J(paramsArray(i-1,1:6),D);
i = i+1;
j = j+1;
end
end
end
end
disp('acceptationt : ratio');
disp(k/nb)
disp('reject : ratio');
disp(j/nb)
% Display mean of different parameters
disp('Parameters with Metropolis-Hastings :');
mean(paramsArray(:,1))
mean(paramsArray(:,2))
mean(paramsArray(:,3))
mean(paramsArray(:,4))
mean(paramsArray(:,5))
mean(paramsArray(:,6))

% Metropolis loop
% Index of main loop
i = 1;
% Reject number
j = 0;
% Accept number
k = 0;
% First solution : generate in once
%pStart = amp*rand(6,1);
% Second solution
l1 = [15, 10];
l2 = [5, 4];
l3 = [10, 5];
l4 = [10, 5];
l5 = [2, 2];
l6 = [2, 2];

% Initial param array
paramsArray = zeros(nb,6);
paramsArray(1,1:6) = [l1(1) + l1(2)*rand(1), l2(1) + l2(2)*rand(1), ...
l3(1) + l3(2)*rand(1), l4(1) + l4(2)*rand(1), ...
l5(1) + l5(2)*rand(1), l6(1) + l6(2)*rand(1)];
% Init parameters vector
pStart = paramsArray(1,1:6);
% variance vectors for proposal distributon
gamFactor = 5;
gamMatrix = gamFactor*diag([randn(6,1)]);
% Starting prior distribution f(p|d)
w_x = Crit_J(pStart,D);
while (i <= nb)
if (i == 1)
% First random
%ptest(1:6) = pStart(1:6) + gamVec*rand(6,1);
ptest = abs(pStart + gamMatrix*randn(6,1));
else
% Other random
ptemp = paramsArray(i,1:6);
ptest = abs(ptemp + gamMatrix*randn(6,1));
end
% Upper term for acceptance ratio
w_y = Crit_J(ptest,D);
% Ratio acceptance
log_prob = log(w_y/w_x);
% Random number for acceptance test
u = rand(1);
% Test acceptation
if (log(u) < log_prob)
% Assing new paramsArray
paramsArray(i,1:6) = ptest(1:6);
w_x = Crit_J(paramsArray(i,1:6),D);
i = i+1;
k = k+1;
else
if (i ~= 1)
% Random values
randNew = rand(1);
% Assing new paramsArray
if (randNew < exp(log_prob))
paramsArray(i,1:6) = ptest(1:6);
i = i+1;
k = k+1;
else
% Assing to previous
paramsArray(i,1:6) = paramsArray(i-1,1:6);
w_x = Crit_J(paramsArray(i,1:6),D);
i = i+1;
j = j+1;
end
end
end
end
disp('acceptationt : ratio');
disp(k/nb)
disp('reject : ratio');
disp(j/nb)

    %%%%%% MCMC - Metropolis Algorithm %%%%%%%
% Declare parameters of array
nb = 1e5;
% Metropolis loop
% Index of main loop
i = 1;
% Reject number
j = 0;
% Accept number
k = 0;
% Starting values
l1 = [10, 3];
l2 = [10, 3];
l3 = [10, 3];
l4 = [10, 3];
l5 = [10, 3];
l6 = [10, 3];

% Initial param array
paramsArray = zeros(nb,6);
paramsArray(1,1:6) = [l1(1) + l1(2)*rand(1), l2(1) + l2(2)*rand(1), ...
l3(1) + l3(2)*rand(1), l4(1) + l4(2)*rand(1), ...
l5(1) + l5(2)*rand(1), l6(1) + l6(2)*rand(1)];
% Init parameters vector
pStart = paramsArray(1,1:6)';
% Standard deviation vectors for proposal distributon
gamMatrix = diag([l1(2), l2(2), l3(2), l4(2), l5(2), l6(2)]);
% Starting prior distribution f(p|d)
w_x = Crit_J(pStart,D)*exp(-((pStart(1)-l1(1))^2/(2*l1(2)^2)+(pStart(2)-l2(1))^2/(2*l2(2)^2)+ ...
(pStart(3)-l3(1))^2/(2*l3(2)^2)+(pStart(4)-l4(1))^2/(2*l4(2)^2)+ ...
(pStart(5)-l5(1))^2/(2*l5(2)^2)+(pStart(6)-l6(1))^2/(2*l6(2)^2)));
while (i <= nb)
if (i == 1)
% First random
ptest = pStart;
else
% Other random
ptest = abs(pStart + gamMatrix*randn(6,1));
end
% Upper term for acceptance ratio
w_y = Crit_J(ptest,D)*exp(-((ptest(1)-l1(1))^2/(2*l1(2)^2)+(ptest(2)-l2(1))^2/(2*l2(2)^2)+ ...
(ptest(3)-l3(1))^2/(2*l3(2)^2)+(ptest(4)-l4(1))^2/(2*l4(2)^2)+ ...
(ptest(5)-l5(1))^2/(2*l5(2)^2)+(ptest(6)-l6(1))^2/(2*l6(2)^2)));
% Generate u uniformly
u = rand(1);
% Ratio acceptance
log_prob = log(w_y/w_x);
% Test acceptation
if (log(u) < log_prob)
% Assing new paramsArray
paramsArray(i,1:6) = ptest(1:6);
w_x = Crit_J(paramsArray(i,1:6),D)*exp(-((ptest(1)-l1(1))^2/(2*l1(2)^2)+(ptest(2)-l2(1))^2/(2*l2(2)^2)+ ...
(ptest(3)-l3(1))^2/(2*l3(2)^2)+(ptest(4)-l4(1))^2/(2*l4(2)^2)+ ...
(ptest(5)-l5(1))^2/(2*l5(2)^2)+(ptest(6)-l6(1))^2/(2*l6(2)^2)));
i = i+1;
k = k+1;
else
% Assing to previous
if (i ~= 1)
paramsArray(i,1:6) = paramsArray(i-1,1:6);
i = i+1;
j = j+1;
end
end
end
disp('acceptationt : ratio');
disp(k/nb)
disp('reject : ratio');
disp(j/nb)
% Display mean of different parameters
disp('Parameters with Metropolis-Hastings :');
mean(paramsArray(:,1))
mean(paramsArray(:,2))
mean(paramsArray(:,3))
mean(paramsArray(:,4))
mean(paramsArray(:,5))
mean(paramsArray(:,6))

6 edited title

and my cost function (assimilated to Likelihood function) :

% Function  of cost
function cost = Crit_J(p,D)

% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = (1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6));
model = Model(:);
d = D(:);
% Introduce H matrix
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
cost = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));

end


If you could see the error ...

If you could see the error ...

and my cost function (assimilated to Likelihood function) :

% Function  of cost
function cost = Crit_J(p,D)

% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = (1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6));
model = Model(:);
d = D(:);
% Introduce H matrix
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
cost = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));

end


If you could see the error ...

5 edited title

# Metropolis Hastings - Acceptance ratio, proposal and lkelihood

Thanks and sorry if answers are obvious.UPDATE 1: Here is below the Metropolis-Hastings I am using :

% Metropolis loop
% Index of main loop
i = 1;
% Reject number
j = 0;
% Accept number
k = 0;
% First solution : generate in once
%pStart = amp*rand(6,1);
% Second solution
l1 = [15, 10];
l2 = [5, 4];
l3 = [10, 5];
l4 = [10, 5];
l5 = [2, 2];
l6 = [2, 2];

% Initial param array
paramsArray = zeros(nb,6);
paramsArray(1,1:6) = [l1(1) + l1(2)*rand(1), l2(1) + l2(2)*rand(1), ...
l3(1) + l3(2)*rand(1), l4(1) + l4(2)*rand(1), ...
l5(1) + l5(2)*rand(1), l6(1) + l6(2)*rand(1)];
% Init parameters vector
pStart = paramsArray(1,1:6);
% variance vectors for proposal distributon
gamFactor = 5;
gamMatrix = gamFactor*diag([randn(6,1)]);
% Starting prior distribution f(p|d)
w_x = Crit_J(pStart,D);
while (i <= nb)
if (i == 1)
% First random
%ptest(1:6) = pStart(1:6) + gamVec*rand(6,1);
ptest = abs(pStart + gamMatrix*randn(6,1));
else
% Other random
ptemp = paramsArray(i,1:6);
ptest = abs(ptemp + gamMatrix*randn(6,1));
end
% Upper term for acceptance ratio
w_y = Crit_J(ptest,D);
% Ratio acceptance
log_prob = log(w_y/w_x);
% Random number for acceptance test
u = rand(1);
% Test acceptation
if (log(u) < log_prob)
% Assing new paramsArray
paramsArray(i,1:6) = ptest(1:6);
w_x = Crit_J(paramsArray(i,1:6),D);
i = i+1;
k = k+1;
else
if (i ~= 1)
% Random values
randNew = rand(1);
% Assing new paramsArray
if (randNew < exp(log_prob))
paramsArray(i,1:6) = ptest(1:6);
i = i+1;
k = k+1;
else
% Assing to previous
paramsArray(i,1:6) = paramsArray(i-1,1:6);
w_x = Crit_J(paramsArray(i,1:6),D);
i = i+1;
j = j+1;
end
end
end
end
disp('acceptationt : ratio');
disp(k/nb)
disp('reject : ratio');
disp(j/nb)


If you could see the error ...

# Acceptance ratio, proposal and lkelihood

Thanks and sorry if answers are obvious.

# Metropolis Hastings - Acceptance ratio, proposal and lkelihood

UPDATE 1: Here is below the Metropolis-Hastings I am using :

% Metropolis loop
% Index of main loop
i = 1;
% Reject number
j = 0;
% Accept number
k = 0;
% First solution : generate in once
%pStart = amp*rand(6,1);
% Second solution
l1 = [15, 10];
l2 = [5, 4];
l3 = [10, 5];
l4 = [10, 5];
l5 = [2, 2];
l6 = [2, 2];

% Initial param array
paramsArray = zeros(nb,6);
paramsArray(1,1:6) = [l1(1) + l1(2)*rand(1), l2(1) + l2(2)*rand(1), ...
l3(1) + l3(2)*rand(1), l4(1) + l4(2)*rand(1), ...
l5(1) + l5(2)*rand(1), l6(1) + l6(2)*rand(1)];
% Init parameters vector
pStart = paramsArray(1,1:6);
% variance vectors for proposal distributon
gamFactor = 5;
gamMatrix = gamFactor*diag([randn(6,1)]);
% Starting prior distribution f(p|d)
w_x = Crit_J(pStart,D);
while (i <= nb)
if (i == 1)
% First random
%ptest(1:6) = pStart(1:6) + gamVec*rand(6,1);
ptest = abs(pStart + gamMatrix*randn(6,1));
else
% Other random
ptemp = paramsArray(i,1:6);
ptest = abs(ptemp + gamMatrix*randn(6,1));
end
% Upper term for acceptance ratio
w_y = Crit_J(ptest,D);
% Ratio acceptance
log_prob = log(w_y/w_x);
% Random number for acceptance test
u = rand(1);
% Test acceptation
if (log(u) < log_prob)
% Assing new paramsArray
paramsArray(i,1:6) = ptest(1:6);
w_x = Crit_J(paramsArray(i,1:6),D);
i = i+1;
k = k+1;
else
if (i ~= 1)
% Random values
randNew = rand(1);
% Assing new paramsArray
if (randNew < exp(log_prob))
paramsArray(i,1:6) = ptest(1:6);
i = i+1;
k = k+1;
else
% Assing to previous
paramsArray(i,1:6) = paramsArray(i-1,1:6);
w_x = Crit_J(paramsArray(i,1:6),D);
i = i+1;
j = j+1;
end
end
end
end
disp('acceptationt : ratio');
disp(k/nb)
disp('reject : ratio');
disp(j/nb)


If you could see the error ...

4 added 1 character in body
3 added 84 characters in body
2 added 16 characters in body
1