#epsilon := 0.1: c := 1.0: a := 1.0: ##### center the data center_xy := proc(xs,ys,n) local i,xbar,ybar,xc,yc; xc := Vector(n); yc := Vector(n); xbar := 0.; ybar := 0.; for i from 1 to n do xbar := xbar + xs[i]; ybar := ybar + ys[i]; end do; xbar := xbar/n; ybar := ybar/n; for i from 1 to n do xc[i] := xs[i] - xbar; yc[i] := ys[i] - ybar; end do; RETURN([xc,yc,xbar,ybar]) end: SVRegQP := proc(x,y,n,epsilon,c,a) local i, j, C, K,K1, H, Aeq, beq, bl,bu,sol; # create C for i from 1 to n do C[i] := y[i] - epsilon: C[n+i] := -y[i] - epsilon: end do: C := Vector(2*n,C,datatype=float): # create K K := Matrix(2*n,2*n); for i from 1 to n do for j from 1 to n do K[i,j] := -k(x[i],x[j],a): K[n+i,j] := -K[i,j]: K[i,n+j] := -K[i,j]: K[n+i,n+j] := K[i,j]: end do: end do: H := Matrix(2*n,2*n,K,datatype = float): # create Aeq for i from 1 to n do Aeq[i] := 1: Aeq[n+i] := -1: end do: Aeq := Matrix(1,2*n,Vector[row](2*n,Aeq,datatype = float)): # create beq, bl and bu beq := Vector([0]): bl := ConstantVector(0,2*n): bu := ConstantVector(c,2*n): sol := Optimization[QPSolve]([C,H],[NoUserValue,NoUserValue,Aeq,beq],[0,c],maximize); K1 := LinearAlgebra[SubMatrix](K,n+1..2*n,1..n); RETURN(sol,K1); end: ### compute dlambda and plambda get_dlplambda := proc(lambda,n) local i,dl,pl; dl := Vector(n); pl := Vector(n); for i from 1 to n do dl[i] := lambda[i]-lambda[n+i]: pl[i] := lambda[i]+lambda[n+i]: end do: RETURN(dl,pl); end: ### optimal function fopt := proc(x,dlambda,xs,n,a,tol) local sm,i; sm := 0.; for i from 1 to n do if (abs(dlambda[i]) > tol) then sm := sm + dlambda[i]*k(xs[i],x,a); end if; end do; RETURN(sm); end: ### compute the indices of the support vectors. get_svs := proc(dlambda,n,tol) local i,j,sv; j := 0; for i from 1 to n do if (abs(dlambda[i]) > tol) then j := j+1; sv[j] := i; end if; end do; RETURN([j,Vector(j,sv)]); end: get_w := proc(dlambda,svs,xs) local i,w; w := 0.; for i in svs do w := w + dlambda[i]*xs[i]; end do; RETURN(w); end: #### compute b get_b := proc(Y,lambda,n,K) local i,dlpl,b,dl,pl; dlpl := get_dlplambda(lambda,n); dl := dlpl[1]; pl := dlpl[2]; b := (LinearAlgebra[DotProduct](pl,Y)- LinearAlgebra[BilinearForm](pl,dl,K))/add(pl[i],i=1..n); RETURN([b,dl,pl]); end: ### sample n pts. with error from a function get_sample := proc(n,f,a,b,sigma) local xs,ys,errors,i,pts,ppts; xs := Vector(n,[stats[random, uniform[a,b]](n)]); ys := Vector(n,[seq(evalf(f(xs[i])),i=1..n)]); errors := Vector(n,[stats[random, normald[0,sigma]](n)]); ys := ys+errors; i := 'i'; pts := [['xs[i]', 'ys[i]'] \$i=1..n]; ppts := plot(pts,x=a..b, style=point,symbol=circle): RETURN([xs,ys,ppts]); end: #### Get the Regression and a plot. SVRegression := proc(xs,ys,n,k,kas,epsilon,c,tol,lower,upper) ## Input: ## xs,ys the data ## n number of data points ## k The kernel function k(x,x',params) ## kas list of kernel parameters ## epsilon epsilon insensitive cost ## c smoothing parameter: c*(emp.risk) + |w|^2/2 ## tol tolerance for supp. vecs. |lambda_i| > tol. ## lower,upper limits for plotting ## ## Output: ## b,dlambda b+fopt(x,dlambda...) is the estimate ## aplot plot of the estimate and data ########################################################### local K,lambdas,Sol,dlambda,plambda,b,aplot,bbs,x; Sol := SVRegQPsolve(xs,ys,n,epsilon,c,k,kas); K := Sol[2]; lambdas := Sol[1][2]; bbs := get_b(ys,lambdas,n,K); b := bbs[1]; dlambda := bbs[2]; plambda := bbs[3]; aplot := plot( b + fopt1(x,dlambda,xs,n,k,kas,tol), x = lower..upper,color=green): RETURN([b,dlambda,aplot]); end: SVRegQPsolve := proc(x,y,n,epsilon,c,k,kas) local i, j, C, K,K1, H, Aeq, beq, bl,bu,sol; # create C for i from 1 to n do C[i] := y[i] - epsilon: C[n+i] := -y[i] - epsilon: end do: C := Vector(2*n,C,datatype=float): # create K K := Matrix(2*n,2*n); for i from 1 to n do for j from 1 to n do K[i,j] := -k(x[i],x[j],op(kas)): K[n+i,j] := -K[i,j]: K[i,n+j] := -K[i,j]: K[n+i,n+j] := K[i,j]: end do: end do: H := Matrix(2*n,2*n,K,datatype = float): # create Aeq for i from 1 to n do Aeq[i] := 1: Aeq[n+i] := -1: end do: Aeq := Matrix(1,2*n,Vector[row](2*n,Aeq,datatype = float)): # create beq, bl and bu beq := Vector([0]): bl := ConstantVector(0,2*n): bu := ConstantVector(c,2*n): sol := Optimization[QPSolve]([C,H],[NoUserValue,NoUserValue,Aeq,beq],[0,c],maximize); K1 := LinearAlgebra[SubMatrix](K,n+1..2*n,1..n); RETURN(sol,K1); end: # define a kernel k := (x,y,a) -> exp(-(x-y)^2/a^2): ### optimal function fopt1 := proc(x,dlambda,xs,n,k,kas,tol) local sm,i; sm := 0.; for i from 1 to n do if (abs(dlambda[i]) > tol) then sm := sm + dlambda[i]*k(xs[i],x,op(kas)); end if; end do; RETURN(sm); end: ### plot the sv points plot_svs := proc(SVs,XYp,a,b) local nsv,px,py,i,svpts, ap1,aplot,v; nsv := SVs[1]; v := SVs[2]; px := XYp[1]; py := XYp[2]; i := 'i'; svpts := [[ 'px[i]', 'py[i]' ] \$i=1..nsv]; ap1 := plot(svpts,x=a..b,style=point,symbol=cross,color=black): RETURN({ap1,XYp[3]}); end: ##### plot the tube plot_tube := proc(SVr,XYp,epsilon,n,k,kas,tol,a,b) local beta, bplus, bminus, fhat,x,p1,p2,SVs; beta := SVr[1]; bplus := beta+epsilon; bminus := beta-epsilon; fhat := x -> fopt1(x,SVr[2],XYp[1],n,k,kas,tol); p1 := plot({bplus+fhat(x),bminus+fhat(x)},x = a..b,color=blue): SVs := get_svs(SVr[2],n,tol); p2 := op(plot_svs(SVs,XYp,a,b)); RETURN({p1,p2}); end: