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
5 个评论
the cyclist
2011-4-1
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
2011-4-1
Sean de Wolski
2011-4-1
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
2011-4-1
Sean de Wolski
2011-4-1
function [outputs] = function_name(inputs)
%save the file as 'function_name'
doc function
%for more info
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Random Number Generation 的更多信息
产品
2011-4-1
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
