Исследование невязок для верификации модели

Можно изучить stats структура, которая возвращается обоими nlmefit и nlmefitsa, чтобы определить качество вашей модели. The stats структура содержит поля с условными взвешенными невязками (cwres поле) и отдельные взвешенные невязки (iwres поле). Поскольку модель принимает, что невязки обычно распределены, можно изучить невязки, чтобы увидеть, насколько хорошо это предположение удерживается.

Этот пример генерирует синтетические данные с использованием нормальных распределений. Он показывает, как выглядит статистика подгонки:

  • Хорошо при тестировании против того же типа модели, что и при генерации данных

  • Плохой при тестировании на неправильные модели данных

  1. Инициализируйте модель 2-D с 100 индивидуумами:

    nGroups = 100; % 100 Individuals
    nlmefun = @(PHI,t)(PHI(:,1)*5 + PHI(:,2)^2.*t); % Regression fcn
    REParamsSelect = [1  2]; % Both Parameters have random effect
    errorParam = .03; 
    beta0 = [ 1.5  5]; % Parameter means
    psi = [ 0.35  0; ...  % Covariance Matrix
           0   0.51 ];
    time =[0.25;0.5;0.75;1;1.25;2;3;4;5;6];
    nParameters = 2;
    rng(0,'twister') % for reproducibility
  2. Сгенерируйте данные для подбора кривой с пропорциональной моделью ошибки:

    b_i = mvnrnd(zeros(1, numel(REParamsSelect)), psi, nGroups);
    individualParameters = zeros(nGroups,nParameters);
    individualParameters(:, REParamsSelect) = ...
         bsxfun(@plus,beta0(REParamsSelect), b_i);
    
    groups = repmat(1:nGroups,numel(time),1);
    groups = vertcat(groups(:));
    
    y = zeros(numel(time)*nGroups,1);
    x = zeros(numel(time)*nGroups,1);
    for i = 1:nGroups
        idx = groups == i;
        f = nlmefun(individualParameters(i,:), time);
        % Make a proportional error model for y:
        y(idx) = f + errorParam*f.*randn(numel(f),1);
        x(idx) = time;
    end
    
    P = [ 1 0 ; 0 1 ];
  3. Подгонка данных с помощью той же функции регрессии и модели ошибки, что и генератор модели:

    [~,~,stats] = nlmefit(x,y,groups, ...
        [],nlmefun,[1 1],'REParamsSelect',REParamsSelect,...
        'ErrorModel','Proportional','CovPattern',P);
  4. Создайте графическое изображение стандартной программы путем копирования следующего определения функции и создания файла plotResiduals.m на вашем MATLAB® путь:

    function plotResiduals(stats)
    pwres = stats.pwres;
    iwres = stats.iwres;
    cwres = stats.cwres;
    figure
    subplot(2,3,1);
    normplot(pwres); title('PWRES')
    subplot(2,3,4);
    createhistplot(pwres);
    
    subplot(2,3,2);
    normplot(cwres); title('CWRES')
    subplot(2,3,5);
    createhistplot(cwres);
    
    subplot(2,3,3);
    normplot(iwres); title('IWRES')
    subplot(2,3,6);
    createhistplot(iwres); title('IWRES')
    
    function createhistplot(pwres)
    h = histogram(pwres);
    
    % x is the probability/height for each bin
    x = h.Values/sum(h.Values*h.BinWidth)
    
    % n is the center of each bin
    n = h.BinEdges + (0.5*h.BinWidth)
    n(end) = [];
    
    bar(n,x);
    ylim([0 max(x)*1.05]);
    hold on;
    x2 = -4:0.1:4;
    f2 = normpdf(x2,0,1);
    plot(x2,f2,'r');
    end
    
    end
  5. Постройте график невязок с помощью plotResiduals функция:

    plotResiduals(stats);

    Графики верхней вероятности выглядят прямо, что означает, что невязки обычно распределены. Графики нижней гистограммы совпадают с наложенным графиком нормальной плотности. Таким образом, можно сделать вывод, что модель ошибки соответствует данным.

  6. Для сравнения подгоните данные с помощью модели постоянной ошибки вместо пропорциональной модели, которая создала данные:

    [~,~,stats] = nlmefit(x,y,groups, ...
        [],nlmefun,[0 0],'REParamsSelect',REParamsSelect,...
        'ErrorModel','Constant','CovPattern',P);
    plotResiduals(stats);

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

  7. Для другого сравнения подбирайте данные к другой структурной модели, чем та, которая создала данные:

    nlmefun2 = @(PHI,t)(PHI(:,1)*5 + PHI(:,2).*t.^4);
    [~,~,stats] = nlmefit(x,y,groups, ...
        [],nlmefun2,[0 0],'REParamsSelect',REParamsSelect,...
        'ErrorModel','constant', 'CovPattern',P);
    plotResiduals(stats);

    Графики верхней вероятности не прямые. Кроме того, графики гистограммы довольно искажены по сравнению с наложенными графиками нормальной плотности. Эти невязки обычно не распределены и не совпадают с моделью.