Nothing
subroutine baklo(x,y,w,npetc,wddnfl,spatol,match,
etal,s,eta,beta,var,dof,
qr,qraux,qpivot,effect,iv,v,iwork,work)
implicit double precision(a-h,o-z)
integer n,p,q,nit,maxit,qrank
integer npetc(7),wddnfl(*),match(*),qpivot(*),iv(*),iwork(*)
##integer which(q),dwhich(q),degree(q),nef(q),liv(q),lv(q),nvmax(q)
double precision x(*),y(*),w(*),spatol(*),
etal(*),s(*),eta(*),beta(*),var(*),dof(*),
qr(*),qraux(*),v(*),effect(*),work(*)
#work size: 4*n + sum( nef(k)*(pj+dj+4)+5+3*pj ) +5*n
# = 9*n + sum( nef(k)*(pj+dj+4)+5+3*pj )
n=npetc(1)
p=npetc(2)
q=npetc(3)
maxit=npetc(5)
qrank=npetc(6)
call baklo0(x,n,p,y,w,q,wddnfl(1),wddnfl(q+1),wddnfl(2*q+1),
spatol(1),wddnfl(3*q+1),dof,match,wddnfl(4*q+1),
etal,s,eta,beta,var,spatol(q+1),
nit,maxit,qr,qraux,qrank,qpivot,effect,
work(1),work(n+1),work(2*n+1),work(3*n+1),
iv,wddnfl(5*q+1),wddnfl(6*q+1),v,wddnfl(7*q+1),
iwork(1),work(4*n+1))
npetc(4)=nit
npetc(6)=qrank
return
end
subroutine baklo0(x,n,p,y,w,q,which,dwhich,pwhich,span,degree,dof,match,nef,
etal,s,eta,beta,var,tol,nit,maxit,
qr,qraux,qrank,qpivot,effect,z,old,sqwt,sqwti,
iv,liv,lv,v,nvmax,iwork,work)
implicit double precision(a-h,o-z)
integer n,p,q,which(q),dwhich(q),pwhich(q),degree(q),match(n,q),nef(q),nit,
maxit,qrank,qpivot(p),iv(*),liv(q),lv(q),nvmax(q),iwork(q)
double precision x(n,p),y(n),w(n),span(q),dof(q),
etal(n),s(n,q),eta(n),beta(p),var(n,q),tol,
qr(n,p),qraux(p),v(*),effect(n),work(*)
#work should be sum( nef(k)*(pj+dj+4)+5+3*pj ) +5*n
double precision z(*),old(*),dwrss,ratio
double precision sqwt(n),sqwti(n)
logical anyzwt
double precision deltaf, normf,onedm7
integer job,info,slv,sliv,iw,j,dj,pj
onedm7=1d-7
job=1101;info=1
if(q==0)maxit=1
ratio=1d0
# fix up sqy's for weighted problems.
anyzwt=.false.
do i=1,n{
if(w(i)>0d0){
sqwt(i)=dsqrt(w(i))
sqwti(i)=1d0/sqwt(i)
}
else{
sqwt(i)=0d0
sqwti(i)=0d0
anyzwt=.true.
}
}
# if qrank > 0 then qr etc contain the qr decomposition
# else baklo computes it.
if(qrank==0){
do i=1,n{
do j=1,p{
qr(i,j)=x(i,j)*sqwt(i)
}
}
do j=1,p{qpivot(j)=j}
call dqrdca(qr,n,n,p,qraux,qpivot,work,qrank,onedm7)
}
do i=1,n{
eta(i)=0d0
for(j=1;j<=q;j=j+1){
eta(i)=eta(i)+s(i,j)
}
}
nit=0
while ((ratio > tol )&(nit < maxit)){
# first the linear fit
deltaf=0d0
nit=nit+1
do i=1,n{
z(i)=(y(i)-eta(i))*sqwt(i)
old(i)=etal(i)
}
# call dqrsl1(qr,dq,qraux,qrank,sqz,one,work(1),etal,two,three)
#job=1101 -- computes fits, effects and beta
call dqrsl(qr,n,n,qrank,qraux,z,work(1),effect(1),beta,
work(1),etal,job,info)
# now unsqrt the fits
#Note: we dont have to fix up the zero weights till the end, since their fits
#are always immaterial to the computation
do i=1,n{
etal(i)=etal(i)*sqwti(i)
}
# now a single non-linear backfitting loop
sliv=1
slv=1
iw=5*n+1
for(k=1;k<=q;k=k+1){
j=which(k)
dj=dwhich(k)
pj=pwhich(k)
do i=1,n{
old(i)=s(i,k)
z(i)=y(i)-etal(i)-eta(i)+old(i)
}
call lo1(x(1,j),z,w,n,dj,pj,nvmax(k),span(k),degree(k),match(1,k),
nef(k),nit,dof(k),s(1,k),var(1,k),work(iw),
# xin,win
work(iw+pj+1),work(iw+nef(k)*dj+pj+1),
# sqwin,sqwini,
work(iw+nef(k)*(dj+1)+pj+2),work(iw + nef(k)*(dj+2)+pj+2),
# xqr,qrank,
work(iw+nef(k)*(dj+3)+pj+2),work(iw+nef(k)*(pj+dj+4)+pj+2),
# qpivot,qraux,
# work(iw+nef(k)*(pj+dj+4)+pj+3),work(iw+nef(k)*(pj+dj+4)+4+2*pj),
iwork(1),work(iw+nef(k)*(pj+dj+4)+4+2*pj),
iv(sliv),liv(k),lv(k),v(slv),
work(1) )
#work should be sum( nef(k)*(pj+dj+4)+5+3*pj ) +5*n
# In the call above I give lo1 pieces of work to use for storing
# the qr decomposition, and it gets the same undisturbed portion
# each time. The fact that it is given a double work word for qrank
# is irrelevant but convenient; it still stores the integer qrank there.
# I do this because there is a partition like this for each lo() term
# in the model, and the number of them is variable
sliv=sliv+liv(k)
slv=slv+lv(k)
iw=iw+nef(k)*(pj+dj+4)+5+3*pj
do i=1,n{
eta(i)=eta(i)+s(i,k)-old(i)
}
deltaf=deltaf+dwrss(n,old,s(1,k),w)
}
normf=0d0
do i=1,n{
normf=normf+w(i)*eta(i)*eta(i)
}
if(normf>0d0){
ratio=dsqrt(deltaf/normf)
}
else {ratio = 0d0}
}
#now package up the results
do j=1,p {work(j)=beta(j)}
do j=1,p {beta(qpivot(j))=work(j)}
if(anyzwt){
do i=1,n {
if(w(i) <= 0d0){
etal(i)=0d0
do j=1,p{
etal(i)=etal(i)+beta(j)*x(i,j)
}
}
}
}
do i=1,n
eta(i)=eta(i)+etal(i)
return
end
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.