function [lcpred, hylc, varlc,nsing]=lcfun(x,y,n,h,mpred,xpred, vary,epsilon) % give the predicted value using local cubic regression nsing=zeros(1,mpred); varlc=zeros(1,mpred); hylc=zeros(1,mpred); for ii=1:mpred w1=zeros(1,n); s0=0; s1=0; s2=0; s3=0; s4=0; s5=0; s6=0; t0=0; t1=0; t2=0; t3=0; for j=1:n tmp=x(j)-xpred(ii); arg=tmp/h; if abs(arg)<1 w1(j)=0.75*(1-arg^2)/h; %nhd=nhd+1; end s0=s0+w1(j); s1=s1+w1(j)*tmp; s2=s2+w1(j)*tmp^2; s3=s3+w1(j)*tmp^3; s4=s4+w1(j)*tmp^4; s5=s5+w1(j)*tmp^5; s6=s6+w1(j)*tmp^6; t0=t0+w1(j)*y(j); t1=t1+w1(j)*y(j)*tmp; t2=t2+w1(j)*y(j)*tmp^2; t3=t3+w1(j)*y(j)*tmp^3; end b0=-s4*s4*s4+ 2*s3*s4*s5 - s2*s5*s5 - s3*s3*s6 + s2*s4*s6; b1=s3*s4*s4 - s3*s3*s5 - s2*s4*s5 + s1*s5*s5 + s2*s3*s6 - s1*s4*s6; b2=-s3*s3*s4 + s2*s4*s4 + s2*s3*s5 - s1*s4*s5 - s2*s2*s6 + s1*s3*s6; b3=s3*s3*s3 - 2*s2*s3*s4 + s1*s4*s4 + s2*s2*s5 - s1*s3*s5; temp1=s3*s3*s3*s3- 3*s2*s3*s3*s4 + s2*s2*s4*s4 + 2*s1*s3*s4*s4 - s0*s4*s4*s4 + 2*s2*s2*s3*s5 - 2*s1*s3*s3*s5 - 2*s1*s2*s4*s5; temp2=2*s0*s3*s4*s5 + s1*s1*s5*s5 - s0*s2*s5*s5 -s2*s2*s2*s6 + 2*s1*s2*s3*s6 - s0*s3*s3*s6 - s1*s1*s4*s6 + s0*s2*s4*s6; delta=temp1+temp2; if delta==0 delta=epsilon; nsing(ii)=1; end for j=1:n tmp=x(j)-xpred(ii); tp0=b0*w1(j); tp1=b1*w1(j)*tmp; tp2=b2*w1(j)*tmp^2; tp3=b3*w1(j)*tmp^3; varlc(ii)=varlc(ii)+(1/delta*(tp0+tp1+tp2+tp3))*(1.0/delta*(tp0+tp1+tp2+tp3))*vary; hylc(ii)=hylc(ii)+(1.0/delta*(tp0+tp1+tp2+tp3))*y(j); end lcpred(ii)= (b0*t0+b1*t1+b2*t2+b3*t3)/delta; end