テストページ

テストページです。
this is test.

%Simple logistic regression model
%y=b0 + b1*x
function main()
%format long
format short
global n
global datax
global datay

hold on
axis([-0.5,2,-0.25,1.3])

x = -1:0.01:2;  %vector for display
datax=[1 0.4150;
       1 0.5797;
       1 0.7076;
       1 0.8865;
       1 1.0086 ];
datay=[0.120  0.333  0.522  0.875  0.880 ]';

n=length(datay);
fprintf('\n The number of data : %d\n\n',n)
disp('      y    :    x_0       x_1')
disp([datay datax])


bfmin = fminsearch(@loglike,[-20 20]);
yfmin = f(x,bfmin);

binit = [10 20];    %initial value of Fisher score method.
bfish = fisherscore(binit);
yfish = f(x,bfish);

loglikefish = loglike(bfish);
loglikefmin = loglike(bfmin);

xlabel('impulse level')
ylabel('response ratio')
plot(datax(:,2),datay(:),'ko')
%The red is the graph with b estimated with fminsearch.
plot(x,yfmin,'r')
%The yellow is the graph with Fisher score method.
plot(x,yfish)
hold off

fprintf('With Fisher score method : [b0 b1] = [%5.10f %5.10f]   \n',bfish(1),bfish(2))
fprintf('With fminsearch :          [b0 b1] = [%5.10f %5.10f]\n\n',bfmin(1),bfmin(2))

fprintf('Maximum loglikelihood with Fisher     : %5.10f\n',loglikefish)
fprintf('Maximum loglikelihood with fminsearch : %5.10f\n',loglikefmin)
%Here are subfunctions.


%The log likelihood function.
function l = loglike(b)
    global datax
    global datay
    global n
    l = 0;
    for i=1:n
        l = l -(datay(i).*b*datax(i,:)' - log(1+exp(b*datax(i,:)')));
    end

%Fisher score method
function b = fisherscore(b)
    global datax
    global datay
    global n

    for j=1:3
         pi = zeros(1,n);
        for i=1:n
            pi(i) = exp(datax(i,:)*b') ./ (1+exp(datax(i,:)*b'));
        end
        
        lambda = eye(n);
        largepi = eye(n);
        for i=1:n
            lambda(i,i) = datay(i)-pi(i);
            largepi(i,i) = pi(i);
        end
        lambda;
        largepi;
        xi = datax*b'+ (largepi*(eye(n)-largepi)) \ (datay-largepi*ones(n,1));
        b = (datax'*largepi*(eye(n)-largepi)*datax \ datax'*largepi*(eye(n)-largepi)*xi)';
        
    end
 
function y = f(x,b)
    y = exp(b(1)+x.*b(2)) ./ (1+exp(b(1)+x.*b(2)));