TP_Chapitre_3.1
2025-01-28
Multivariate case
Fonctions
library(sm)Package 'sm', version 2.2-6.0: type help(sm) for summary informationCompute_CV=function(x,y,method,weights=NA)
{
hopt=sm.regression(x,y,method=method,display='none')$h
n=nrow(x)
yest=rep(NA,n)
if (is.na(weights)) weights=rep(1,n)
CV=0
for (i in 1:n)
{
yest[i]=(sm.regression(x[-i,],y[-i],h=hopt,eval.points=x,eval.grid=F,display='none',nbins=0)$estimate)[i]
CV=CV+(y[i]-yest[i])^2*weights[i]
}
return(list('CV'=CV,'y_est'=yest))
}
####
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))
}
#######
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], 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)]
M_table=table(Y,Class)
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=M_table, err=1-sum(diag(M_table))/n, Class0=Class0,Prob0=P0))}
else {return(list(Class=Class, Prob=Prob, M_table=M_table, err=1-sum(diag(M_table))/n))}
}Density
Simulation
library(sm)
#library(plotly)
library(plot3D)
#library(rpanel)
n=1000
alpha=0.4
mu11=7
mu12=0
sig11=1
sig12=3
mu21=0
mu22=3
sig21=2
sig22=2
x1=matrix(c(rnorm(n,mu11,sig11),rnorm(n,mu12,sig12)),ncol=2)
x2=matrix(c(rnorm(n,mu21,sig21),rnorm(n,mu22,sig22)),ncol=2)
z=rbinom(n,1,alpha)
x=x1
x[z==1,]=x2[z==1,]
M=300
r1=seq(-10,10,len=M)
r2=seq(-10,10,len=M)
x_1=x[,1]
x_2=x[,2]Estimation de la densité
par(mfrow=c(1,3))
sm.density(x,method='cv',structure='common',nbins=0)
sm.density(x,method='cv',structure="scaled",nbins=0)
sm.density(x,method='cv',structure='separate',nbins=0)
par(mfrow=c(1,2))
rx=sm.density(x,method='cv',eval.points=x
,eval.grid=F,
nbins=0)
r=sm.density(x,method='cv',eval.points=cbind(r1,r2)
,eval.grid=T,
nbins=0)
f_k=r$estimate
#fra=data.frame(x_1,x_2,f_k)
Comparaison avec la vraie densité
f=matrix(NA,M,M)
for (j in 1:M)
{
for (k in 1:M)
{
f[j,k]=(1-alpha)*dnorm(r1[j],mu11,sig11)*dnorm(r2[k],mu12,sig12)+(1-alpha)*dnorm(r1[j],mu21,sig21)*dnorm(r2[k],mu22,sig22)
}
}
v1=rep(r1,M)
v2=rep(r2,rep(M,M))
fv=as.vector(f)
fra=data.frame(v1,v2,fv)
par(mfrow=c(1,3))
image(f,zlim=range(c(f,r$estimate)))
image(r$estimate,zlim=range(c(f,r$estimate)))
image(abs(f-r$estimate),zlim=range(c(f,r$estimate)))
surf3D(matrix(rep(r1,M),M,M),matrix(rep(r2,rep(M,M)),M,M),f,clim=range(c(f,f_k)))
surf3D(matrix(rep(r1,M),M,M),matrix(rep(r2,rep(M,M)),M,M),f_k,clim=range(c(f,f_k)))
surf3D(matrix(rep(r1,M),M,M),matrix(rep(r2,rep(M,M)),M,M),abs(f-f_k),clim=range(c(f,f_k)))
Estimation du mode
p_mo=which(f==max(f),arr.ind=T)
(mode=c(r1[p_mo[1]],r2[p_mo[2]]))[1] 6.98996656 0.03344482p_mok=which(f_k==max(f_k),arr.ind=T)
(modek=c(r1[p_mok[1]],r2[p_mok[2]]))[1] 7.25752508 0.03344482p_mo1=which(f[1:(M/2),]==max(f[1:(M/2),]),arr.ind=T)
(modek1=c(r1[p_mo1[1]],r2[p_mo1[2]]))[1] -0.03344482 2.97658863p_mok1=which(f_k[1:(M/2),]==max(f_k[1:(M/2),]),arr.ind=T)
(modek1=c(r1[p_mok1[1]],r2[p_mok1[2]]))[1] -0.03344482 2.57525084#plot_ly(fra,x=~v1,y=~v2,z = ~fv, type = "surface")Real world data
Estimation de la densité
data(faithful)
attach(faithful)
par(mfrow=c(1,3))
sm.density(faithful,method='cv',structure='common',nbins=0)
sm.density(faithful,method='cv',structure="scaled",nbins=0)
sm.density(faithful,method='cv',structure='separate',nbins=0)
Estimation du mode
x=eruptions
y=waiting
x0=seq(min(x),max(x),len=100)
y0=seq(min(y),max(y),len=100)
re=sm.density(cbind(x,y),method='cv',eval.points=cbind(x0,y0),
eval.grid=T,structure='common',nbins=0)
p_mo=which(re$estimate==max(re$estimate),arr.ind=T)
(mode=c(x0[p_mo[1]],y0[p_mo[2]]))[1] 4.428283 81.010101p_mo1=which(re$estimate[1:30,1:30]==max(re$estimate[1:30,1:30]),arr.ind=T)
(mode1=c(x0[p_mo1[1]],y0[p_mo1[2]]))[1] 1.918182 53.171717Regression
Real world data
load('Randonnee.RData')
x=long
y=lat
z=alti # variable of interest
x0=seq(min(x),max(x),len=100)
y0=seq(min(y),max(y),len=100)
sm.regression(cbind(x,y),z,method='cv',eval.points=cbind(x0,y0),
eval.grid=T,structure='common',nbins=0)
alti_profile=sm.regression(cbind(x,y),z,method='cv',eval.points=Position,
eval.grid=F,structure='common',nbins=0)$estimate
km=rep(0,nrow(Position))
for (j in 1:(nrow(Position)-1)) {km[j+1]=km[j]+sqrt((Position[j+1,1]-Position[j,1])^2+(Position[j+1,2]-Position[j,2])^2)}
#km=cumsum(km)/1000
plot(km/10,alti_profile,type='l')
Qualité
CV with Kernel Estimator
RE=Compute_CV(cbind(x,y),z,method='cv')
RE$CV/length(alti)[1] 128.3074plot(RE$y_est,z)
abline(0,1)
CV avec moyenne Z
n=length(z)
CV_z=0
for (j in 1:n)
{
CV_z=CV_z+(z[j]-mean(z[-j]))^2
}
CV_z/n[1] 1582.402Classification
data(iris)
v=combn(1:4,2)
for(j in 1:ncol(v))
{
print(colnames(iris)[v[,j]])
C_NP=Classif_NP(iris[,v[,j]],iris[,5])
print(C_NP$err)
print(C_NP$M_table)
}[1] "Sepal.Length" "Sepal.Width"
[1] 0.22
Class
Y setosa versicolor virginica
setosa 49 1 0
versicolor 2 32 16
virginica 0 14 36
[1] "Sepal.Length" "Petal.Length"
[1] 0.05333333
Class
Y setosa versicolor virginica
setosa 50 0 0
versicolor 0 47 3
virginica 0 5 45
[1] "Sepal.Length" "Petal.Width"
[1] 0.04666667
Class
Y setosa versicolor virginica
setosa 50 0 0
versicolor 0 48 2
virginica 0 5 45
[1] "Sepal.Width" "Petal.Length"
[1] 0.05333333
Class
Y setosa versicolor virginica
setosa 50 0 0
versicolor 0 46 4
virginica 0 4 46
[1] "Sepal.Width" "Petal.Width"
[1] 0.06666667
Class
Y setosa versicolor virginica
setosa 49 1 0
versicolor 0 45 5
virginica 0 4 46
[1] "Petal.Length" "Petal.Width"
[1] 0.04
Class
Y setosa versicolor virginica
setosa 50 0 0
versicolor 0 47 3
virginica 0 3 47