Maximum Likelihood

17 views (last 30 days)
Dominic Lawson
Dominic Lawson on 1 Apr 2011
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
  5 Comments
Dominic Lawson
Dominic Lawson on 1 Apr 2011
Thanks,
How does one write a program as a function?
Sean de Wolski
Sean de Wolski on 1 Apr 2011
function [outputs] = function_name(inputs)
%save the file as 'function_name'
doc function
%for more info

Sign in to comment.

Answers (0)

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!