fminsearch подвергнут цензуре MATLAB

РЕДАКТИРОВАТЬ Я обнаружил, что y(c) в моем скрипте дает эту ошибку:

Subscript indices must either be real positive integers or logicals.

но в примере сценария y(c) печатает цензурированные значения y.

Я пытаюсь использовать скрипт, который я нашел для анализа цензурированных данных. Сам скрипт работает нормально, но когда я пытаюсь использовать его со своими данными, я получаю сообщение об ошибке.

Это пример кода:

% I have some (x,y) data
n = 100;
x = 10*rand(n,1);
y = 5 + .5*x + randn(n,1);
plot(x,y,'o','Color',[.8 .8 .8]);

% But it's censored, I can't observe values larger than 8
c = y>8
o = min(y,8);
line(x,o,'Marker','o','Color','b','LineStyle','none')

% If I fit a line to the data I observe, no good
b = polyfit(x,o,1)
s = norm(o-polyval(b,x))/sqrt(n)
xx = linspace(0,10);
line(xx,polyval(b,xx),'Color','r')

% Instead I need a likelihood function that taket it censoring into account
nloglik = @(p) - sum(log(normpdf(o(~c),p(1)*x(~c)+p(2),p(3)))) ...
               - sum(log(1-normcdf(o(c),p(1)*x(c)+p(2),p(3))));

nloglik = @(p) - sum(log(normpdf(tof(~c),p(1)*z(~c)+p(2),p(3)))) ...
           - sum(log(1-normcdf(tof(c),p(1)*z(c)+p(2),p(3))));
p = fminsearch(nloglik,[b,s])

Вот мой код:

load('UV.mat') % is an 18 column array

for i = 1 : length(UV{1,7})
    if UV{1,7}(i) ~= 0
        x(i)=log10(UV{1,4}(i));
        y(i) = UV{1,7}(i);
    end
end
c=zeros(length(y),1); % c stands for censored
for i=1:length(y)
    if UV{1,8}{i} == 'u' % u stands for upper limit
        c(i)=1;
    end
end

b = polyfit(x,y,1)
s = norm(y-polyval(b,x))/sqrt(length(x))
nloglik = @(p) - sum(log(normpdf(y(~c),p(1)*x(~c)+p(2),p(3)))) ...
               - sum(log(1-normcdf(y(c),p(1)*x(c)+p(2),p(3))));
p = fminsearch(nloglik,[b,s])

Ошибка:

Subscript indices must either be real positive integers or logicals.

Error in @(p)-sum(log(normpdf(y(~c),p(1)*x(~c)+p(2),p(3))))-sum(log(1-normcdf(y(c),p(1)*x(c)+p(2),p(3))))


Error in fminsearch (line 191)
fv(:,1) = funfcn(x,varargin{:});

Error in zWithCensoring (line 22)  %zWithCensoring is the name of my script
p = fminsearch(nloglik,[b,s])

Я пытался отладить его, но p, похоже, вызывается до того, как он определен. Я обнаружил, что nloglik — это функция с входным параметром p.

Как мой скрипт выдает эту ошибку, а не пример скрипта? И как мне обойти эту ошибку?


person Matt Majic    schedule 15.04.2015    source источник
comment
Я не могу запустить Matlab прямо сейчас, но я думаю, что это происходит, если x и y - это массивы строк, но должны быть массивами столбцов. Выполните транспонирование на них с помощью .' перед выполнением подгонки.   -  person Buck Thorn    schedule 15.04.2015
comment
Спасибо, я пытался их переставить, но все равно выдает ту же ошибку.   -  person Matt Majic    schedule 15.04.2015


Ответы (2)


Измените переменную c на логическую, т.е.:

c = logical(c)

прежде чем использовать его для индексации.

В качестве альтернативы создайте 'c' как логичное для начала:

c= false(length(y),1); % c stands for censored
for i=1:length(y)
    if UV{1,8}{i} == 'u' % u stands for upper limit
        c(i)=true;
    end
end

Если это не сработает, публикация рабочего кода поможет решить проблему.

person bendervader    schedule 15.04.2015

y(c) выдавал ошибку, не уверен, почему это сработало в примере скрипта. Вручную сделал ya и yc, xa и xc равными y(~c), y(c), x(~c) и x(c) соответственно, и теперь это работает

РЕДАКТИРОВАТЬ:

Я сделал ya, yc, xa, xc вот так:

for i =1:length(y)
    if UV{1,8}{i} == 'u'
        yc(j)=y(i);
        xc(j)=x(i);
        j=j+1;
    else
        ya(k)=y(i);
        xa(k)=x(i);
        k=k+1;
    end
end

и подставил их для y(~c), y(c), x(~c) и x(c) в функцию nloglik

person Matt Majic    schedule 15.04.2015