function [para,standard_deviation,fv]=my_mle(fun,para0,varargin)

%estimate parameters and standard errors when using maximium likelihood estimation(MLE)

%input

%fun: a function defined by users for calculating log probability density function (pdf) and negative sum of logarithm of pdf

%para0 : given initial parameters

% varargin: other needed parameters required by fun

%output

%para: estimated parameters

% standard_deviation: standard deviations of estimated parameters

% fv: maximized likelihood function value

%%%%%%%%%%%

%example1: estimate mean and standard deviation by realizations of a

%random variable which is normally distributed

%function f=mynormpdfsum(x,num,y)

%yy=1/sqrt(2*pi)/x(2)*exp(-(y-x(1)).^2/2/x(2)^2);

%if num==1 %(note: it must be set to 1)

%f=log(yy);

%else f=-sum(log(yy));end

%%%%%%%%%

%y=2+3*randn(5000,1);

%[para,standard_deviation]=my_mle('mynormpdfsum',[0;2],y)

%%%%%%%%%

% example2:estimate coefficients in a linear regression

% clear

% x=randn(500,1);

% y=2+3*x+randn(500,1);

% [para,standard_deviation,fv]=my_mle('mynormpdfsum001',[1;2;3],y,x)

% %%%%%%%%%%%

% function f=mynormpdfsum001(para,num,y,x)

% yy=1/sqrt(2*pi)/para(3)*exp(-(y-para(1)-para(2)*x).^2/2/para(3)^2);

% if num==1 %(note: it must be set to 1)

% f=log(yy);

% else f=-sum(log(yy));end

para0=para0(:);

[para,fv]=fminsearch(fun,para0,[],2,varargin{:});

fv=-fv;

d=numericalfirstderivative(fun,para,1,varargin{:});

standard_deviation=sqrt(diag(pinv(d'*d)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f=numericalfirstderivative(fun,parameter,varargin)

% input:

% fun: the name of a function

% parameter: given parameter with respect to which first-order derivative is calculated

% varargin: other needed parameters required by fun

% output:

% f: numerical first order derivative of fun at parameter

n=length(parameter);

for i=1:n

    a=zeros(n,1);

    a(i)=min(parameter(i)*1e-6,1e-5);

    y1(:,i)=feval_r(fun,parameter+a,varargin{:});

    y2(:,i)=feval_r(fun,parameter-a,varargin{:});

    f(:,i)=(y1(:,i)-y2(:,i))/2/a(i);

end