使用 fminsearch 进行参数估计
Using fminsearch for parameter estimation
我试图找到高斯分布的对数最大似然估计,以便估计参数。
我知道 Matlab 有一个内置函数,它通过拟合高斯分布来执行此操作,但我需要使用 logMLE 来执行此操作,以便稍后将此方法扩展到其他分布。
所以这里是 gaussian dist 的对数似然函数:
Gaussian Log MLE
我使用这段代码通过 fminsearch 估计一组变量 (r) 的参数。但我的搜索没有涵盖,我也不完全明白问题出在哪里:
clear
clc
close all
%make random numbers with gaussian dist
r=[2.39587291079469
1.57478022109723
-0.442284350603745
4.39661178526569
7.94034385633171
7.52208574723178
5.80673144943155
-3.11338531920164
6.64267230284774
-2.02996003947964];
% mu=2 sigma=3
%introduce f
f=@(x,r)-(sum((-0.5.*log(2*3.14.*(x(2))))-(((r-(x(2))).^2)./(2.*(x(1))))))
fun = @(x)f(x,r);
% starting point
x0 = [0,0];
[y,fval,exitflag,output] = fminsearch(fun,x0)
f =
@(x,r)-(sum((-0.5.*log(2*3.14.*(x(2))))-(((r-(x(2))).^2)./(2.*(x(1))))))
Exiting: Maximum number of function evaluations has been exceeded
- increase MaxFunEvals option.
Current function value: 477814.233176
y = 1×2
1.0e+-3 *
0.2501 -0.0000
fval = 4.7781e+05 + 1.5708e+01i
exitflag = 0
output =
iterations: 183
funcCount: 400
algorithm: 'Nelder-Mead simplex direct search'
message: 'Exiting: Maximum number of function evaluations has been exceeded↵ - increase MaxFunEvals option.↵ Current function value: 477814.233176 ↵'
Rewrite f as follows:
function y = g(x, r)
n = length(r);
log_part = 0.5.*n.*log(x(2).^2);
sum_part = ((sum(r-x(1))).^2)./(2.*x(2).^2);
y = log_part + sum_part;
end
Use fmincon
instead of fminsearch
because standard deviation is
always a positif number.
Set standard deviation lower bound to zero 0
整个代码如下:
%make random numbers with gaussian dist
r=[2.39587291079469
1.57478022109723
-0.442284350603745
4.39661178526569
7.94034385633171
7.52208574723178
5.80673144943155
-3.11338531920164
6.64267230284774
-2.02996003947964];
% mu=2 sigma=3
fun = @(x)g(x, r);
% starting point
x0 = [0,0];
% borns
lb = [-inf, 0];
ub = [inf, inf];
[y, fval] = fmincon(fun,x0,[],[],[],[],lb,ub, []);
function y = g(x, r)
n = length(r);
log_part = 0.5.*n.*log(x(2).^2);
sum_part = ((sum(r-x(1))).^2)./(2.*x(2).^2);
y = log_part + sum_part;
end
解决方案
y = [3.0693 0.0000]
为了更好的估计直接使用mle()
代码非常简单:
y = mle(r,'distribution','normal')
解决方案
y = [3.0693 3.8056]
我试图找到高斯分布的对数最大似然估计,以便估计参数。 我知道 Matlab 有一个内置函数,它通过拟合高斯分布来执行此操作,但我需要使用 logMLE 来执行此操作,以便稍后将此方法扩展到其他分布。 所以这里是 gaussian dist 的对数似然函数: Gaussian Log MLE
我使用这段代码通过 fminsearch 估计一组变量 (r) 的参数。但我的搜索没有涵盖,我也不完全明白问题出在哪里:
clear
clc
close all
%make random numbers with gaussian dist
r=[2.39587291079469
1.57478022109723
-0.442284350603745
4.39661178526569
7.94034385633171
7.52208574723178
5.80673144943155
-3.11338531920164
6.64267230284774
-2.02996003947964];
% mu=2 sigma=3
%introduce f
f=@(x,r)-(sum((-0.5.*log(2*3.14.*(x(2))))-(((r-(x(2))).^2)./(2.*(x(1))))))
fun = @(x)f(x,r);
% starting point
x0 = [0,0];
[y,fval,exitflag,output] = fminsearch(fun,x0)
f =
@(x,r)-(sum((-0.5.*log(2*3.14.*(x(2))))-(((r-(x(2))).^2)./(2.*(x(1))))))
Exiting: Maximum number of function evaluations has been exceeded
- increase MaxFunEvals option.
Current function value: 477814.233176
y = 1×2
1.0e+-3 *
0.2501 -0.0000
fval = 4.7781e+05 + 1.5708e+01i
exitflag = 0
output =
iterations: 183
funcCount: 400
algorithm: 'Nelder-Mead simplex direct search'
message: 'Exiting: Maximum number of function evaluations has been exceeded↵ - increase MaxFunEvals option.↵ Current function value: 477814.233176 ↵'
Rewrite f as follows:
function y = g(x, r)
n = length(r);
log_part = 0.5.*n.*log(x(2).^2);
sum_part = ((sum(r-x(1))).^2)./(2.*x(2).^2);
y = log_part + sum_part;
end
Use
fmincon
instead offminsearch
because standard deviation is always a positif number.Set standard deviation lower bound to zero
0
整个代码如下:
%make random numbers with gaussian dist
r=[2.39587291079469
1.57478022109723
-0.442284350603745
4.39661178526569
7.94034385633171
7.52208574723178
5.80673144943155
-3.11338531920164
6.64267230284774
-2.02996003947964];
% mu=2 sigma=3
fun = @(x)g(x, r);
% starting point
x0 = [0,0];
% borns
lb = [-inf, 0];
ub = [inf, inf];
[y, fval] = fmincon(fun,x0,[],[],[],[],lb,ub, []);
function y = g(x, r)
n = length(r);
log_part = 0.5.*n.*log(x(2).^2);
sum_part = ((sum(r-x(1))).^2)./(2.*x(2).^2);
y = log_part + sum_part;
end
解决方案
y = [3.0693 0.0000]
为了更好的估计直接使用mle()
代码非常简单:
y = mle(r,'distribution','normal')
解决方案
y = [3.0693 3.8056]