Загрузка данных

Загрузите carsmall набор данных, который содержит измерения автомобилей, сделанных в 1970-х и в начале 1980-х.

load carsmall

Этот пример использует автомобильные измерения веса в Weight переменная.

Соответствуйте распределению 2D параметра Weibull

Во-первых, соответствуйте 2D параметру распределение Weibull к Weight.

pd = fitdist(Weight,'Weibull')

pd = 
  WeibullDistribution

  Weibull distribution
    A = 3321.64   [3157.65, 3494.15]
    B = 4.10083   [3.52497, 4.77076]

Постройте подгонку с гистограммой.

figure
histogram(Weight,8,'Normalization','pdf')
hold on
x = linspace(0,6000);
plot(x,pdf(pd,x),'LineWidth',2)
hold off

Подходящий график распределения не совпадает с гистограммой хорошо. Гистограмма не показывает выборок в области где Weight <1500. Соответствуйте распределению Weibull снова после вычитания 1500 от Weight.

pd = fitdist(Weight-1500,'Weibull')

pd = 
  WeibullDistribution

  Weibull distribution
    A = 1711.75   [1543.58, 1898.23]
    B = 1.99963   [1.70954, 2.33895]

figure
histogram(Weight-1500,8,'Normalization','pdf')
hold on
plot(x,pdf(pd,x),'LineWidth',2)
hold off

Подходящий график распределения совпадает с гистограммой лучше.

Вместо того, чтобы задать произвольное значение для предела распределения, можно задать пользовательскую функцию для распределения Weibull с тремя параметрами и оценить предел (параметр положения $$c$$).

Задайте Пользовательскую PDF для Распределения Weibull С тремя параметрами

Задайте функцию плотности вероятности для распределения Weibull с тремя параметрами.

f_def = @(x,a,b,c) (x>c).*(b/a).*(((x-c)/a).^(b-1)).*exp(-((x-c)/a).^b);

В качестве альтернативы можно использовать wblpdf функция, чтобы задать распределение Weibull с тремя параметрами.

f = @(x,a,b,c) wblpdf(x-c,a,b);

Соответствуйте распределению Weibull с тремя параметрами

Найдите MLEs для этих трех параметров при помощи mle функция. Необходимо также задать начальные значения параметров (Start аргумент значения имени) для пользовательского дистрибутива.

try
    mle(Weight,'pdf',f,'Start',[1700 2 1500])
catch ME
    disp(ME)
end

  MException with properties:

    identifier: 'stats:mle:NonpositivePdfVal'
       message: 'Custom probability function returned negative or zero values.'
         cause: {}
         stack: [12x1 struct]
    Correction: []

mle возвращает ошибку, потому что пользовательская функция возвращает неположительные значения. Эта ошибка является типичной проблемой, когда вы соответствуете нижнему пределу распределения или соответствуете распределению областью, которая имеет нулевую плотность вероятности. mle попытки некоторые значения параметров, которые имеют нулевую плотность, и затем не удаются оценить параметры. В предыдущем вызове функции, mle значения попыток $$c$$ это выше, чем минимальное значение Weight, который приводит к нулевой плотности для некоторых точек и возвращает ошибку.

Чтобы избежать этой проблемы, можно выключить опцию, которая проверяет на недопустимые значения функции, и задайте границы параметра, когда вы вызываете mle функция.

Отобразите опции по умолчанию для итеративного процесса оценки mle функция.

statset('mlecustom')

ans = struct with fields:
          Display: 'off'
      MaxFunEvals: 400
          MaxIter: 200
           TolBnd: 1.0000e-06
           TolFun: 1.0000e-06
       TolTypeFun: []
             TolX: 1.0000e-06
         TolTypeX: []
          GradObj: 'off'
         Jacobian: []
        DerivStep: 6.0555e-06
      FunValCheck: 'on'
           Robust: []
     RobustWgtFun: []
           WgtFun: []
             Tune: []
      UseParallel: []
    UseSubstreams: []
          Streams: {}
        OutputFcn: []

Замените то значение по умолчанию, с помощью структуры опций, созданной с statset функция. Задайте FunValCheck значение поля как 'off' выключить проверку достоверности для пользовательских значений функции.

opt = statset('FunValCheck','off');

Найдите MLEs этих трех параметров снова. Задайте опцию итеративного процесса (Options аргумент значения имени) и границы параметра (LowerBound и UpperBound аргументы name-value). Шкала и параметры формы должны быть положительными, и параметр положения должен быть меньшим, чем минимум выборочных данных.

params = mle(Weight,'pdf',f,'Start',[1700 2 1500],'Options',opt, ...
    'LowerBound',[0 0 -Inf],'UpperBound',[Inf Inf min(Weight)])

params = 1×3
103 ×

    1.3874    0.0015    1.7581

Постройте подгонку с гистограммой.

figure
histogram(Weight,8,'Normalization','pdf')
hold on
plot(x,f(x,params(1),params(2),params(3)),'LineWidth',2)
hold off

Подходящий график распределения совпадает с гистограммой хорошо.

Соответствуйте распределению Weibull с тремя параметрами для $$b<1$$

Если масштабный коэффициент $$b$$ меньше 1, PDF бесконечности подходов распределения Weibull около нижнего предела $$c$$ (параметр положения). Можно избежать этой проблемы путем определения подвергнутых цензуре интервалом данных, в подходящих случаях.

Загрузите cities набор данных. Данные включают оценки для девяти различных индикаторов качества жизни в 329 городах США: климат, корпус, здоровье, преступление, транспортировка, образование, искусства, воссоздание и экономика. Для каждого индикатора более высокая оценка лучше.

load cities

Найдите MLEs для седьмого индикатора (искусства).

Y = ratings(:,7);
params1 = mle(Y,'pdf',f,'Start',[median(Y) 1 0],'Options',opt)

Warning: Maximum likelihood estimation did not converge.  Iteration limit exceeded.

params1 = 1×3
103 ×

    2.7584    0.0008    0.0520

Предупреждающее сообщение указывает, что оценка не сходится. Измените опции оценки и найдите MLEs снова. Увеличьте максимальное число итераций (MaxIter) и максимальное количество оценок целевой функции (MaxFunEvals).

opt.MaxIter = 1e3;
opt.MaxFunEvals = 1e3;
params2 = mle(Y,'pdf',f,'Start',params1,'Options',opt)

Warning: Maximum likelihood estimation did not converge.  Function evaluation limit exceeded.

params2 = 1×3
103 ×

    2.7407    0.0008    0.0520

Итерация все еще не сходится, потому что PDF приближается к бесконечности около нижнего предела.

Примите что индикаторы в Y значения, округленные до ближайшего целого числа. Затем можно обработать значения в Y как подвергнутые цензуре интервалом наблюдения. Наблюдение y в Y указывает, что фактическая оценка между y–0.5 и y+0.5. Создайте матрицу, в которой каждая строка представляет интервал, окружающий каждое целое число в Y.

intervalY = [Y-0.5, Y+0.5];

Найдите MLEs снова с помощью intervalY. Чтобы соответствовать пользовательскому дистрибутиву к подвергнутому цензуре набору данных, необходимо передать и PDF и cdf к mle функция.

F = @(x,a,b,c) wblcdf(x-c,a,b);
params = mle(intervalY,'pdf',f,'cdf',F,'Start',params2,'Options',opt)

params = 1×3
103 ×

    2.7949    0.0008    0.0515

Функция находит MLEs без любых проблем сходимости. Эта подгонка основана на подходящих вероятностях к интервалам, таким образом, она не сталкивается с проблемой плотности приближающаяся бесконечность в одной точке. Можно использовать этот подход только, когда преобразование данных к подвергнутой цензуре интервалом версии является соответствующим.

Постройте график результатов.

figure
histogram(Y,'Normalization','pdf')
hold on
x = linspace(0,max(Y));
plot(x,f(x,params(1),params(2),params(3)),'LineWidth',2)
hold off

Подходящий график распределения совпадает с гистограммой хорошо.