Kernel estimation

1. Density estimation

1.1. Histogram

1.1.1. Real world data

1.1.1.1. Galaxies

library('sm')

Package 'sm', version 2.2-5.7: type help(sm) for summary information

library(MASS)

Attachement du package : 'MASS'

L'objet suivant est masqué depuis 'package:sm':

    muscle

data(galaxies)
?galaxies
hist(galaxies)

hist(galaxies,freq=F)

hist(galaxies,freq=F,nclass=20)

hist(galaxies,breaks=quantile(galaxies,seq(0,1,len=20)))

1.1.1.2. Faithful

data(faithful)
attach(faithful)
?faithful
hist(waiting,freq=F)

hist(waiting,freq=F,nclass=20)

hist(waiting,breaks=quantile(waiting,seq(0,1,len=20)))

hist(eruptions,freq=F)

hist(eruptions,freq=F,nclass=20)

hist(eruptions,breaks=quantile(eruptions,seq(0,1,len=30)))

1.1.2. Simulation of a mixing model

rmixing=function(n,alpha,l0,l1,p0,p1)
# Generate data from a mixing model 
{
  z=rbinom(n,1,alpha)
  f1=eval(parse(text=paste('r',l1,'(',paste(c(n,p1),collapse=','),')',sep='')))
  f0=eval(parse(text=paste('r',l0,'(',paste(c(n,p0),collapse=','),')',sep='')))
  x=z*f1+(1-z)*f0
  return(x=x)
}

dmixing=function(t,alpha,l0,l1,p0,p1)
# draw the density of the mixing model
{
  res=alpha*eval(parse(text=paste('d',l1,'(t,',paste(p1,collapse=','),')',sep='')))+(1-alpha)*eval(parse(text=paste('d',l0,'(t,',paste(p0,collapse=','),')',sep='')))
}

#Example  
n=300
alpha=0.3
l0='norm'
p0=c(8,1)
l1='norm'
p1=c(0,2)
s=seq(-10,10,0.001)

x=rmixing(n,alpha,l0,l1,p0,p1)

#### histogram
par(mfrow=c(1,3))
hist(x,freq=F,ylim=c(0,0.4))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
hist(x,freq=F,ylim=c(0,0.4),nclass=20)
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
hist(x,breaks=quantile(x,seq(0,1,len=20)),ylim=c(0,0.4))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')

1.2. Moving window estimator

1.2.1. On the simulated mixing model

par(mfrow=c(2,3))
plot(density(x,bw=0.001,kernel='rectangular'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=0.01,kernel='rectangular'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=0.1,kernel='rectangular'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,kernel='rectangular'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=10,kernel='rectangular'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,kernel='rectangular',bw=100),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')

par(mfrow=c(1,1))
hist(x,freq=F,ylim=c(0,0.4),xlim=c(-7,12))
lines(density(x,kernel='rectangular'),col='blue')
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')

1.2.2. Come back to real world data

# Galaxies 
hist(galaxies,freq=F,ylim=range(density(galaxies,kernel='rectangular')$y))
lines(density(galaxies,kernel='rectangular'),col='blue')

# Faithful 
hist(waiting,freq=F,ylim=range(density(waiting,kernel='rectangular')$y))
lines(density(waiting,kernel='rectangular'),col='blue')

hist(eruptions,freq=F,ylim=range(density(eruptions,kernel='rectangular')$y))
lines(density(eruptions,kernel='rectangular'),col='blue')

1.3. Kernel estimator

1.3.1. On the simulated mixing model

Effect of h value

par(mfrow=c(2,3))
plot(density(x,bw=0.001,kernel='g'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=0.01,kernel='g'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=0.1,kernel='g'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,kernel='g'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=10,kernel='g'),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,kernel='g',bw=100),ylim=c(0,0.4),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')

par(mfrow=c(1,1))
hist(x,freq=F,ylim=c(0,0.4),xlim=c(-7,12))
lines(density(x,kernel='g'),col='blue')
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')

Effect of K value

par(mfrow=c(2,3))
plot(density(x,bw=1,kernel='r'),main='Uniform',ylim=c(0,0.3),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=1,kernel='g'),main='Gaussian',ylim=c(0,0.3),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=1,kernel='e'),main='Epanechnikov',ylim=c(0,0.3),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=1,kernel='triangular'),main='Triangular',ylim=c(0,0.3),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,bw=1,kernel='b'),main='Biweight',ylim=c(0,0.3),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')
plot(density(x,kernel='cosine',bw=1),main='Cosine',ylim=c(0,0.3),xlim=c(-7,12))
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')

######## Quadratic loss
# number of simulations
J=100
hs=(1:20)/20
s=seq(-10,10,0.01)

h0=5
s0=1001


QUAD_LOSS=function(s,hs,J,n,alpha,l0,l1,p0,p1,h0,s0)
{
ls= length(s)
lh=length(hs)  
EST=array(NA,c(J,ls,lh))
for (j in 1:J)
{
  x=rmixing(n,alpha,l0,l1,p0,p1)
  for (h in 1:lh) 
     {
      EST[j,,h]=sm.density(x,h=hs[h],display='none',ylim=c(0,0.4),nbins=0,eval.points=s)$estimate
     }
}
BIAS=apply(EST,c(2,3),mean)-dmixing(s,alpha,l0,l1,p0,p1)
VAR=apply(EST,c(2,3),var)
EQ=BIAS^2+VAR

nl=2
if (!is.null(s0)) nl=nl+1
  
layout(matrix(c(1:3,rep(4,3),5:(3*nl+1)),byrow=TRUE, ncol=3))
plot(hs,abs(apply(BIAS,2,mean)),type='l',ylab='|BIAS|')
plot(hs,apply(VAR,2,mean),type='l',ylab='VAR')
plot(hs,apply(EQ,2,mean),type='l',ylab='EQ')
abline(v=sm.density(x,method='normal',display="none")$h,col='blue')
abline(v=sm.density(x,method='sj',display="none")$h,col='green')
abline(v=sm.density(x,method='cv',display="none")$h,col='red')


hopt=which(apply(EQ,2,mean)==min(apply(EQ,2,mean)))
if (is.null(h0)) h0=hopt
EST2=EST[,,h0]

plot(s,EST2[1,],type='l',ylab='Estimates',main=paste('h=',hs[h0],sep=''))
for (j in 1:J) lines(s,EST2[j,])
lines(s,dmixing(s,alpha,l0,l1,p0,p1),col='red')

if (!is.null(h0)) 
  {
  plot(s,abs(BIAS[,h0]),type='l',ylab='|BIAS|',main=paste('h=',hs[h0],sep=''))
  plot(s,VAR[,h0],type='l',ylab='VAR',main=paste('h=',hs[h0],sep=''))
  plot(s,EQ[,h0],type='l',ylab='EQ',main=paste('h=',hs[h0],sep=''))
  }
if (!is.null(s0))
  {
  plot(hs,abs(BIAS[s0,]),type='l',ylab='|BIAS|',main=paste('s=',s[s0],sep=''))
  plot(hs,VAR[s0,],type='l',ylab='VAR',main=paste('s=',s[s0],sep=''))
  plot(hs,EQ[s0,],type='l',ylab='EQ',main=paste('s=',s[s0],sep=''))
  }
return(list(BIAS=BIAS,VAR=VAR,EQ=EQ,hopt=hs[hopt]))
}

RES=QUAD_LOSS(s,hs,J,n,alpha,l0,l1,p0,p1,NULL,NULL)

Automatic choice of h

plot(s,dmixing(s,alpha,l0,l1,p0,p1),col='red',type='l')
lines(density(x,kernel='e'),col='blue')
lines(density(x,bw='nrd',kernel='e'))
lines(density(x,bw='SJ',kernel='e'),col='green')
lines(density(x,bw='ucv',kernel='e'),col='orange')

plot(s,dmixing(s,alpha,l0,l1,p0,p1),col='red',type='l')
sm.density(x,method='normal',kernel='e',add=T)
sm.density(x,method='sj',kernel='e',col='green',add=T)
sm.density(x,method='cv',kernel='e',col='orange',add=T)

1.2.2. Come back to real world data

# Galaxies 
hist(galaxies,freq=F,ylim=range(density(galaxies,kernel='e',bw='ucv')$y))
lines(density(galaxies,kernel='rectangular',bw='ucv'),col='blue')
lines(density(galaxies,kernel='e',bw='ucv'),col='orange')

# Faithful 
hist(waiting,freq=F,ylim=range(density(waiting,kernel='e',bw='ucv')$y))
lines(density(waiting,kernel='rectangular',bw='ucv'),col='blue')
lines(density(waiting,kernel='e',bw='ucv'),col='orange')

hist(eruptions,freq=F,ylim=range(density(eruptions,kernel='e',bw='ucv')$y))
lines(density(eruptions,kernel='rectangular',bw='ucv'),col='blue')
lines(density(eruptions,kernel='e',bw='ucv'),col='orange')

## 1.4. Applications
### 1.4.1. Mode estimation

density.mode=function(x,a,b,M,bw='ucv',kernel='e',plot=T)
{
  disc=seq(a,b,length.out=M)
  dens=density(x,from=a,to=b,n=M,bw=bw,kernel=kernel)$y
  mod=disc[(order(dens))[M]]
  max=max(dens)
  if (plot) {plot(disc,dens,type='l')}
  return(list(mode=mod,max=max))
}

sm.mode=function(x,a,b,M,method='cv',plot=T)
{
  disc=seq(a,b,length.out=M)
  display="line"
  if (plot) {display="none"}
  dens=sm.density(x,eval.points=disc,method=method,nbins=0)$estimate
  mod=disc[(order(dens))[M]]
  max=max(dens)
  return(list(mode=mod,max=max))
}

plot(s,dmixing(s,alpha,l0,l1,p0,p1),col='red',type='l')
lines(density(x,bw='ucv',kernel='e'),col='orange')
re=density.mode(x,-10,10,1000,bw='ucv',kernel='e',F)
segments(re$mod,0,re$mod,re$max)
re1=density.mode(x,-10,3,1000,bw='ucv',kernel='e',F)
segments(re1$mod,0,re1$mod,re1$max)

1.4.1. Clustering from density level set

2. Regression estimation

1.1 Experiments

data("pressure")
t=pressure$temperature
p=pressure$pressure
plot(t,p)

res=sm.regression(t,p,method='cv')

names(res)

 [1] "eval.points" "estimate"    "model.y"     "se"          "sigma"      
 [6] "h"           "hweights"    "weights"     "data"        "call"       

res$h

[1] 7.034694

res$eval.points

 [1]   0.000000   7.346939  14.693878  22.040816  29.387755  36.734694
 [7]  44.081633  51.428571  58.775510  66.122449  73.469388  80.816327
[13]  88.163265  95.510204 102.857143 110.204082 117.551020 124.897959
[19] 132.244898 139.591837 146.938776 154.285714 161.632653 168.979592
[25] 176.326531 183.673469 191.020408 198.367347 205.714286 213.061224
[31] 220.408163 227.755102 235.102041 242.448980 249.795918 257.142857
[37] 264.489796 271.836735 279.183673 286.530612 293.877551 301.224490
[43] 308.571429 315.918367 323.265306 330.612245 337.959184 345.306122
[49] 352.653061 360.000000

res$estimate

 [1] 1.999996e-04 5.675549e-04 9.324045e-04 1.671898e-03 3.456180e-03
 [6] 5.149423e-03 1.085884e-02 1.972357e-02 2.846339e-02 4.836680e-02
[11] 7.041253e-02 9.758724e-02 1.635183e-01 2.291547e-01 3.373181e-01
[16] 5.150429e-01 6.883220e-01 1.018800e+00 1.423762e+00 1.837755e+00
[21] 2.665554e+00 3.527883e+00 4.566661e+00 6.266176e+00 7.944415e+00
[26] 1.035057e+01 1.348534e+01 1.658134e+01 2.152680e+01 2.696677e+01
[31] 3.269061e+01 4.175806e+01 5.088840e+01 6.170955e+01 7.610754e+01
[36] 9.033653e+01 1.096621e+02 1.321099e+02 1.545683e+02 1.863899e+02
[41] 2.194453e+02 2.548252e+02 3.022988e+02 3.495649e+02 4.055299e+02
[46] 4.725901e+02 5.391172e+02 6.237562e+02 7.149016e+02 8.060000e+02

n=300
x=runif(n,.1,3)
y=10*exp(-3*x)+7*cos(2*pi*x)/sqrt(x)+rnorm(n,0,2*exp(x/7))
plot(x,y)

res=sm.regression(x,y,h=.007)

res=sm.regression(x,y,h=.07)

res=sm.regression(x,y,h=.7)

res=sm.regression(x,y,method='cv')
z=seq(.1,3,len=1000)
lines(z,10*exp(-3*z)+7*cos(2*pi*z)/sqrt(z),col='red')
sm.regression(x,y,add=T,col='green')
sm.regression(x,y,add=T,col='blue',method='aicc')

names(res)

 [1] "eval.points" "estimate"    "model.y"     "se"          "sigma"      
 [6] "h"           "hweights"    "weights"     "data"        "call"       

res$h

[1] 0.0521145

res$eval.points

 [1] 0.1068731 0.1654310 0.2239890 0.2825469 0.3411048 0.3996627 0.4582207
 [8] 0.5167786 0.5753365 0.6338944 0.6924524 0.7510103 0.8095682 0.8681261
[15] 0.9266841 0.9852420 1.0437999 1.1023578 1.1609158 1.2194737 1.2780316
[22] 1.3365895 1.3951475 1.4537054 1.5122633 1.5708212 1.6293791 1.6879371
[29] 1.7464950 1.8050529 1.8636108 1.9221688 1.9807267 2.0392846 2.0978425
[36] 2.1564005 2.2149584 2.2735163 2.3320742 2.3906322 2.4491901 2.5077480
[43] 2.5663059 2.6248639 2.6834218 2.7419797 2.8005376 2.8590956 2.9176535
[50] 2.9762114

res$estimate

 [1] 24.68945985 14.82203744  7.72429763  2.63854621 -1.67429354 -5.11908185
 [7] -7.25383176 -7.53048540 -6.49775434 -4.90419285 -2.29352556  0.47548620
[13]  3.14650350  5.68113063  7.69170389  8.10760863  7.04094221  5.97199988
[19]  4.52558437  2.50297704 -0.01664411 -2.10043708 -3.77965774 -5.17589053
[25] -5.69361928 -4.93151709 -3.51109380 -1.80304178 -0.42086308  1.06481448
[31]  2.98008281  3.64459870  4.07106946  4.02615227  2.98825097  1.87340426
[37]  0.56020271 -0.45896738 -1.44966844 -2.81222319 -4.02016684 -4.37714947
[43] -3.01109025 -1.63778007 -1.77138612 -1.43005603  0.61823807  2.67657678
[49]  3.57667221  4.79411211

res=sm.regression(x,y,method='cv',eval.points=sort(x))

plot(res$estimate,y[order(x)])
abline(0,1)

1.2 Compute CV

Compute_CV=function(x,y,method,weights=NA)
{
  hopt=sm.regression(x,y,method=method,display='none')$h
  n=length(x)
  if (is.na(weights)) weights=rep(1,n)
  CV=0
  for (i in 1:n)
  {
    CV=CV+(y[i]-sm.regression(x[-i],y[-i],h=hopt,eval.points=x[i],display='none')$estimate)^2*weights[i] 
  }
  return(CV)
}
Compute_CV(x,y,'cv')

[1] 2001.398

1.3 Discrimination

load('Dopage.RData')
hema

 [1] 35.81801 57.69442 54.68696 45.73142 57.02237 57.03312 46.94246 51.21867
 [9] 37.23734 44.22872 43.19019 55.46140 48.31791 43.52074 44.31600 46.83335
[17] 43.98864 45.58534 40.67590 41.64725 46.58286 45.40957 46.14132 41.32432
[25] 51.92825 57.57791 45.10522 49.81766 54.59558 59.94339 46.68741 45.03787
[33] 43.78030 44.32881 47.85452 55.01740 49.64086 44.75317 53.80269 41.54966
[41] 52.33685 43.40389 57.02227 45.73978 56.13190 56.19256 45.33341 56.38840
[49] 45.36690 55.29222 42.93230 44.22211 48.11121 44.84885 46.35119 49.09645
[57] 46.76331 50.31146 44.09170 59.08955 42.67297 47.43126 54.06605 56.46029
[65] 39.68902 52.62125 45.54387 45.69921 53.13572 43.71261 34.72460 48.33640
[73] 51.22310 39.82190 45.93546

test

 [1] "negatif" "positif" "positif" "negatif" "positif" "positif" "negatif"
 [8] "negatif" "negatif" "negatif" "negatif" "positif" "positif" "negatif"
[15] "negatif" "positif" "negatif" "negatif" "negatif" "negatif" "negatif"
[22] "negatif" "negatif" "negatif" "negatif" "positif" "negatif" "positif"
[29] "positif" "positif" "negatif" "negatif" "negatif" "negatif" "negatif"
[36] "positif" "negatif" "negatif" "positif" "negatif" "negatif" "negatif"
[43] "positif" "negatif" "positif" "positif" "negatif" "positif" "negatif"
[50] "positif" "negatif" "negatif" "positif" "negatif" "negatif" "negatif"
[57] "negatif" "negatif" "negatif" "positif" "negatif" "negatif" "positif"
[64] "positif" "negatif" "positif" "negatif" "negatif" "positif" "negatif"
[71] "negatif" "negatif" "positif" "negatif" "negatif"

ROC Curve

ROC=function(z,p1)
{
  p=sort(unique(p1))
  lp=length(p)
  ROC=matrix(NA,lp,2)
  for (s in 1:lp)
  {
    ROC[s,1]=mean(p1[z==0]>p[s])
    ROC[s,2]=mean(p1[z==1]>p[s])  
  }
  AUC=sum(ROC[-lp,2]*(ROC[-lp,1]-ROC[-1,1]))
  
  colnames(ROC)=c("FPR","TPR")
  plot(ROC,type='l',main=paste('AUC =', AUC))
  return(list(ROC=ROC,AUC=AUC))
  }

Estimated Bayes estimator

Classif_NP=function(X,Y,X0=NA)
{
  X=as.matrix(X)
  n=length(Y)
  V=sort(unique(Y))
  n_V=length(V)
  Prob=matrix(NA,n,n_V)
  colnames(Prob)=V
  Class=rep(NA,n)
  if (!is.na(max(X0)))
  {
    X0=as.matrix(X0)
    P0=matrix(NA,nrow(X0),n_V)
    Class0=rep(NA,n)
  }
  for (v in 1:n_V)
  {
    z=as.numeric(Y==V[v])
    for (i in 1:n )
    {
      Prob[i,v]=sm.regression(X[-i,],z[-i],nbins=0, eval.points=X,eval.grid=FALSE,display='none')$estimate[i]
    }
    if (!is.na(max(X0))) {P0[,v]=sm.regression(X,z,eval.points=X0,eval.grid=FALSE,display='none')$estimate}
  }
  if (n_V==2) {ROC(Y==V[2],Prob[,2])}
  Class=V[apply(Prob,1,which.max)]
  
  if (!is.na(max(X0))) {Class0=V[apply(P0,1,which.max)]}
  if (!is.na(max(X0))) {return(list(Class=Class, Prob=Prob, M_table=table(Y,Class), Class0=Class0,Prob0=P0))}
  else {return(list(Class=Class, Prob=Prob, M_table=table(Y,Class)))}
}

(R=Classif_NP(hema,test))

$Class
 [1] "negatif" "positif" "positif" "negatif" "positif" "positif" "negatif"
 [8] "positif" "negatif" "negatif" "negatif" "positif" "negatif" "negatif"
[15] "negatif" "negatif" "negatif" "negatif" "negatif" "negatif" "negatif"
[22] "negatif" "negatif" "negatif" "positif" "positif" "negatif" "negatif"
[29] "positif" "positif" "negatif" "negatif" "negatif" "negatif" "negatif"
[36] "positif" "negatif" "negatif" "positif" "negatif" "positif" "negatif"
[43] "positif" "negatif" "positif" "positif" "negatif" "positif" "negatif"
[50] "positif" "negatif" "negatif" "negatif" "negatif" "negatif" "negatif"
[57] "negatif" "negatif" "negatif" "positif" "negatif" "negatif" "positif"
[64] "positif" "negatif" "positif" "negatif" "negatif" "positif" "negatif"
[71] "negatif" "negatif" "negatif" "negatif" "negatif"

$Prob
           negatif       positif
 [1,]  1.000016213 -1.621285e-05
 [2,] -0.015399975  1.015400e+00
 [3,]  0.131798629  8.682014e-01
 [4,]  0.913610113  8.638989e-02
 [5,] -0.007226770  1.007227e+00
 [6,] -0.007434607  1.007435e+00
 [7,]  0.822317428  1.776826e-01
 [8,]  0.452064700  5.479353e-01
 [9,]  1.000037060 -3.706033e-05
[10,]  0.984403703  1.559630e-02
[11,]  1.002079889 -2.079889e-03
[12,]  0.063393131  9.366069e-01
[13,]  0.766288803  2.337112e-01
[14,]  0.998567461  1.432539e-03
[15,]  0.981894532  1.810547e-02
[16,]  0.870295392  1.297046e-01
[17,]  0.990382634  9.617366e-03
[18,]  0.923010966  7.698903e-02
[19,]  1.002976953 -2.976953e-03
[20,]  1.004961146 -4.961146e-03
[21,]  0.851112730  1.488873e-01
[22,]  0.933663838  6.633616e-02
[23,]  0.884913425  1.150866e-01
[24,]  1.004358978 -4.358978e-03
[25,]  0.376311492  6.236885e-01
[26,] -0.014582062  1.014582e+00
[27,]  0.950271761  4.972824e-02
[28,]  0.653951737  3.460483e-01
[29,]  0.141132940  8.588671e-01
[30,] -0.009277185  1.009277e+00
[31,]  0.842818837  1.571812e-01
[32,]  0.953618386  4.638161e-02
[33,]  0.994551141  5.448859e-03
[34,]  0.981510520  1.848948e-02
[35,]  0.748594649  2.514054e-01
[36,]  0.100112801  8.998872e-01
[37,]  0.603680357  3.963196e-01
[38,]  0.966419808  3.358019e-02
[39,]  0.229610874  7.703891e-01
[40,]  1.004793373 -4.793373e-03
[41,]  0.331009187  6.689908e-01
[42,]  0.999996640  3.359742e-06
[43,] -0.007224877  1.007225e+00
[44,]  0.913058111  8.694189e-02
[45,]  0.022202558  9.777974e-01
[46,]  0.019332617  9.806674e-01
[47,]  0.938043481  6.195652e-02
[48,]  0.010994495  9.890055e-01
[49,]  0.936136018  6.386398e-02
[50,]  0.076533578  9.234664e-01
[51,]  1.003810698 -3.810698e-03
[52,]  0.984586259  1.541374e-02
[53,]  0.780585234  2.194148e-01
[54,]  0.962360848  3.763915e-02
[55,]  0.869140063  1.308599e-01
[56,]  0.649572310  3.504277e-01
[57,]  0.836751458  1.632485e-01
[58,]  0.542609205  4.573908e-01
[59,]  0.987977417  1.202258e-02
[60,] -0.015773466  1.015773e+00
[61,]  1.004855560 -4.855560e-03
[62,]  0.782654530  2.173455e-01
[63,]  0.199051930  8.009481e-01
[64,]  0.008278253  9.917217e-01
[65,]  1.001237547 -1.237547e-03
[66,]  0.369378542  6.306215e-01
[67,]  0.925590318  7.440968e-02
[68,]  0.915724000  8.427600e-02
[69,]  0.308905764  6.910942e-01
[70,]  0.995717811  4.282189e-03
[71,]  1.000003364 -3.363682e-06
[72,]  0.710345424  2.896546e-01
[73,]  0.522081041  4.779190e-01
[74,]  1.001423773 -1.423773e-03
[75,]  0.899719315  1.002807e-01

$M_table
         Class
Y         negatif positif
  negatif      47       3
  positif       5      20

plot(hema,R$Prob[,2])

Classif_NP(hema,test,c(37,42,58,57))

$Class
 [1] "negatif" "positif" "positif" "negatif" "positif" "positif" "negatif"
 [8] "positif" "negatif" "negatif" "negatif" "positif" "negatif" "negatif"
[15] "negatif" "negatif" "negatif" "negatif" "negatif" "negatif" "negatif"
[22] "negatif" "negatif" "negatif" "positif" "positif" "negatif" "negatif"
[29] "positif" "positif" "negatif" "negatif" "negatif" "negatif" "negatif"
[36] "positif" "negatif" "negatif" "positif" "negatif" "positif" "negatif"
[43] "positif" "negatif" "positif" "positif" "negatif" "positif" "negatif"
[50] "positif" "negatif" "negatif" "negatif" "negatif" "negatif" "negatif"
[57] "negatif" "negatif" "negatif" "positif" "negatif" "negatif" "positif"
[64] "positif" "negatif" "positif" "negatif" "negatif" "positif" "negatif"
[71] "negatif" "negatif" "negatif" "negatif" "negatif"

$Prob
           negatif       positif
 [1,]  1.000016213 -1.621285e-05
 [2,] -0.015399975  1.015400e+00
 [3,]  0.131798629  8.682014e-01
 [4,]  0.913610113  8.638989e-02
 [5,] -0.007226770  1.007227e+00
 [6,] -0.007434607  1.007435e+00
 [7,]  0.822317428  1.776826e-01
 [8,]  0.452064700  5.479353e-01
 [9,]  1.000037060 -3.706033e-05
[10,]  0.984403703  1.559630e-02
[11,]  1.002079889 -2.079889e-03
[12,]  0.063393131  9.366069e-01
[13,]  0.766288803  2.337112e-01
[14,]  0.998567461  1.432539e-03
[15,]  0.981894532  1.810547e-02
[16,]  0.870295392  1.297046e-01
[17,]  0.990382634  9.617366e-03
[18,]  0.923010966  7.698903e-02
[19,]  1.002976953 -2.976953e-03
[20,]  1.004961146 -4.961146e-03
[21,]  0.851112730  1.488873e-01
[22,]  0.933663838  6.633616e-02
[23,]  0.884913425  1.150866e-01
[24,]  1.004358978 -4.358978e-03
[25,]  0.376311492  6.236885e-01
[26,] -0.014582062  1.014582e+00
[27,]  0.950271761  4.972824e-02
[28,]  0.653951737  3.460483e-01
[29,]  0.141132940  8.588671e-01
[30,] -0.009277185  1.009277e+00
[31,]  0.842818837  1.571812e-01
[32,]  0.953618386  4.638161e-02
[33,]  0.994551141  5.448859e-03
[34,]  0.981510520  1.848948e-02
[35,]  0.748594649  2.514054e-01
[36,]  0.100112801  8.998872e-01
[37,]  0.603680357  3.963196e-01
[38,]  0.966419808  3.358019e-02
[39,]  0.229610874  7.703891e-01
[40,]  1.004793373 -4.793373e-03
[41,]  0.331009187  6.689908e-01
[42,]  0.999996640  3.359742e-06
[43,] -0.007224877  1.007225e+00
[44,]  0.913058111  8.694189e-02
[45,]  0.022202558  9.777974e-01
[46,]  0.019332617  9.806674e-01
[47,]  0.938043481  6.195652e-02
[48,]  0.010994495  9.890055e-01
[49,]  0.936136018  6.386398e-02
[50,]  0.076533578  9.234664e-01
[51,]  1.003810698 -3.810698e-03
[52,]  0.984586259  1.541374e-02
[53,]  0.780585234  2.194148e-01
[54,]  0.962360848  3.763915e-02
[55,]  0.869140063  1.308599e-01
[56,]  0.649572310  3.504277e-01
[57,]  0.836751458  1.632485e-01
[58,]  0.542609205  4.573908e-01
[59,]  0.987977417  1.202258e-02
[60,] -0.015773466  1.015773e+00
[61,]  1.004855560 -4.855560e-03
[62,]  0.782654530  2.173455e-01
[63,]  0.199051930  8.009481e-01
[64,]  0.008278253  9.917217e-01
[65,]  1.001237547 -1.237547e-03
[66,]  0.369378542  6.306215e-01
[67,]  0.925590318  7.440968e-02
[68,]  0.915724000  8.427600e-02
[69,]  0.308905764  6.910942e-01
[70,]  0.995717811  4.282189e-03
[71,]  1.000003364 -3.363682e-06
[72,]  0.710345424  2.896546e-01
[73,]  0.522081041  4.779190e-01
[74,]  1.001423773 -1.423773e-03
[75,]  0.899719315  1.002807e-01

$M_table
         Class
Y         negatif positif
  negatif      47       3
  positif       5      20

$Class0
[1] "negatif" "negatif" "positif" "positif"

$Prob0
             [,1]          [,2]
[1,]  1.000022248 -2.224845e-05
[2,]  1.004992295 -4.992295e-03
[3,] -0.014065877  1.014066e+00
[4,] -0.006230051  1.006230e+00

3. Multivariate and Functional case

fdaplot=function(X,ylim=NULL)
{
  matplot(t(X),type='l',ylim=ylim)
}

library("fda")

Le chargement a nécessité le package : splines

Le chargement a nécessité le package : fds

Le chargement a nécessité le package : rainbow

Le chargement a nécessité le package : pcaPP

Le chargement a nécessité le package : RCurl

Le chargement a nécessité le package : deSolve

Attachement du package : 'fda'

L'objet suivant est masqué depuis 'package:graphics':

    matplot

TECATOR=read.table('npfda-spectrometric.dat')
CURVES=as.matrix(TECATOR[,-101])
Y=TECATOR[,101]
fdaplot(CURVES)

source('npfda.R')
R0=funopare.kernel.cv(Y,CURVES, CURVES,0,10,c(0,1))
R0$Mse

[1] 117.2638

plot(R0$Estimated.values,Y)
abline(0,1)

R1=funopare.kernel.cv(Y,CURVES, CURVES,1,10,c(0,1))
R1$Mse

[1] 44.27052

plot(R1$Estimated.values,Y)
abline(0,1)

R2=funopare.kernel.cv(Y,CURVES, CURVES,2,10,c(0,1))
R2$Mse

[1] 5.784178

plot(R2$Estimated.values,Y)
abline(0,1)

R3=funopare.kernel.cv(Y,CURVES, CURVES,3,10,c(0,1))
R3$Mse

[1] 9.118275

plot(R3$Estimated.values,Y)
abline(0,1)

Rpca=funopare.kernel.cv(Y,CURVES, CURVES,semimetric='pca',10)
Rpca$Mse

[1] 122.7028

plot(Rpca$Estimated.values,Y)
abline(0,1)

load("PHONEMES.RData")
par(mfrow=c(2,3))
Phon=sort(unique(PHONEME))
for (j in Phon)
{
fdaplot(CURVES[PHONEME==j,],ylim=range(CURVES))
title(j)
}