partials.rar_Partials_least squares

function[beta,phi,s,llMax]=pls(data,m,beta,ind,j,retninger) m; convcrit=0.00001; [d X S A n_deaths]=var_maker(data); [n p]=size(X); if m==0 beta=zeros(p,1); phi=logbreslow(X*beta,A,n_deaths); my=exp(X*beta).*(A*exp(phi)); s=1; ll=d'*log(my)-sum(my); else eta=zeros(n,1); phi=logbreslow(eta,A,n_deaths); my=exp(eta).*(A*exp(phi)); myny=my; betas=[]; T=zeros(n,1); P=zeros(p,1); q=zeros(1,1); W=zeros(p,1); s=0; myDiffs=[]; ll=[-100000 d'*log(my)-sum(my)]; c=0; while ((max(abs(my-myny))>convcrit | s==0) & s<10000) %[m s] [mymaxdiff whichmy]=max(abs(my-myny)); s=s+1; my=myny; V=diag(my); myDiffs=[myDiffs mymaxdiff]; en=ones(n,1); xhat=(X'*my)/sum(my); E=X-en*xhat'; f=eta+d./my-en; for k=1:m, wk=E'*(my.*f); tk=E*wk; qk=f'*(my.*tk)/(tk'*(my.*tk)); pk=E'*(my.*tk)/(tk'*(my.*tk)); E=E-tk*pk'; f=f-tk*qk; q(k,1)=qk; T(:,k)=tk; W(:,k)=wk; P(:,k)=pk; end %Linjes�k: if c==0 llc=[]; cvec=linspace(0,1,10); for c=cvec etac=c*T*q+(1-c)*eta; etac=etac-mean(etac); phic=logbreslow(etac,A,n_deaths); mynyc=exp(etac).*(A*exp(phic)); llc=[llc d'*log(mynyc)-sum(mynyc)]; end [lmax cind]=max(llc); cmax=cvec(cind); %Pr�ver et finere grid dersom det f�rste ikke %gav et positivt resultat: if(cind==1) llc=[]; cvec=linspace(0,cvec(2),50); for c=cvec etac=c*T*q+(1-c)*eta; etac=etac-mean(etac); phic=logbreslow(etac,A,n_deaths); mynyc=exp(etac).*(A*exp(phic)); llc=[llc d'*log(mynyc)-sum(mynyc)]; end [lmax cind]=max(llc); cmax=cvec(cind); if(cind==1) % figure % plot(cvec,llc) % title(['Log likelihood iter ',num2str(s)]) % xlabel('c_m') cmax; ll=[ll lmax]; break end end % figure % plot(cvec,llc) % title(['Log likelihood iter ',num2str(s)]) % xlabel('c_m') ll=[ll lmax]; c=0; else cmax=c; end eta=cmax*T*q+(1-cmax)*eta; eta=eta-mean(eta); phi=logbreslow(eta,A,n_deaths); myny=exp(eta).*(A*exp(phi)); end beta=W*inv(P'*W)*q; %beta=(W*inv(P'*W)*q)./sqrt(var(data(:,3:p+2)))'; % figure % plot(ll) % title('Log likelihood iters ') % xlabel('iter no.') llMax=ll(end); end;