auc.b <- function(x=dat.sin,xt=dat.sin,label1=TRUE,T=5,lambda=0.001){ N=T probs=1:30/30 y=x[,dim(x)[2]]==label1 # training data xのラベル y=Tがpresent(kyphosis(後彎症)あり),y=Fがabsent x=x[,1:2] n0=sum(y==F) # y=Fのサンプル数 n1=sum(y==TRUE) #y=Tのサンプル数   n = dim(x)[1]  #n=n0+n1 m = dim(x)[2] #マーカーの数(3つ) sd.n=1 T.auc=0 F0=rep(0,n0) F1=rep(0,n1) B=list() Z=list() Fz=list() for(i in 1:m){ B[[i]]=quantile(x[,i],probs=probs) Z[[i]]=c(min(x[,i]),min(x[,i]),B[[i]],max(x[,i])) Fz[[i]]=rep(0,length(Z[[i]])) } nb=unlist(lapply(B,length)) M=NULL for(i in 1:m){ Md=matrix(x[,i],n,nb[i]) Md=t(t(Md)>=B[[i]]) M=cbind(M,Md) } M0=M[y==F,] M1=M[y==TRUE,] a=numeric(0) b=numeric(0) p=numeric(0) alpha=numeric(0) ff=function(x){ nb=length(x)-3 2*(x[1:nb]-3*x[2:(nb+1)]+3*x[3:(nb+2)]-x[4:(nb+3)]) } nn=100 eps=0.1 grid1=seq(min(x[,1])-eps,max(x[,1])+eps,length=nn) grid2=seq(min(x[,2])-eps,max(x[,2])+eps,length=nn) Grid=expand.grid(grid1,grid2) windows(width=12,height=8) xx=-100:800/100 for(k in 1:N){#N:ブースティングの繰り返し数 T.auc=T.auc+1 w1=unlist(lapply(Fz,ff)) d=0 for(i in 1:n0) # 弱判別機f_tの選択のための計算(p969のstep a.) d=d+apply(1/(n0*n1)*dnorm(F1-F0[i],sd=sd.n)*(t(t(M1)-M0[i,])),2,sum) d=d-lambda*w1 maxd=d[order(d)[length(d)]] mind=d[order(d)[1]] if(maxd>=-mind){ a[k]=1 l=order(d)[length(d)] } else{ a[k]=-1 l=order(d)[1] } ii=0 ld=l while(ld>0){ ii=ii+1 ld=ld-nb[ii] } p[k]=ii #選ばれたf_tに対応するマーカーのid ld=ld+nb[ii] b[k]=B[[p[k]]][ld] #選ばれたf_t(stump) if(k<4) al=1/sd.n else al=median(alpha[(k-3):(k-1)])/sd.n al0=0 nz=length(Fz[[p[k]]]) AA=0 while(abs(al-al0)>0.000001&AA<5){#弱判別機の係数\alphaの更新 AA=AA+1 al0=al d1=0 d2=0 for(i in 1:n0){ d1=d1+sum(1/(n0*n1)*dnorm(F1-F0[i]+al0*a[k]*(M1[,l]-M0[i,l]),sd=sd.n)*a[k]*(M1[,l]-M0[i,l])) d2=d2-sum(1/(n0*n1)*(F1-F0[i]+al0*a[k]*(M1[,l]-M0[i,l]))*dnorm(F1-F0[i]+al0*a[k]*(M1[,l]-M0[i,l]),sd=sd.n)/sd.n^2*(M1[,l]-M0[i,l])^2) } dw1=2*(a[k]*(Fz[[p[k]]][ld]-3*Fz[[p[k]]][ld+1]+3*Fz[[p[k]]][ld+2] -Fz[[p[k]]][ld+3])+2*al0) dw2=4 al=al0-(d1-lambda*dw1)/(d2-lambda*dw2) } alpha[k]=al if(a[k]==1){ F0=F0+al*M0[,l] F1=F1+al*M1[,l] Fz[[p[k]]][(ld+2):nz]=Fz[[p[k]]][(ld+2):nz]+al } else{ F0=F0+al*(!M0[,l]) F1=F1+al*(!M1[,l]) Fz[[p[k]]][1:(ld+1)]=Fz[[p[k]]][1:(ld+1)]+al } max.fz=0 for(i in 1:m){ if(max.fz