Maximum Likelihood

Hello,

Can anyone please give me some tips on a code I have written so far? I am trying to write a MCMC simulation (my first one) that will calculate the maximum likelihood from a chi squared value at any P = (a,b). I then randomly jump to a new point P' = (a',b') and calculate the maximum likelihood there and compare these two values labelled R. R is then compared to a uniform number between 0 and 1 to see which way I progress. I want this process to be repeated till I converge on the maximum likelihood.

clear all % clears all previous data

s = load('domain\data.txt');

% loads in data file

% s(:,1) data number

% s(:,2) x variables

% s(:,3) y variables

x = s(:,2);

y = s(:,3);

for z = 1:10000

chi2PN = zeros(100,100);

% STEP 1 Sampling an initial random point P1 = (a1,b1)

a = 9.50;

b = 10.49;

c = a + (b-a) * rand; % random a value

d = 4.50;

e = 5.49;

f = d + (d-e) * rand; % random b value

% rand command for uniform random number generation

% STEP 2 Evaluating the value of chi squared at point P1 = (a1,b1)

chi2P1 = sum((y - (c+(f.*x))).^2); % a is equivalent to c and b is equivalent to f

% STEP 3 Sampling tentative new point P' = (a',b') % Sampling random "jump" to add to P1. This time randn % for sampling from a normal distribution

g = 0;

h = 1;

i = abs(g + (h-g)*randn);

j = c+i;

k = f+i;

% STEP 4 Evaluating the value of chi squared at point P' = (a',b')

chi2PN = sum((y - (j+(k.*x))).^2); % a is equivalent to j and b is equivalent to k

% STEP 5 Calculating R a ratio of the ML which is as follows

R = exp((-1/2)*(chi2PN-chi2P1));

l = 0;

m = 1;

n = l + (m-l) * rand; % random number between 0 and 1 drawn from a uniform distribution

if (n<R);

    chi2PN = chi2P1 + 1;
else (n>=R);
    chi2PN = chi2P1

end

end

I know this is a complicated question but if anyone has any example code it would be greatly appreciated...Thanks

[the cyclist]

Can you please give a little more detail about what info you are seeking? Is your code producing error messages, or the wrong answer? Or is it correct but slow? I don't understand why you need "sample code" when you seem to have the code all written.

[Dominic Lawson]

Hello,

Thanks for getting back to me. I don't think I have written the if else loop correctly and I need to work on making my whole programme loop. Any pointers please?

[Sean de Wolski]

Don't call:

>>clear all

It does bad things behind the scenes.

Use:

>>clearvars;

if you want to get rid of variables. Or better yet, write a your program as a function!

[Dominic Lawson]

Thanks,

How does one write a program as a function?

[Sean de Wolski]

function [outputs] = function_name(inputs)

%save the file as 'function_name'

doc function

%for more info