Nothing
# Revised Version of Fitzenberger's Steepest Descent Algorithm for the Powell Estimator
#
# Roger Koenker: last revision: 6 February 2008
#
# Problem: || min{Ax,c} - b ||_tau = min!
#
subroutine powell(n,p,p2,A,b,c,x,tau,h,f,U,s,g,d,xh,maxit,nflag)
integer n,p,p2
double precision x(p),A(n,p),b(n),c(n)
double precision f(n),U(p,p),s(n),g(p2),d(p),xh(p)
double precision zero, one, mone, step, tau,pow
integer h(p),hin,hout,k,it,inset,maxit,nflag
PARAMETER(zero = 0.0d0, one = 1.d0, mone = -1.d0)
it = 0
repeat {
it = it + 1
if(it > 1) # Otherwise, assume U is A(h)^{-1}
call pivot(n,p,h,hin,hout,A,U,d,xh,nflag)
if(nflag > 0)
{nflag = nflag + 2; return}
do i = 1,p{
xh(i) = b(h(i))
}
call dgemv('N',p,p,one,U,p,xh,1,zero,x,1)
call dgemv('N',n,p,one,A,n,x,1,zero,f,1)
do i = 1,n{
if(inset(p,i,h) > 0 | f(i) > c(i))
s(i) = zero
else if(b(i) < f(i))
s(i) = one - tau
else
s(i) = - tau
}
call dgemv('T',n,p,one,A,n,s,1,zero,xh,1)
call dgemv('T',p,p,one,U,p,xh,1,zero,g,1)
do i = 1,p {
if(f(h(i)) < c(h(i)))
if(b(h(i)) < c(h(i)))
g(i + p) = - g(i) + one - tau
else
g(i + p) = - g(i) - tau
else
g(i + p) = - g(i) + tau
g(i) = g(i) + one - tau
}
k = idmin(p2,g,1)
if(g(k) >= 0 | it > maxit)
break
call dscal(p,zero,d,1)
if(k <= p)
call daxpy(p,one,U(1,k),1,d,1)
else{
k = k - p
call daxpy(p,mone,U(1,k),1,d,1)
}
call dgemv('N',n,p,one,A,n,d,1,zero,s,1)
do i = 1,n {
call dcopy(p,x,1,xh,1)
step = (b(i) - f(i))/s(i)
call daxpy(p,step,d,1,xh,1)
s(i) = pow(n,p,xh,A,b,c,tau)
}
hin = idmin(n,s,1)
if(inset(p,hin,h) > 0)
{nflag = 2; break}
hout = h(k)
}
if(it > maxit) nflag = 1
return
end
#####################################################
subroutine pivot(n,p,h,hin,hout,A,B,u,v,eflag)
integer n,p,h(p),hin,hout,inset,k,eflag
double precision A(n,p),B(p,p),u(p),v(p)
double precision zero,one
PARAMETER(zero = 0.d0, one = 1.d0)
eflag = 0
k = inset(p,hout,h)
if(k == 0)
{eflag = 1; return}
if(inset(p,hin,h) > 0)
{eflag = 2; return}
if(hin < 1 | hin > n)
{eflag = 3; return}
call dcopy(p,A(hin,1),n,v,1)
call dgemv('T',p,p,one,B,p,v,1,zero,u,1)
call dcopy(p,B(1,k),1,v,1)
do j = 1,p{
do i = 1,p{
if(j == k)
B(i,j) = B(i,j)/u(k)
else
B(i,j) = B(i,j) - (u(j)/u(k)) * v(i)
}
}
h(k) = hin
return
end
#####################################################
integer function inset(p,k,h)
integer p,k,h(p)
do inset = 1,p{
if(h(inset) == k) return
}
inset = 0
return
end
#####################################################
double precision function pow(n,p,x,A,b,c,tau)
integer n,p
double precision x(p),A(n,p),b(n),c(n)
double precision tau,u,zero,rho,fit,ddot
PARAMETER(zero= 0.d0)
pow = zero
do i = 1,n{
fit = ddot(p,A(i,1),n,x,1)
u = b(i) - min(fit,c(i))
pow = pow + rho(u, tau)
}
return
end
#####################################################
double precision function rho(u,tau)
double precision u,tau,one
PARAMETER(one = 1.d0)
if(u < 0)
rho = u * (tau - one)
else
rho = u * tau
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.