function B= teach_log_reg(X,Y,R,Lam)
% this function fits binary logistic regresion parameters to data
X(isnan(X))=0;
N=size(X,1);
D=size(X,2);
w=.01*ones(D,1);
v=.01*ones(D,1);
b=0.01;
x0=[w;v;b]; % 2xD +1 parameter vector
%Lambda=0.2;
%x=maxlikelihood(X,Y,R,x0);
x=cg2(x0);
% x_=[-5:1:5];
% y_=[-5:1:5];
% z_=[-5:1:5];
% % x_val=[-3:.2:3];
% % y_val=[-50:5:50];
% %
% % VAL=zeros(length(x_val),length(y_val));
% % for x_=1:length(x_val)
% % for y_=1:length(y_val)
% % VAL(x_,y_)=F(X,Y,R,[x_val(x_),0,y_val(y_)]);
% %
% % VAL(x_,y_)=objFunc([x_val(x_),0,y_val(y_)]);
% %
% % end
% % end
% % mesh(x_val,y_val,VAL')
% % xlabel('x')
% %
% % b=F(X,Y,R,x);
%x = fminsearch(@F,x0)
% x2 = fminsearch(@(x) F(X,Y,R,x),x0);
pause(.1)
% no gradient used
tic
[x2,fval,exitflag] = fminsearch(@(x) objFunc(x),x0,optimset('TolX',1e-7,'MaxFunEvals',1e5,'MaxIter',1e5))
B=prob_dens_f(x);
t(1)=toc
% % with cradient
%options = optimset('GradObj','on','TolX',1e-5,'MaxFunEvals',1e4,'MaxIter',1e4);
[x,fval] = fminunc(@myfun,x0,optimset('GradObj','on','TolX',1e-5,'MaxFunEvals',1e4,'MaxIter',1e4))
% B=prob_dens_f(x);
% t(2)=toc
% x1=x0;
% x1(30)=x1(30)+0.01;
% dx_test=(objFunc(x0)-objFunc(x1))./(x0(30)-x1(30));
% dx_=gradFunc(x0);
function [f,g] = myfun(x)
f = objFunc(x); % Cost function
if nargout > 1
g=gradFunc(x);
end
end
% conjugate gradient minimization
function B=prob_dens_f(x)
D_=size(X,2);
N_=size(X,1);
B=zeros(N_,2);
w_=x(1:D_);
v_=x(D_+1:2*D_);
b_=x(2*D_+1);
tot_summa=0;
for n=1:N_
summa=0;
for d=1:D_
if R(n,d) ~= 0
summa=summa+R(n,d)*(w_(d)*X(n,d)+v_(d));
end
end
B(n,2)=1/(1+exp(-(b_+summa)));
B(n,1)=1-B(n,2);
% tot_summa=tot_summa+log(1+exp(-Y(n)*(b_+summa)));
end
% f=tot_summa;
end
function x= cg2 (x0)
x=x0;
visLineSearch = 0;
% note in rasmussen uses uses polak update with inaxact line search for
% minimize -> http://www.kyb.tuebingen.mpg.de/bs/people/carl/code/minimize/
mode = 'polak'; % 'fletcher', 'hestenes' or 'polak', on 'none'
linesearch = 'test2'; % 'armijo', 'quadratic' or 'exact'
objFuncValue = objFunc(x);
oldObjFuncValue = objFuncValue*0.1+1;
% implementation according to nocedal 5.2 algorithm 5.4 (FR-CG)
% polak and hestenes updates for beta can be found in the same section
dx = gradFunc(x);
dir = -dx/norm(dx);
% iterate
iter = 0;
numOfIter = 1000;
prec = 1e-8;
% convergence if gradient smaller than prec, change in objective function
% smaller than prec or maximum number of iteration reached...
while iter < numOfIter && abs((oldObjFuncValue-objFuncValue)/objFuncValue)>prec && norm(dx)>prec
% iteration counter
iter = iter + 1;
if strcmp(linesearch,'armijo')
alpha = mb_backtrackingLineSearch(objFunc,objFuncValue,x,dx,dir);
elseif strcmp(linesearch,'exact')
maxStepLength = 2;
linesearchFunc = @(delta) objFunc(x+delta*dir);
alpha = fminbnd(linesearchFunc,0,maxStepLength); % matlab inherent function to find 1d min
elseif strcmp(linesearch,'quadratic')
alpha = mb_quadraticApproximationLineSearch(objFunc,objFuncValue,x,dx,dir);
elseif strcmp(linesearch,'test')
alphas=1e-10:(.1*10^(-2))/(2*iter):(.5*10^(-1))/(2*iter);
value = zeros(size(alphas));
for alpha_idx=1:length(alphas)
delta=alphas(alpha_idx);
value(alpha_idx)=objFunc(x+delta*dir);
end
% figure(1)
% plot(alphas,value)
[C,I]=min(value);
alpha=alphas(I);
elseif strcmp(linesearch,'test2')
maxStepLength = (1*10^(-1))/(2*iter);
linesearchFunc = @(delta) objFunc(x+delta*dir);
alpha = fminbnd(linesearchFunc,0,maxStepLength);
if alpha==maxStepLength
disp('warning1')
end
end
% update x
x = x + alpha*dir;
% update obj func values
oldObjFuncValue = objFuncValue;
objFuncValue = objFunc(x);
% update dx
oldDx = dx;
dx = gradFunc(x);
dx=dx./norm(dx);
if strcmp(mode,'fletcher')
beta = (dx'*dx)/(oldDx'*oldDx);
elseif strcmp(mode,'hestenes')
beta = (dx'*(dx-oldDx))/((dx-oldDx)'*dir);
elseif strcmp(mode,'polak')
beta = (dx'*(dx-oldDx))/(oldDx'*oldDx);
elseif strcmp(mode,'none')
beta=0;
end
% update search direction
dir = -dx + beta*dir;
fprintf(['Iteration ' num2str(iter) ' - Obj Func = ' num2str(objFuncValue) ' @ x = [' num2str(x') ']\n']);
end
end
% function [x] = maxlikelihood(X,Y,R,x)
%
% k=1;
% MINVAL=0;
% MAX_ITER=1000;
% alpha=.1;
% values=zeros(MAX_ITER,1);
%
% d_k=zeros(size(x));
%
% while abs(F(X,Y,R,x))>MINVAL && k<MAX_ITER
%
%
%
% if k>1
% g_k=G(X,Y,R,x);
% beta=(g_k'*g_k)/(g_k1'*g_k1);
% g_k1=g_k;
%
% d_k=-G(X,Y,R,x)+beta*d_k;
% else % first itertion
% d_k=-G(X,Y,R,x);
% g_k1=-d_k;
% g_k=-d_k;
% beta=0;
% end
%
%
%
% [x,alpha]=linesearch(x,d_k,X,Y,R);
%
% % x=x+alpha*d_k;
%
% values(k)= F(X,Y,R,x);
% k=k+1;
% end
%
% figure(1000)
% plot(values)
% pause(1)
%
% end
%
% function [x_best,best_alpha]=linesearch(x,d_k,X,Y,R)
% min_val=1e50;
% for alpha=.1:.1:10
% x=x+alpha*d_k;
% value= F(X,Y,R,x);
% if value < min_val || alpha == .1
% best_alpha=alpha;
% x_best=x;
% end
% end
% end
% objective function
% function [f]=F(X,Y,R,x)
% D_=size(X,2);
% N_=size(X,1);
% w_=x(1:D_);
% v_=x(D_+1:2*D_);
% b_=x(2*D_+1);
% Lambda=0;
%
% tot_summa=0;
% for n=1:N_
%
% summa=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa=summa+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
%
% tot_summa=tot_summa+log(1+exp(-Y(n)*(b_+summa)))+Lambda*(w_'*w_+v_'*v_);
%
% end
%
% f=tot_summa;
%
% end
% objective function gradient
% function [g]=G(X,Y,R,x)
% D_=size(X,2);
% N_=size(X,1);
% w_=x(1:D_);
% v_=x(D_+1:2*D_);
% b_=x(2*D_+1);
% g=zeros(size(x));
%
% Lambda=0.1;
%
% sum_w=sum(w_);
% sum_v=sum(v_);
%
% for g_idx=1:length(x)
% if g_idx<=length(w_)
%
% tot_summa=0;
% for n=1:N_
%
% summa1=0;
% summa2=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa1=summa1+R(n,d)*X(n,d);
% summa2=summa2+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
%
% tot_summa=tot_summa+(-Y(n)*summa1*exp(-Y(n)*(b_+summa2)))/(1+exp(-Y(n)*(b_+summa2)));
%
% end
% g(g_idx)=tot_summa+2*sum_w;
%
% elseif g_idx <= 2*Lambda*length(w_)
%
% tot_summa=0;
% for n=1:N_
%
% summa1=0;
% summa2=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa1=summa1+R(n,d);
% summa2=summa2+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
%
% tot_summa=tot_summa+(-Y(n)*summa1*exp(-Y(n)*(b_+summa2)))/(1+exp(-Y(n)*(b_+summa2)));
%
% end
% g(g_idx)=tot_summa+2*sum_v;
% else
%
% tot_summa=0;
% for n=1:N_
%
%
% summa2=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa2=summa2+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
%
% tot_summa=tot_summa+(-Y(n)*exp(-Y(n)*(b_+summa2)))/(1+exp(-Y(n)*(b_+summa2)));
%
% end
%
% g(g_idx)=tot_summa;
% end
% end
%
%
%
% end
function f=objFunc(x)
D_=size(X,2);
N_=size(X,1);
w_=x(1:D_);
v_=x(D_+1:2*D_);
b_=x(2*D_+1);
% Lambda=0.1;
tot_summa=0;
for n=1:N_
% summa=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa=summa+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
summa=R(n,:)*(X(n,:)'.*w_+v_);
tot_summa=tot_summa+log(1+exp(-Y(n)*(b_+summa)));
end
f=tot_summa+Lam*(w_'*w_+v_'*v_);
%f=tot_summa+Lam*N_*(norm(w_,1)+norm(v_,1));
end
function g=gradFunc(x)
D_=size(X,2);
N_=size(X,1);
w_=x(1:D_);
v_=x(D_+1:2*D_);
b_=x(2*D_+1);
g=zeros(size(x));
% Lambda=0.1;
sum_w=sum(w_);
sum_v=sum(v_);
for g_idx=1:length(x)
if g_idx<=length(w_)
tot_summa=0;
for n=1:N_
% summa1=0;
% summa2=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa1=summa1+R(n,d)*X(n,d);
% summa2=summa2+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
summa1=R(n,g_idx)*X(n,g_idx)';
summa2=R(n,:)*(X(n,:)'.*w_+v_);
tot_summa=tot_summa+(-Y(n)*summa1*exp(-Y(n)*(b_+summa2)))/(1+exp(-Y(n)*(b_+summa2)));
end
g(g_idx)=tot_summa+2*Lam*x(g_idx);
elseif g_idx <= 2*length(w_)
tot_summa=0;
for n=1:N_
% summa1=0;
% summa2=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa1=summa1+R(n,d);
% summa2=summa2+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
summa1=R(n,g_idx-length(w_))*1;
summa2=R(n,:)*(X(n,:)'.*w_+v_);
tot_summa=tot_summa+(-Y(n)*summa1*exp(-Y(n)*(b_+summa2)))/(1+exp(-Y(n)*(b_+summa2)));
end
g(g_idx)=tot_summa+2*Lam*x(g_idx);
else
tot_summa=0;
for n=1:N_
% summa2=0;
% for d=1:D_
% if R(n,d) ~= 0
% summa2=summa2+R(n,d)*(w_(d)*X(n,d)+v_(d));
% end
% end
summa2=R(n,:)*(X(n,:)'.*w_+v_);
tot_summa=tot_summa+(-Y(n)*exp(-Y(n)*(b_+summa2)))/(1+exp(-Y(n)*(b_+summa2)));
end
g(g_idx)=tot_summa;
end
end
end
end