7 added 928 characters in body
source | link
    %%%%%% 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
source | link

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
source | link

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
source | link
3 added 84 characters in body
source | link
2 added 16 characters in body
source | link
1
source | link