Использование функции MLE для оценки параметров пользовательского дистрибутива - PullRequest
2 голосов
/ 10 июня 2019

Я пытаюсь использовать функцию mle() в MATLAB для оценки параметров 6-параметрического пользовательского распределения.

PDF пользовательского дистрибутива:

enter image description here

, а CDF

enter image description here

, где Γ (x, y) и Γ (x) являются верхняя неполная гамма-функция и гамма-функция соответственно. α , θ , β , a , b и c являютсяпараметры нестандартного дистрибутива. K задается как

enter image description here

Учитывая вектор данных 'data', я хочу оценить параметры α , θ , β , a, b и c.

Итак, пока я придумал этот код:

data        =  rand(20000,1); % Since I cannot upload the acutal data, we may use this
t           =  0:0.0001:0.5;    
fun         =  @(w,a,b,c) w^(a-1)*(1-w)^(b-1)*exp^(-c*w);

% to estimate the parameters
custpdf     =  @(data,myalpha,mybeta,mytheta,a,b,c)...
                ((integral(@(t)fun(t,a,b,c),0,1)^-1)*...
                mybeta*...
                igamma(myalpha,((mytheta/t)^mybeta)^(a-1))*...
                (mytheta/t)^(myalpha*mybeta+1)*...
                exp(-(mytheta/t)^mybeta-(c*(igamma(myalpha,(mytheta/t)^mybeta)/gamma(myalpha)))))...
                /...
                (mytheta*...
                gamma(myalpha)^(a+b-1)*...
                (gamma(myalpha)-igamma(myalpha,(mytheta/t)^mybeta))^(1-b));

custcdf     =  @(data,myalpha,mybeta,mytheta,a,b,c)...
                (integral(@(t)fun(t,a,b,c),0,1)^-1)*...
                integral(@(t)fun(t,a,b,c),0,igamma(myalpha,(mytheta/t)^mybeta)^mybeta/gamma(myalpha));

phat        =  mle(data,'pdf',custpdf,'cdf',custcdf,'start',0.0);

Но я получаю следующую ошибку:

Error using mlecustom (line 166)
Error evaluating the user-supplied pdf function
'@(data,myalpha,mybeta,mytheta,a,b,c)((integral(@(t)fun(t,a,b,c),0,1)^-1)*mybeta*igamma(myalpha,((mytheta/t)^mybeta)^(a-1))*(mytheta/t)^(myalpha*mybeta+1)*exp(-(mytheta/t)^mybeta-(c*(igamma(myalpha,(mytheta/t)^mybeta)/gamma(myalpha)))))/(mytheta*gamma(myalpha)^(a+b-1)*(gamma(myalpha)-igamma(myalpha,(mytheta/t)^mybeta))^(1-b))'.

Error in mle (line 245)
            phat = mlecustom(data,varargin{:});

Caused by:
    Not enough input arguments.

Я попытался просмотреть строки ошибок, но не могу понять, где на самом деле ошибка.

В какой функции отсутствуетменьше входов?Это относится к fun?Почему mle не хватает входных данных при попытке оценить параметры?

Может ли кто-нибудь помочь мне отладить ошибку?

Заранее спасибо.

1 Ответ

1 голос
/ 10 июня 2019
  • exp() - это функция, а не переменная, точный аргумент
exp^(-c*w) ---> exp(-c*w)
  • Отправная точка касается 6 parameters, а не только одного 0.1*ones(1,6)
  • В custcdf mle требуется, чтобы верхняя граница интеграла была скаляр, я сделал несколько проб и ошибок и диапазон [2~9]. для пробные некоторые значения приводят к отрицательному cdf или меньше 1 отбрасывают их.
  • Затем используйте правильный, чтобы вычислить верхнюю границу. такой же, как тот, который вы предустановили.

Я переписываю все функции, проверяю их

Код следующий

Censored = ones(5,1);% All data could be trusted 

data        =  rand(5,1); % Since I cannot upload the acutal data, we may use this

f         =  @(w,a,b,c) (w.^(a-1)).*((1-w).^(b-1)).*exp(-c.*w);
% to estimate the parameters
custpdf     =  @(t,alpha,theta,beta, a,b,c)...
                (((integral(@(w)f(w,a,b,c), 0,1)).^-1).*...
                beta.*...
                ((igamma(alpha, (theta./t).^beta)).^(a-1)).*...
                ((theta./t).^(alpha.*beta + 1 )).*...
                exp(-(((theta./t).^beta)+...
                c.*igamma(alpha, (theta./t).^beta)./gamma(alpha))))./...
                (theta.*...
                ((gamma(alpha)).^(a+b-1)).*...
                 ((gamma(alpha)-...
                 igamma(alpha, (theta./t).^beta)).^(1-b)));


custcdf = @(t,alpha,theta,beta, a,b,c)...
         ((integral(@(w)f(w,a,b,c), 0,1)).^-1).*...         
     (integral(@(w)f(w,a,b,c), 0,2));



phat = mle(data,'pdf',custpdf,'cdf',custcdf,'start', 0.1.*ones(1,6),'Censoring',Censored);

Результат

    phat = 0.1017    0.1223    0.1153    0.1493   -0.0377    0.0902
...