Kernel estimation
1. Density estimation
1.1. Histogram
1.1.1. Real world data
1.1.1.1. Galaxies
library('sm')
library(MASS)
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=400
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",nbins=0)$h,col='blue')
abline(v=sm.density(x,method='sj',display="none",nbins=0)$h,col='green')
abline(v=sm.density(x,method='cv',display="none",nbins=0)$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
par(oma=c(0,0,0,7))
plot(s,dmixing(s,alpha,l0,l1,p0,p1),col='red',type='l')
lines(density(x,kernel='e'),col='blue',lty=2)
lines(density(x,bw='nrd',kernel='e'),lty=2)
lines(density(x,bw='SJ',kernel='e'),col='green',lty=2)
lines(density(x,bw='ucv',kernel='e'),col='orange',lty=2)
legend(par('usr')[2], par('usr')[4],xpd=NA,
legend = c('f','nrd0','nrd','SJ','ucv'),
col = c('red','blue','black','green','orange'), lty=c(1,rep(2,4)))
par(oma=c(0,0,0,7))
plot(s,dmixing(s,alpha,l0,l1,p0,p1),col='red',type='l')
sm.density(x,method='normal',add=T,nbins=0)
sm.density(x,method='sj',col='green',add=T,nbins=0)
sm.density(x,method='cv',col='orange',add=T,nbins=0)
legend(par('usr')[2], par('usr')[4],xpd=NA,
legend = c('f','normal','sj','ucv'),
col = c('red','black','green','orange'), lty=c(1,rep(2,3)))
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,display=display)$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)
plot(s,dmixing(s,alpha,l0,l1,p0,p1),col='red',type='l')
sm.density(x,method='cv',col='orange',add=T)
re=sm.mode(x,-10,10,1000,method='cv',plot=F)
segments(re$mod,0,re$mod,re$max)
re1=sm.mode(x,-10,3,1000,method='cv',plot=F)
segments(re1$mod,0,re1$mod,re1$max)
Finding extrema
finding.extrema=function(x,tau=10,a=min(x),b=max(x),M=1000,method='cv',plot=T)
{
disc=seq(a,b,length.out=M)
n=length(x)
dens=sm.density(x,eval.points= disc,method=method,nbins=0,display='none')$estimate
inc=rep(0,M)
valleys=rep(0,M)
hills=rep(0,M)
kv=1
kh=1
for (j in 1:M)
{
i1=max(1,j-tau)
i2=min(M,j+tau)
ratio=(dens[i2]-dens[i1])/(disc[i2]-disc[i1])
if (ratio>0) {inc[j]=1}
if (j>1) {if (inc[j]==0 & inc[j-1]==1 ) {kv=kv+1}}
if (j>1) {if (inc[j]==1 & inc[j-1]==0 ) {kh=kh+1}}
valleys[j]=kv
hills[j]=kh
}
if (plot) {plot(disc,dens,type='l')}
kv=max(valleys)
argmin=rep(NA,kv)
min=rep(NA,kv)
kh=max(hills)
argmax=rep(NA,kh)
max=rep(NA,kh)
for (j in 1:kv)
{
candidates=disc[valleys==j]
argmin[j]=candidates[which.min(dens[valleys==j])]
min[j]=min(dens[valleys==j])
if (plot) {segments(argmin[j],0,argmin[j],min[j],col='darkgreen')}
}
for (j in 1:kh)
{
candidates=disc[hills==j]
argmax[j]=candidates[which.max(dens[hills==j])]
max[j]=max(dens[hills==j])
if (plot) {segments(argmax[j],0,argmax[j],max[j],col='orange')}
}
return(list(argmax=argmax,max=max,argmin=argmin,min=min))
}
finding.extrema(x)
$argmax
[1] -5.1750653 -1.5475745 -0.3159148 2.1305325 3.8177376 7.8164135
$max
[1] 0.004975100 0.064995576 0.075113138 0.032013548 0.006760535 0.271930468
$argmin
[1] -6.035540 -4.297719 -1.024541 1.472523 3.429680 4.138307 10.819638
$min
[1] 0.003358290 0.001010069 0.060101470 0.025647813 0.005913254 0.006390660
[7] 0.007943269finding.extrema(x,100)
$argmax
[1] -5.6137386 -0.3159148 7.8164135
$max
[1] 0.00391354 0.07511314 0.27193047
$argmin
[1] -4.297719 3.429680 10.819638
$min
[1] 0.001010069 0.005913254 0.007943269On Galaxies and Faithful
$argmax
[1] 9649.511 16133.601 19878.288 22592.559 26663.964 32444.354 34153.339
$max
[1] 4.468359e-05 1.569090e-05 1.855198e-04 1.231445e-04 1.690490e-05
[6] 1.333955e-05 8.416955e-06
$argmin
[1] 9172.00 13243.41 17189.15 21461.61 26186.45 29529.03 33675.83 34279.00
$min
[1] 3.660048e-05 5.959553e-10 6.418702e-06 1.045646e-04 1.623745e-05
[6] 3.827367e-09 7.972951e-06 8.304511e-06
$argmax
[1] 1.873273 2.861261 4.483383
$max
[1] 0.62140065 0.04320936 0.61947705
$argmin
[1] 1.600000 2.728128 3.110010 5.100000
$min
[1] 0.10429578 0.03800498 0.02169283 0.06940804
$argmax
[1] 53.23924 80.03103
$max
[1] 0.02149506 0.04052094
$argmin
[1] 43.00000 66.87387 96.00000
$min
[1] 0.004585450 0.008997603 0.0021666371.4.2. Clustering from density level set
sm.clustering.level.sets=function(x,level=0.20,a=min(x),b=max(x),M=1000,method='cv',plot=T)
{
disc=seq(a,b,length.out=M)
n=length(x)
o=order(x)
x.and.disc=c(x,disc)
#type=
names(x.and.disc)=rep(c('data','disc'),c(n,M))
x1=sort(x.and.disc)
type1=names(x1)
n1=length(x1)
dens=sm.density(x,eval.points=x1,method=method,nbins=0,display='none')$estimate
adjusted.level=quantile(dens,level)
is.over.level=dens>adjusted.level
cluster=rep(0,M+n)
k=1
for (j in 1:(M+n))
{
if (j>1) {if (is.over.level[j]==0 & is.over.level[j-1]==1 ) {k=k+1}}
if (is.over.level[j]==1) {cluster[j]=k}
}
if (plot) {par(oma=c(0, 0, 0, 7))
plot(x1,dens,type='l')
abline(adjusted.level,0,col='orange',lty=2)}
k=max(cluster)
cluster.bounds=matrix(NA,2,k)
palette=rainbow(k)
for (j in 1:k)
{
cluster.bounds[,j]=range(x1[cluster==j])
if (plot) { segments(cluster.bounds[1,j],adjusted.level,cluster.bounds[2,j],adjusted.level,col=palette[j])}
}
if (plot) { legend(par('usr')[2], par('usr')[4],xpd=NA,#"topright", inset=c(-0.2,0),
legend = c('No Cluster',paste('Cluster ',1:k)),
col = c('orange',palette), lty=c(2,rep(1,length((unique(cluster))))))}
cluster.on.data=(cluster[type1=='data'])
cluster.on.data[cluster.on.data==0]=NA
return(list(cluster=cluster.on.data[o],cluster.bounds=cluster.bounds))
}
sm.clustering.level.sets(x,0.3)
$cluster
[1] 2 1 2 2 2 1 1 2 2 2 2 1 1 NA NA 2 2 2 2 2 1 2 2 2 2
[26] 2 1 2 1 2 2 NA 2 1 1 2 NA 2 2 1 2 2 2 1 2 2 2 2 2 2
[51] 2 2 2 NA 2 2 2 1 NA 2 2 2 2 2 2 2 2 2 2 2 2 2 1 NA 1
[76] 2 NA 2 2 2 2 1 1 2 2 2 1 2 1 2 2 2 1 1 2 1 1 2 1 2
[101] 2 2 2 2 1 2 2 2 1 2 2 2 2 2 2 2 1 2 1 2 NA 2 1 2 2
[126] 1 2 2 2 1 2 1 2 2 2 1 2 2 2 2 1 2 2 2 2 1 2 2 2 1
[151] 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 1 2 2 2 1 2 2 NA 2 2
[176] NA 2 2 2 2 2 1 2 1 2 1 1 2 2 1 1 2 2 2 2 1 NA 2 1 2
[201] 2 2 2 NA 2 2 2 1 2 2 1 2 1 1 1 2 2 2 2 1 2 2 1 2 2
[226] 2 2 1 2 2 1 2 2 2 2 2 2 1 2 1 2 1 1 2 2 2 2 2 1 2
[251] 2 2 2 2 2 2 2 2 2 1 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2
[276] 2 2 2 2 2 2 2 1 2 1 2 1 1 2 2 1 1 2 1 1 2 2 2 1 2
[301] 2 2 2 1 2 2 2 2 2 2 2 2 1 2 1 2 2 2 NA 2 2 2 2 NA 2
[326] 1 1 1 2 1 2 2 2 2 NA 2 2 1 1 1 2 1 2 1 2 1 2 NA 1 2
[351] 2 2 2 1 NA 2 2 2 1 2 1 1 2 2 2 2 2 1 2 2 2 1 2 1 1
[376] NA 1 2 2 1 1 1 1 2 1 1 2 2 2 NA 2 2 2 1 2 2 2 2 2 2
$cluster.bounds
[,1] [,2]
[1,] -2.711746 5.521815
[2,] 2.653566 10.245989sm.clustering.level.sets(x,0.7)
$cluster
[1] 1 NA 1 1 1 NA NA 1 1 NA 1 NA NA NA NA 1 1 1 1 1 NA 1 1 NA 1
[26] 1 NA 1 NA 1 1 NA 1 NA NA 1 NA 1 1 NA 1 1 1 NA 1 1 1 1 1 1
[51] 1 NA 1 NA 1 1 1 NA NA 1 NA 1 1 1 1 1 1 1 1 1 1 1 NA NA NA
[76] NA NA 1 1 1 1 NA NA NA 1 1 NA NA NA 1 NA 1 NA NA 1 NA NA 1 NA 1
[101] 1 1 1 NA NA 1 1 1 NA 1 NA NA NA 1 1 1 NA 1 NA 1 NA NA NA 1 NA
[126] NA 1 1 1 NA 1 NA 1 1 1 NA 1 1 1 1 NA 1 NA 1 1 NA 1 1 1 NA
[151] 1 1 1 1 NA 1 1 NA NA NA 1 1 1 1 1 NA 1 1 1 NA NA 1 NA 1 NA
[176] NA 1 NA 1 1 1 NA 1 NA 1 NA NA 1 1 NA NA 1 NA 1 1 NA NA 1 NA 1
[201] 1 1 NA NA 1 1 1 NA 1 1 NA 1 NA NA NA 1 1 NA 1 NA 1 1 NA 1 1
[226] 1 1 NA 1 1 NA NA 1 1 1 1 NA NA 1 NA 1 NA NA NA NA 1 1 1 NA 1
[251] 1 NA 1 1 1 1 1 1 1 NA NA 1 1 1 NA 1 NA 1 1 1 1 1 1 1 1
[276] 1 NA NA 1 1 1 1 NA 1 NA 1 NA NA 1 1 NA NA 1 NA NA 1 1 1 NA 1
[301] 1 1 1 NA 1 1 1 1 1 1 1 1 NA 1 NA 1 1 1 NA NA 1 1 1 NA 1
[326] NA NA NA 1 NA 1 1 1 NA NA 1 1 NA NA NA 1 NA 1 NA 1 NA 1 NA NA 1
[351] 1 1 1 NA NA 1 1 1 NA 1 NA NA 1 1 1 NA 1 NA 1 1 1 NA 1 NA NA
[376] NA NA 1 1 NA NA NA NA 1 NA NA 1 1 1 NA 1 1 1 NA 1 1 1 1 1 1
$cluster.bounds
[,1]
[1,] 6.567882
[2,] 9.554235sm.clustering.level.sets(x,0.1)
$cluster
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 NA 1 1 1 1 1 1 1 1 1 1 1
[26] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[51] 1 1 1 NA 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[76] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[101] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[126] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[151] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[176] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[201] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[226] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[251] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[276] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[301] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[326] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[351] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[376] NA 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
$cluster.bounds
[,1]
[1,] -3.724069
[2,] 10.819638On Galaxies and Faithful
$cluster
[1] 1 1 1 1 1 1 1 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[39] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[77] 3 3 3 4 4 4
$cluster.bounds
[,1] [,2] [,3] [,4]
[1,] 9172.0 15329.37 17289.68 31564.73
[2,] 11207.7 17088.62 27644.12 34279.00
$cluster
[1] 1 1 2 1 1 1 1 1 1 2 2 2 1 1 2 2 1 2 NA 1 1 1 NA 2 2
[26] 2 2 2 1 NA NA NA 2 2 2 2 1 2 2 2 2 2 2 2 1 1 1 2 1 2
[51] 2 1 2 2 2 1 1 2 2 2 1 2 2 1 NA 2 2 NA 1 1 2 1 2 2 2
[76] 1 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 NA 2 1 2 1 2 1 1 1
[101] 2 1 2 1 2 2 2 2 NA 1 1 1 2 1 2 NA 2 2 2 1 1 2 2 2 1
[126] 2 1 1 2 2 1 2 2 2 1 1 2 2 2 2 1 1 NA 1 2 1 2 2 2 2
[151] 1 2 2 2 1 2 2 2 2 2 1 2 2 2 2 1 2 2 2 1 2 2 2 2 2
[176] 2 2 1 NA 1 1 NA 2 2 2 1 1 2 2 2 NA 2 2 2 1 2 2 2 NA 1
[201] 2 2 2 2 1 NA 2 1 1 NA 2 2 2 2 2 1 1 1 1 2 2 2 1 1 1
[226] 2 2 2 2 2 2 2 2 2 1 2 2 2 NA 2 1 1 1 1 2 2 1 2 2 2
[251] 1 1 1 1 2 2 2 NA 2 1 1 2 1 NA 2 1 2 2 2 2 1 2
$cluster.bounds
[,1] [,2]
[1,] 1.631532 3.567000
[2,] 2.437337 4.991391
$cluster
[1] NA 2 2 2 2 2 2 NA NA 1 1 NA 2 1 1 NA NA 2 2 2 2 NA 2 2 2
[26] 2 NA 1 NA 2 2 2 1 1 1 2 2 1 2 2 2 2 2 2 NA 1 1 1 2 2
[51] 2 2 2 NA 1 2 2 2 2 1 2 2 2 2 2 2 1 1 2 NA 1 1 1 2 2
[76] 2 2 1 1 NA 2 2 NA NA 1 1 NA NA 2 2 1 1 2 2 1 1 2 1 2 2
[101] 2 1 2 1 2 2 2 1 1 2 2 2 NA 1 1 1 1 2 2 2 NA 1 1 2 2
[126] NA 1 1 NA 2 2 2 2 NA 1 1 1 1 2 2 2 2 2 1 1 2 2 2 2 2
[151] 2 2 2 2 2 1 1 1 1 1 1 NA 2 2 2 2 2 2 2 NA NA 1 1 1 2
[176] 2 2 2 2 2 1 1 1 1 NA 2 2 2 1 NA 2 2 2 2 2 2 2 2 2 2
[201] NA 1 1 1 1 NA 2 2 2 2 2 2 2 1 1 1 1 1 NA 2 2 2 2 2 2
[226] 2 NA 1 1 1 1 NA NA 2 2 2 NA NA NA 1 2 2 2 NA 2 2 2 2 2 2
[251] 2 NA 1 NA 2 2 2 2 2 2 1 1 NA 2 2 NA 2 1 2 2 2 2
$cluster.bounds
[,1] [,2]
[1,] 46.71371 70.16316
[2,] 62.04605 89.686691.4.3. Clustering from density extrema
sm.clustering.extrema=function(x,tau=10,a=min(x),b=max(x),M=1000,method='sj',plot=T)
{
disc=seq(a,b,length.out=M)
n=length(x)
dens=sm.density(x,eval.points= disc,method=method,nbins=0,display='none')$estimate
inc=rep(0,M)
valleys=rep(0,M)
k=1
for (j in 1:M)
{
i1=max(1,j-tau)
i2=min(M,j+tau)
ratio=(dens[i2]-dens[i1])/(disc[i2]-disc[i1])
if (ratio>0) {inc[j]=1}
if (j>1) {if (inc[j]==0 & inc[j-1]==1 ) {k=k+1}}
valleys[j]=k
}
if (plot) {par(oma=c(0, 0, 0, 4))
plot(disc,dens,type='l')}
k=max(valleys)
thresholds=rep(NA,k)
mins=rep(NA,k)
for (j in 1:k)
{
candidates=disc[valleys==j]
thresholds[j]=candidates[which.min(dens[valleys==j])]
mins[j]=min(dens[valleys==j])
#if (plot) {segments(thresholds[j],0,thresholds[j],mins[j],col='orange')}
}
n=length(x)
cluster=rep(0,n)
for (j in 1:k) {if (j==1) {cluster[x==thresholds[j]]=j}
cluster[x>thresholds[j]]=j}
cluster.disc=rep(0,M)
for (j in 1:k) {if (j==1) {cluster.disc[disc==thresholds[j]]=j}
cluster.disc[disc>thresholds[j]]=j}
K=max(cluster.disc)
palette=rainbow(K+1)
if (plot) {for (j in 0:K) {if (sum(cluster.disc==j)!=0) {
disc.j=disc[cluster.disc==j]
dens.j=dens[cluster.disc==j]
n.j=sum(cluster.disc==j)
polygon(c(disc.j[1],disc.j,disc.j[n.j]),c(0,dens.j,0),
col = palette[j+1])}
}
#par(xpd=TRUE)
legend(par('usr')[2], par('usr')[4],xpd=NA,#"topright", inset=c(-0.2,0),
legend = paste('Cluster ',sort(unique(cluster.disc))),
fill = palette[sort(unique(cluster.disc))+1])
}
#cluster.on.data[cluster.on.data==0]=NA
return(list(cluster=cluster,thresholds=thresholds))
}
sm.clustering.extrema(x,10)
$cluster
[1] 2 5 3 5 2 2 5 5 5 5 5 5 5 5 5 5 5 5 3 5 5 5 3 2 5 5 5 5 5 2 5 5 4 5 2 5 5
[38] 5 5 5 3 5 3 2 5 5 5 3 3 5 5 4 5 5 5 5 3 5 5 5 3 2 5 5 5 5 2 5 5 5 2 5 5 5
[75] 5 3 5 3 5 5 5 5 3 5 5 2 5 5 5 5 5 4 5 1 5 2 5 5 5 5 5 5 5 2 5 5 5 5 5 3 5
[112] 5 5 5 5 5 5 3 5 5 5 2 5 5 5 5 4 3 5 4 2 5 5 5 5 3 2 5 5 5 4 5 3 5 5 2 3 5
[149] 3 5 5 5 5 5 5 5 5 5 5 5 3 5 5 5 5 5 5 5 5 5 5 3 5 5 2 5 5 5 5 5 3 5 2 3 2
[186] 3 5 4 5 5 5 5 5 5 5 5 5 5 2 4 4 5 5 5 4 5 5 5 5 5 5 2 3 5 5 5 5 3 5 5 5 5
[223] 2 5 5 5 3 5 3 5 2 5 5 5 3 3 5 3 5 5 3 5 3 5 5 5 2 5 5 5 5 1 5 5 5 5 3 5 5
[260] 5 3 5 5 3 5 2 5 5 5 5 3 5 5 5 5 5 3 5 5 5 4 5 3 5 5 5 5 5 5 5 5 5 5 5 5 5
[297] 4 3 5 5 3 5 2 5 5 5 2 5 5 5 4 5 5 4 2 3 5 5 5 3 5 5 5 5 2 4 5 5 3 5 4 5 5
[334] 5 5 5 5 5 3 5 5 5 5 5 2 4 3 5 2 2 5 5 5 2 5 2 3 5 3 2 5 5 5 5 5 5 5 3 3 5
[371] 5 5 5 5 2 3 5 5 5 1 5 2 5 2 5 3 5 5 3 5 3 5 5 4 5 5 5 5 3 3
$thresholds
[1] -6.035540 -4.348335 -1.159517 1.590627 3.733377 10.819638sm.clustering.extrema(x,30)
$cluster
[1] 2 3 2 3 2 2 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 2 2 3 3 3 3 3 2 3 3 2 3 2 3 3
[38] 3 3 3 2 3 2 2 3 3 3 2 2 3 3 2 3 3 3 3 2 3 3 3 2 2 3 3 3 3 2 3 3 3 2 3 3 3
[75] 3 2 3 2 3 3 3 3 2 3 3 2 3 3 3 3 3 2 3 1 3 2 3 3 3 3 3 3 3 2 3 3 3 3 3 2 3
[112] 3 3 3 3 3 3 2 3 3 3 2 3 3 3 3 2 2 3 2 2 3 3 3 3 2 2 3 3 3 2 3 2 3 3 2 2 3
[149] 2 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 2 3 3 2 3 3 3 3 3 2 3 2 2 2
[186] 2 3 2 3 3 3 3 3 3 3 3 3 3 2 2 2 3 3 3 2 3 3 3 3 3 3 2 2 3 3 3 3 2 3 3 3 3
[223] 2 3 3 3 2 3 2 3 2 3 3 3 2 2 3 2 3 3 2 3 2 3 3 3 2 3 3 3 3 1 3 3 3 3 2 3 3
[260] 3 2 3 3 2 3 2 3 3 3 3 2 3 3 3 3 3 2 3 3 3 2 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3
[297] 2 2 3 3 2 3 2 3 3 3 2 3 3 3 2 3 3 2 2 2 3 3 3 2 3 3 3 3 2 2 3 3 2 3 2 3 3
[334] 3 3 3 3 3 2 3 3 3 3 3 2 2 2 3 2 2 3 3 3 2 3 2 2 3 2 2 3 3 3 3 3 3 3 2 2 3
[371] 3 3 3 3 2 2 3 3 3 1 3 2 3 2 3 2 3 3 2 3 2 3 3 2 3 3 3 3 2 2
$thresholds
[1] -6.035540 -4.348335 3.733377 10.819638sm.clustering.extrema(x,70)$cluster
[1] 2 3 2 3 2 2 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 2 2 3 3 3 3 3 2 3 3 2 3 2 3 3
[38] 3 3 3 2 3 2 2 3 3 3 2 2 3 3 2 3 3 3 3 2 3 3 3 2 2 3 3 3 3 2 3 3 3 2 3 3 3
[75] 3 2 3 2 3 3 3 3 2 3 3 2 3 3 3 3 3 2 3 1 3 2 3 3 3 3 3 3 3 2 3 3 3 3 3 2 3
[112] 3 3 3 3 3 3 2 3 3 3 2 3 3 3 3 2 2 3 2 2 3 3 3 3 2 2 3 3 3 2 3 2 3 3 2 2 3
[149] 2 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 2 3 3 2 3 3 3 3 3 2 3 2 2 2
[186] 2 3 2 3 3 3 3 3 3 3 3 3 3 2 2 2 3 3 3 2 3 3 3 3 3 3 2 2 3 3 3 3 2 3 3 3 3
[223] 2 3 3 3 2 3 2 3 2 3 3 3 2 2 3 2 3 3 2 3 2 3 3 3 2 3 3 3 3 1 3 3 3 3 2 3 3
[260] 3 2 3 3 2 3 2 3 3 3 3 2 3 3 3 3 3 2 3 3 3 2 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3
[297] 2 2 3 3 2 3 2 3 3 3 2 3 3 3 2 3 3 2 2 2 3 3 3 2 3 3 3 3 2 2 3 3 2 3 2 3 3
[334] 3 3 3 3 3 2 3 3 3 3 3 2 2 2 3 2 2 3 3 3 2 3 2 2 3 2 2 3 3 3 3 3 3 3 2 2 3
[371] 3 3 3 3 2 2 3 3 3 1 3 2 3 2 3 2 3 3 2 3 2 3 3 2 3 3 3 3 2 2
$thresholds
[1] -6.035540 -4.348335 3.733377 10.819638On Galaxies and Faithful
$cluster
[1] 1 1 1 1 1 1 1 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[39] 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[77] 4 5 5 6 6 7
$thresholds
[1] 9172.00 13243.41 17164.02 21486.74 26236.72 29529.03 33726.09 34279.00
$cluster
[1] 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 1 2 1 2 1 1 2 2 2 2 1 2 2 2 2 2 2 2 2 1 1
[38] 2 1 2 2 1 2 1 2 2 2 1 2 1 2 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 2 2 1 2 2 1 2 2
[75] 1 2 1 2 2 2 2 2 2 1 2 2 2 2 1 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2 2 1 2 1 2 2 2
[112] 1 2 2 1 2 1 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 2 2 1 2 1
[149] 2 1 2 2 1 2 2 2 2 2 1 2 1 2 1 2 2 2 1 2 1 2 1 1 2 2 2 2 2 1 2 2 1 2 2 2 1
[186] 2 2 1 2 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2
[223] 1 2 2 2 2 2 2 2 2 1 2 1 2 1 1 2 2 1 2 1 2 1 2 2 1 2 1 2 1 2 2 2 2 2 2 2 1
[260] 2 2 2 1 2 1 1 2 2 1 2 1 2
$thresholds
[1] 1.600000 3.053954 5.100000
$cluster
[1] 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 1 2 1 2 1 1 2 2 2 2 1 2 2 2 2 2 1 2 2 1 1
[38] 2 1 2 2 1 2 1 2 2 1 1 2 1 2 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 2 2 1 2 2 1 2 2
[75] 1 2 1 2 2 2 2 2 2 1 2 2 2 2 1 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2 2 1 2 1 2 2 2
[112] 1 2 2 1 2 1 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 2 2 1 2 1
[149] 2 1 2 2 1 2 2 2 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 2 2 2 2 2 1 2 2 1 2 2 2 1
[186] 2 2 1 2 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1 2 2 1 2 2 2 1 2 1 2 1 2 1 2 1 2
[223] 1 2 2 2 2 2 2 2 2 1 2 1 2 1 1 2 2 1 2 1 2 1 2 2 1 2 2 2 1 2 2 2 2 2 2 2 1
[260] 2 2 2 1 2 1 1 2 2 1 2 1 2
$thresholds
[1] 43.00000 66.92693 96.000002. Regression estimation
2.1 Experiments
data("pressure")
t=pressure$temperature
p=pressure$pressure
plot(t,p)
res=sm.regression(t,p,method='cv',nbins=0)
names(res) [1] "eval.points" "estimate" "model.y" "se" "sigma"
[6] "h" "hweights" "weights" "data" "call" res$h[1] 7.034694res$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.000000res$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+02n=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,nbins=0)
res=sm.regression(x,y,h=.07,nbins=0)
res=sm.regression(x,y,h=.7,nbins=0)
par(oma=c(0,0,0,7))
res=sm.regression(x,y,method='cv',nbins=0,col='orange',lty=2)
z=seq(.1,3,len=1000)
lines(z,10*exp(-3*z)+7*cos(2*pi*z)/sqrt(z),col='red',lty=1)
sm.regression(x,y,add=T,col='green',nbins=0,lty=2)
sm.regression(x,y,add=T,col='blue',method='aicc',nbins=0,lty=2)
legend(par('usr')[2], par('usr')[4],xpd=NA,#"topright", inset=c(-0.2,0),
legend = c('r','df','cv','aicc'),
col = c('red','green','orange','blue'), lty=c(1,rep(2,3)))
names(res) [1] "eval.points" "estimate" "model.y" "se" "sigma"
[6] "h" "hweights" "weights" "data" "call" res$h[1] 0.05355961res$eval.points [1] 0.1010678 0.1602235 0.2193792 0.2785350 0.3376907 0.3968464 0.4560022
[8] 0.5151579 0.5743137 0.6334694 0.6926251 0.7517809 0.8109366 0.8700923
[15] 0.9292481 0.9884038 1.0475595 1.1067153 1.1658710 1.2250267 1.2841825
[22] 1.3433382 1.4024939 1.4616497 1.5208054 1.5799611 1.6391169 1.6982726
[29] 1.7574283 1.8165841 1.8757398 1.9348955 1.9940513 2.0532070 2.1123627
[36] 2.1715185 2.2306742 2.2898299 2.3489857 2.4081414 2.4672971 2.5264529
[43] 2.5856086 2.6447643 2.7039201 2.7630758 2.8222315 2.8813873 2.9405430
[50] 2.9996987res$estimate [1] 25.76365944 16.41472439 8.75007624 2.60032010 -2.44060198 -5.25041925
[7] -7.19543126 -8.05636656 -6.45478626 -3.99495107 -2.09935943 0.12676566
[13] 2.41614738 4.76854043 7.20270490 8.14290601 7.89032455 6.31457105
[19] 3.76977648 1.08492548 -0.80589254 -2.47060584 -4.01434706 -4.82256280
[25] -4.84263442 -4.37380002 -3.14254616 -1.39608035 0.04764027 1.58887859
[31] 2.94870347 4.09447862 3.90115457 3.50094976 3.42317728 2.38727191
[37] 0.14876600 -2.10504013 -4.16458105 -5.43560815 -4.87167886 -4.57263781
[43] -4.33724342 -3.42142928 -0.99378872 1.06338075 2.25123545 2.59044737
[49] 3.54531340 4.95102849res=sm.regression(x,y,method='cv',eval.points=sort(x),nbins=0)
plot(res$estimate,y[order(x)])
abline(0,1)
2.2 Compute CV
Compute_CV=function(x,y,method,weights=NA)
{
hopt=sm.regression(x,y,method=method,display='none',nbins=0)$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',nbins=0)$estimate)^2*weights[i]
}
return(CV)
}
Compute_CV(x,y,'cv')[1] 1825.422.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.93546test [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
ROC=function(z,p1,plot=FALSE)
{
p1=pmin(pmax(p1,0),1)
p=c(-0.0001,sort(unique(p1)),1.0001)#
#p=seq(-0.01,1.01,0.01)#
#c(-0.0001,sort(unique(p1))-0.00000000001,1.0001)
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]))
AUC=sum((ROC[-(lp),2]+ROC[-1,2])*(ROC[-(lp),1]-ROC[-(1),1])/2)
colnames(ROC)=c("FPR","TPR")
if (plot)
{
plot(ROC,type='l',main=paste('AUC =', AUC),xlim=c(0,1),ylim=c(0,1))
abline(0,1,col='orange')
}
return(list(ROC=ROC,AUC=AUC))
}Estimated Bayes estimator
#######
Classif_NP=function(X,Y,X0=NA,plot=FALSE)
{
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)
colnames(P0)=V
Class0=rep(NA,n)
}
for (v in 1:n_V)
{
z=as.numeric(Y==V[v])
for (i in 1:n )
{
value=sm.regression(X[-i,],z[-i], eval.points=X,eval.grid=FALSE,display='none',method='cv',nbins=0)$estimate[i]
if (!is.na(value)) {Prob[i,v]=value}
if (is.na(value)) {Prob[i,v]=mean(z[-i])}
}
if (!is.na(max(X0))) {P0[,v]=sm.regression(X,z,eval.points=X0,eval.grid=FALSE,display='none',method='cv',nbins=0)$estimate
P0[is.na(P0[,v]),v]=mean(z)}
}
if (n_V==2) {Roc=ROC(Y==V[2],Prob[,2],plot)}
Class=V[apply(Prob,1,which.max)]
V_est=sort(unique(Class))
if (length(V_est)==n_V){M_table=table(Y,Class)}
else {
M_table=matrix(0,n_V,n_V)
M_table0=table(Y,Class)
for (j in 1:length(V_est)) {M_table[,which(V==V_est[j])]=M_table0[,j]}
}
Err=1-(sum(diag(M_table))/sum(M_table))
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=Err, Class0=Class0,Prob0=P0,Auc=ifelse(n_V==2,Roc$AUC,NA)))}
else {return(list(Class=Class, Prob=Prob, M_table=M_table, Err=Err,Auc=ifelse(n_V==2,Roc$AUC,NA)))}
}
(R=Classif_NP(hema,test,plot=TRUE))
$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.00154507 -0.0015450745
[2,] -0.04382121 1.0438212062
[3,] 0.15434486 0.8456551429
[4,] 0.90091255 0.0990874512
[5,] -0.01638242 1.0163824152
[6,] -0.01690106 1.0169010600
[7,] 0.81826853 0.1817314683
[8,] 0.44026871 0.5597312871
[9,] 1.00271923 -0.0027192346
[10,] 0.97256267 0.0274373255
[11,] 0.99912319 0.0008768104
[12,] 0.08722628 0.9127737158
[13,] 0.74966874 0.2503312641
[14,] 0.99257122 0.0074287785
[15,] 0.96949129 0.0305087136
[16,] 0.85954431 0.1404556941
[17,] 0.98031247 0.0196875314
[18,] 0.90954426 0.0904557440
[19,] 1.01159987 -0.0115998736
[20,] 1.01219351 -0.0121935133
[21,] 0.84465421 0.1553457920
[22,] 0.91948677 0.0805132308
[23,] 0.87502350 0.1249764977
[24,] 1.01244526 -0.0124452622
[25,] 0.37195543 0.6280445651
[26,] -0.03976347 1.0397634697
[27,] 0.93550456 0.0644954388
[28,] 0.70593566 0.2940643421
[29,] 0.16263064 0.8373693607
[30,] -0.10987017 1.1098701736
[31,] 0.83712560 0.1628744014
[32,] 0.93883629 0.0611637069
[33,] 0.98622358 0.0137764204
[34,] 0.96902890 0.0309710958
[35,] 0.74592644 0.2540735614
[36,] 0.12494026 0.8750597443
[37,] 0.58845052 0.4115494755
[38,] 0.95204114 0.0479588551
[39,] 0.23653296 0.7634670375
[40,] 1.01233017 -0.0123301720
[41,] 0.33199241 0.6680075853
[42,] 0.99507654 0.0049234647
[43,] -0.01637771 1.0163777104
[44,] 0.90040899 0.0995910059
[45,] 0.03615177 0.9638482341
[46,] 0.03197847 0.9680215341
[47,] 0.92364019 0.0763598122
[48,] 0.01908133 0.9809186730
[49,] 0.92182592 0.0781740837
[50,] 0.10130621 0.8986937855
[51,] 1.00314294 -0.0031429360
[52,] 0.97278976 0.0272102447
[53,] 0.76563576 0.2343642391
[54,] 0.94776350 0.0522365039
[55,] 0.86088772 0.1391122774
[56,] 0.63806219 0.3619378147
[57,] 0.83158518 0.1684148177
[58,] 0.52611805 0.4738819468
[59,] 0.97711055 0.0228894489
[60,] -0.07255244 1.0725524359
[61,] 1.00631785 -0.0063178519
[62,] 0.78037033 0.2196296735
[63,] 0.21164314 0.7883568648
[64,] 0.01457321 0.9854267891
[65,] 1.00865833 -0.0086583277
[66,] 0.35111770 0.6488822958
[67,] 0.91193400 0.0880660046
[68,] 0.90284398 0.0971560181
[69,] 0.30071516 0.6992848429
[70,] 0.98798550 0.0120145043
[71,] 1.00131155 -0.0013115545
[72,] 0.70514409 0.2948559066
[73,] 0.69715038 0.3028496153
[74,] 1.00909475 -0.0090947457
[75,] 0.88832274 0.1116772601
$M_table
Class
Y negatif positif
negatif 47 3
positif 5 20
$Err
[1] 0.1066667
$Auc
[1] 0.9656plot(hema,R$Prob[,2])
Classif_NP(hema,test,c(37,42,58,57),TRUE)
$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.00154507 -0.0015450745
[2,] -0.04382121 1.0438212062
[3,] 0.15434486 0.8456551429
[4,] 0.90091255 0.0990874512
[5,] -0.01638242 1.0163824152
[6,] -0.01690106 1.0169010600
[7,] 0.81826853 0.1817314683
[8,] 0.44026871 0.5597312871
[9,] 1.00271923 -0.0027192346
[10,] 0.97256267 0.0274373255
[11,] 0.99912319 0.0008768104
[12,] 0.08722628 0.9127737158
[13,] 0.74966874 0.2503312641
[14,] 0.99257122 0.0074287785
[15,] 0.96949129 0.0305087136
[16,] 0.85954431 0.1404556941
[17,] 0.98031247 0.0196875314
[18,] 0.90954426 0.0904557440
[19,] 1.01159987 -0.0115998736
[20,] 1.01219351 -0.0121935133
[21,] 0.84465421 0.1553457920
[22,] 0.91948677 0.0805132308
[23,] 0.87502350 0.1249764977
[24,] 1.01244526 -0.0124452622
[25,] 0.37195543 0.6280445651
[26,] -0.03976347 1.0397634697
[27,] 0.93550456 0.0644954388
[28,] 0.70593566 0.2940643421
[29,] 0.16263064 0.8373693607
[30,] -0.10987017 1.1098701736
[31,] 0.83712560 0.1628744014
[32,] 0.93883629 0.0611637069
[33,] 0.98622358 0.0137764204
[34,] 0.96902890 0.0309710958
[35,] 0.74592644 0.2540735614
[36,] 0.12494026 0.8750597443
[37,] 0.58845052 0.4115494755
[38,] 0.95204114 0.0479588551
[39,] 0.23653296 0.7634670375
[40,] 1.01233017 -0.0123301720
[41,] 0.33199241 0.6680075853
[42,] 0.99507654 0.0049234647
[43,] -0.01637771 1.0163777104
[44,] 0.90040899 0.0995910059
[45,] 0.03615177 0.9638482341
[46,] 0.03197847 0.9680215341
[47,] 0.92364019 0.0763598122
[48,] 0.01908133 0.9809186730
[49,] 0.92182592 0.0781740837
[50,] 0.10130621 0.8986937855
[51,] 1.00314294 -0.0031429360
[52,] 0.97278976 0.0272102447
[53,] 0.76563576 0.2343642391
[54,] 0.94776350 0.0522365039
[55,] 0.86088772 0.1391122774
[56,] 0.63806219 0.3619378147
[57,] 0.83158518 0.1684148177
[58,] 0.52611805 0.4738819468
[59,] 0.97711055 0.0228894489
[60,] -0.07255244 1.0725524359
[61,] 1.00631785 -0.0063178519
[62,] 0.78037033 0.2196296735
[63,] 0.21164314 0.7883568648
[64,] 0.01457321 0.9854267891
[65,] 1.00865833 -0.0086583277
[66,] 0.35111770 0.6488822958
[67,] 0.91193400 0.0880660046
[68,] 0.90284398 0.0971560181
[69,] 0.30071516 0.6992848429
[70,] 0.98798550 0.0120145043
[71,] 1.00131155 -0.0013115545
[72,] 0.70514409 0.2948559066
[73,] 0.69715038 0.3028496153
[74,] 1.00909475 -0.0090947457
[75,] 0.88832274 0.1116772601
$M_table
Class
Y negatif positif
negatif 47 3
positif 5 20
$Err
[1] 0.1066667
$Class0
[1] "negatif" "negatif" "positif" "positif"
$Prob0
negatif positif
[1,] 1.00181472 -0.001814722
[2,] 1.01045534 -0.010455342
[3,] -0.04880557 1.048805567
[4,] -0.01445100 1.014450998
$Auc
[1] 0.96563. Multivariate and Functional case
3.1 Multivariate case
3.1.1 Regression
load('Randonnee.RData')
library(knitr)
library(kableExtra)
kable(cbind(long,lat,alti),'html',caption="Data given in Randonnee.RData") %>%
kable_styling() %>%
scroll_box(width = "100%", height='7cm')| long | lat | alti |
|---|---|---|
| 99.0624679 | 53.3076627 | 1503.552 |
| -6.8700422 | 28.2190488 | 1565.064 |
| -3.9917884 | 92.2394410 | 1523.288 |
| 31.8887100 | -9.1710346 | 1495.187 |
| 20.8955966 | 25.1860652 | 1612.269 |
| 15.3682732 | 80.5044866 | 1511.161 |
| 11.4690357 | 88.5174201 | 1503.051 |
| 79.8319344 | 102.4300348 | 1535.229 |
| 56.0075673 | 107.6692620 | 1503.432 |
| 29.4260498 | 16.7352382 | 1573.805 |
| -4.3391908 | 102.0056042 | 1478.104 |
| 72.7352515 | 32.1038612 | 1557.071 |
| 80.2107291 | 0.1943822 | 1509.054 |
| 13.0087572 | 93.5295195 | 1519.488 |
| 17.4792537 | -7.7909738 | 1517.155 |
| 35.5514487 | 5.6133995 | 1532.295 |
| 13.5712102 | 25.8517935 | 1614.147 |
| 2.0481796 | -7.4330306 | 1508.660 |
| 22.9454398 | 27.9979667 | 1606.543 |
| 36.9931913 | 31.7602686 | 1551.006 |
| 17.0726327 | 80.8224924 | 1498.533 |
| 97.9897970 | 77.7829623 | 1547.869 |
| 100.2870167 | 61.8244704 | 1510.328 |
| 105.5622748 | 106.6027542 | 1507.775 |
| 21.9907795 | 12.9867447 | 1559.080 |
| 46.8235776 | 29.8245296 | 1565.675 |
| 45.4129347 | 61.9143907 | 1490.260 |
| 34.3569008 | 93.7161598 | 1505.586 |
| 35.6348959 | -3.0641687 | 1517.431 |
| 42.1465047 | 90.6220438 | 1496.808 |
| 54.9246438 | 71.5075653 | 1522.724 |
| -0.6032085 | 107.4425999 | 1503.936 |
| 38.9608041 | 79.0489369 | 1504.162 |
| 107.5811816 | 23.9735269 | 1519.739 |
| 52.8454450 | 100.7978796 | 1512.662 |
| 50.1340882 | -0.9997788 | 1531.310 |
| 7.7555722 | 32.0321676 | 1648.128 |
| -5.8852304 | 3.6869886 | 1522.290 |
| 105.2549906 | 78.7011588 | 1518.597 |
| 47.1467618 | 66.8373312 | 1522.303 |
| -4.6626037 | 92.5977188 | 1479.269 |
| 8.7576174 | 34.3697844 | 1641.704 |
| 52.2500823 | 13.6763844 | 1536.718 |
| -7.6250490 | -0.9289558 | 1501.063 |
| 36.5993508 | 54.0141936 | 1499.579 |
| 92.0424269 | 84.5729018 | 1555.299 |
| 103.4379530 | 33.0318845 | 1510.353 |
| 50.7240071 | 2.2425787 | 1509.026 |
| 97.4716941 | 77.9736253 | 1515.126 |
| 73.7016949 | 99.7308117 | 1536.054 |
| 93.0219061 | 67.1765613 | 1547.277 |
| 65.6176874 | 98.6371979 | 1526.580 |
| -1.3699311 | 73.6878124 | 1495.551 |
| 52.3760663 | 78.5963519 | 1530.826 |
| 84.7114132 | 42.2642248 | 1495.388 |
| 107.5771911 | 48.0871036 | 1506.657 |
| 64.4604991 | 27.6326412 | 1591.445 |
| 59.9686836 | 51.0937702 | 1505.332 |
| 8.4859424 | 75.3498196 | 1493.466 |
| 60.8453952 | 81.6616607 | 1542.620 |
| 32.6244900 | 26.3138900 | 1555.251 |
| 109.0847392 | 29.1766504 | 1516.603 |
| 72.0450547 | 19.7782183 | 1586.880 |
| 62.9426072 | 27.6257777 | 1591.079 |
| 49.7946765 | 51.4461836 | 1512.832 |
| -3.7844638 | 108.6813870 | 1496.539 |
| 14.6025297 | 102.8902643 | 1511.974 |
| 68.9487006 | 98.7188453 | 1524.002 |
| 19.4557040 | 103.9766379 | 1498.990 |
| 23.4424702 | 91.6052843 | 1507.005 |
| 96.9559850 | 17.2084007 | 1540.036 |
| 17.4085176 | 16.5421165 | 1569.178 |
| 7.6359436 | 93.0814601 | 1499.439 |
| 105.9557061 | 88.8294518 | 1510.976 |
| -0.6395778 | 27.0807526 | 1583.589 |
| 40.0386443 | 95.0945388 | 1494.300 |
| 108.1165074 | 91.5913831 | 1499.777 |
| 37.4706876 | 36.5362536 | 1548.957 |
| 85.6727815 | 15.5733686 | 1548.695 |
| 83.4007558 | 32.6764199 | 1546.176 |
| 34.5829461 | 80.8233364 | 1503.971 |
| 74.1981223 | 56.3901347 | 1552.938 |
| 65.7956928 | 21.5336391 | 1590.259 |
| 63.1200929 | 44.2666265 | 1501.393 |
| 54.6831032 | 75.6658447 | 1534.508 |
| 25.2570909 | 8.4122134 | 1540.408 |
| 109.0367705 | -4.0753386 | 1503.623 |
| 45.4518286 | 7.7297729 | 1527.511 |
| 77.6834594 | 45.9806080 | 1499.680 |
| 40.0775024 | 78.3242334 | 1501.483 |
| 46.8420351 | 9.7027086 | 1520.665 |
| 67.0032970 | 47.9183023 | 1531.010 |
| 62.9415996 | 28.0067210 | 1580.939 |
| 4.6844503 | 59.6427374 | 1497.629 |
| 12.4587184 | 55.0600018 | 1535.396 |
| 97.3825485 | 59.4960008 | 1513.649 |
| 0.6584921 | 90.8339782 | 1488.929 |
| 34.1886764 | 36.2218561 | 1548.870 |
| -3.9022830 | -6.3884918 | 1494.956 |
| 15.7665207 | 64.0273335 | 1494.208 |
| 109.2956715 | 107.3541429 | 1498.673 |
| 103.5520772 | 6.5037241 | 1505.485 |
| 73.3759137 | 22.8530218 | 1614.789 |
| 108.1589985 | 64.7802229 | 1505.992 |
| 38.7293840 | 58.6801093 | 1505.095 |
| 26.9991154 | 52.6071416 | 1543.953 |
| 93.2730551 | 14.1047208 | 1539.745 |
| 74.6520478 | 74.5988890 | 1719.653 |
| -3.8600043 | 32.2399566 | 1569.702 |
| 91.0284812 | 61.8065544 | 1530.437 |
| 39.7020352 | 44.2541240 | 1507.814 |
| 89.3797159 | -0.0059296 | 1513.673 |
| 8.9815342 | 12.6989354 | 1558.054 |
| 37.3474934 | 83.3984831 | 1521.442 |
| 82.5426909 | 72.2339545 | 1730.754 |
| 41.5408641 | -3.0519692 | 1507.140 |
| 69.1754561 | 49.5537411 | 1536.695 |
| 107.2183426 | 59.4791128 | 1494.653 |
| 95.5583194 | 36.3866017 | 1522.613 |
| 45.1385816 | 32.8278674 | 1529.149 |
| 71.4557269 | 53.9506069 | 1527.618 |
| 0.6201439 | 52.5379462 | 1534.111 |
| 86.4482111 | -1.9838803 | 1485.289 |
| 49.4328418 | 57.9250941 | 1511.933 |
| 8.9252843 | -6.5767055 | 1510.112 |
| 48.3882847 | 15.9442132 | 1544.139 |
| 101.8168975 | 85.5506530 | 1505.108 |
| 73.1332982 | 36.1948201 | 1537.994 |
| 78.2796414 | 75.6409429 | 1763.630 |
| 33.9245340 | 69.5515850 | 1504.757 |
| 64.2372093 | 101.3554996 | 1507.480 |
| 54.9992699 | 56.6553987 | 1510.998 |
| 34.5305953 | 9.1060737 | 1533.052 |
| 86.1025130 | 53.3294582 | 1527.883 |
| 17.8944715 | 18.5229879 | 1591.969 |
| 65.3099775 | 92.7764085 | 1540.220 |
| 68.3107806 | 78.4274628 | 1609.276 |
| -9.1280097 | 74.9568642 | 1500.164 |
| -8.4378068 | 7.7611201 | 1515.473 |
| -9.4375850 | 58.2351528 | 1500.782 |
| 72.9383652 | 9.1428750 | 1545.672 |
| 3.1518398 | 44.8344709 | 1593.204 |
| -4.5201100 | 82.0173986 | 1506.172 |
| 85.1011407 | 94.7242777 | 1540.134 |
| 58.4218492 | 85.5438200 | 1514.645 |
| 21.5381986 | 84.0624910 | 1490.193 |
| 79.5683257 | 107.7563988 | 1519.539 |
| 16.8094861 | 0.1065805 | 1518.925 |
| 13.7168195 | 52.0209065 | 1557.457 |
| 61.6293227 | 22.2596313 | 1603.621 |
| 68.8874839 | 72.1662832 | 1591.796 |
| 31.9630115 | 53.6803757 | 1522.681 |
| -2.4530409 | 104.8597823 | 1496.549 |
| 67.5898129 | 64.8423330 | 1574.086 |
| -7.8298760 | 86.5592916 | 1505.518 |
| 106.1900164 | 109.1123272 | 1505.576 |
| 30.3883853 | 70.1951086 | 1501.087 |
| 37.8230056 | 88.7199266 | 1506.887 |
| 85.1603546 | 60.2615274 | 1588.606 |
| 11.7026249 | 76.6645402 | 1495.859 |
| -1.5148323 | 36.7445473 | 1582.416 |
| 80.6039050 | 108.0380414 | 1492.412 |
| 42.9912725 | 23.7753807 | 1542.679 |
| 65.6447334 | 51.4841565 | 1528.102 |
| 31.0267306 | 82.7948433 | 1496.188 |
| 15.2377335 | 60.5768411 | 1516.101 |
| 35.1865419 | 32.5634608 | 1557.464 |
| 59.5301555 | 44.2504127 | 1521.483 |
| 93.0327959 | 66.2506154 | 1540.498 |
| 19.7201658 | 20.2465274 | 1598.686 |
| 4.4496954 | 60.1801283 | 1517.926 |
| 43.5310138 | 89.3777994 | 1512.035 |
| 57.3573328 | 66.7790963 | 1515.571 |
| 96.1369763 | 99.0289473 | 1507.688 |
| 53.9646099 | 38.8948529 | 1522.677 |
| 6.4655638 | 55.8074369 | 1521.504 |
| 109.9032347 | 14.6825246 | 1513.087 |
| 57.0984624 | 80.8084841 | 1528.505 |
| 108.2334729 | -7.0929847 | 1508.686 |
| 66.5735029 | 16.4661541 | 1569.290 |
| 26.3783911 | 74.6614091 | 1489.932 |
| 95.4501844 | 92.9089064 | 1525.671 |
| 12.1398498 | -1.6929547 | 1516.118 |
| 61.3491014 | -3.2044426 | 1497.463 |
| 97.0415713 | 92.2378615 | 1507.057 |
| 67.8564440 | 57.2680875 | 1536.328 |
| 94.5919541 | 17.4978882 | 1538.884 |
| 107.7176833 | 48.5589143 | 1500.362 |
| 86.2763029 | 52.1070050 | 1527.315 |
| 109.9107141 | 85.6347466 | 1509.714 |
| 82.0665086 | 102.9683853 | 1507.577 |
| 74.3138162 | 72.2160242 | 1685.401 |
| 32.8671800 | 107.6644860 | 1503.280 |
| 74.6963442 | 26.8117389 | 1586.010 |
| 24.5364777 | 75.7511033 | 1495.289 |
| 64.5714931 | 40.3632850 | 1511.456 |
| 29.4825232 | 36.9292619 | 1589.894 |
| 42.2030134 | 95.1776168 | 1505.560 |
| 32.9716050 | 12.2117467 | 1553.082 |
| 47.3855472 | 74.9843602 | 1537.770 |
| 32.9537145 | 16.3360995 | 1574.133 |
| 19.7007396 | 86.4568328 | 1487.686 |
| 82.3053686 | 85.2028623 | 1698.334 |
| 29.7235724 | 66.4521574 | 1509.959 |
| 41.6190261 | 99.0089405 | 1511.693 |
| 0.4811028 | 77.2090162 | 1503.002 |
| 33.8853365 | 31.3755255 | 1569.472 |
| 69.1828809 | -0.2142083 | 1484.409 |
| -8.5955022 | 68.7260890 | 1498.053 |
| 83.8354765 | 80.6439016 | 1710.081 |
| 38.8785215 | 81.4825360 | 1509.792 |
| 101.6520975 | 95.0324957 | 1500.674 |
| 38.2653881 | 15.0770843 | 1554.486 |
| 92.0097578 | 79.4081350 | 1580.155 |
| 48.5082542 | 7.1348003 | 1499.173 |
| 37.4501140 | 0.8510836 | 1508.732 |
| 97.5671286 | 36.7153758 | 1514.219 |
| 32.2246511 | 47.0856551 | 1547.203 |
| 27.6389411 | -4.2890589 | 1510.215 |
| -9.0511564 | 93.6098017 | 1493.465 |
| 77.5580233 | 39.0210219 | 1550.051 |
| 94.9180675 | 40.9643174 | 1498.312 |
| 58.0855414 | 86.3287291 | 1528.816 |
| 41.4742304 | 85.9292547 | 1518.739 |
| -4.9285801 | 19.5426376 | 1537.191 |
| 12.5128935 | 54.1810474 | 1522.310 |
| 98.8614764 | 29.0301088 | 1512.817 |
| -2.2395305 | 13.0767943 | 1530.683 |
| -7.8879634 | 96.4831019 | 1502.368 |
| 106.0577077 | 29.4490712 | 1521.374 |
| 56.8031475 | 103.1553699 | 1519.785 |
| 3.2581064 | 16.1983463 | 1553.751 |
| 55.3891077 | 81.7767231 | 1514.695 |
| 109.5003712 | 34.9192071 | 1507.356 |
| 103.9531725 | 2.6824025 | 1515.738 |
| 90.9163884 | 64.6651965 | 1565.437 |
| 17.5675969 | 40.4335865 | 1605.220 |
| 56.6706462 | 53.7304172 | 1530.561 |
| 103.7406261 | 108.0123845 | 1510.335 |
| 83.7976525 | 10.6557069 | 1517.980 |
| 36.3899025 | 58.6566456 | 1503.187 |
| 41.1020887 | 37.8390720 | 1530.025 |
| -8.1780918 | 100.2869428 | 1519.485 |
| 11.0855115 | 82.6831211 | 1510.155 |
| 78.5711419 | 67.2791163 | 1704.259 |
| 10.5236860 | 19.0974620 | 1572.168 |
| 87.2162916 | 51.8274711 | 1526.432 |
| 88.9293157 | 85.9404894 | 1574.387 |
| 10.9053526 | 68.4225406 | 1505.185 |
| 9.4655780 | 81.5465188 | 1491.694 |
| 49.8283107 | 108.2866624 | 1508.596 |
| -9.5947153 | 92.2925736 | 1499.622 |
| 84.1258040 | 96.2783463 | 1542.454 |
| 43.7018672 | 44.9876156 | 1516.929 |
| 39.2791223 | 80.1971819 | 1513.394 |
| 34.7361794 | 19.1695937 | 1561.411 |
| 59.4643598 | 34.0714970 | 1570.119 |
| 54.2132162 | -0.7885425 | 1527.899 |
| 95.8608577 | 69.9028278 | 1516.370 |
| 79.3937709 | 8.9799858 | 1507.082 |
| -3.4154369 | 74.7767093 | 1503.140 |
| 102.3425618 | 18.2671573 | 1507.488 |
| 42.2433895 | 75.0631941 | 1512.044 |
| 80.3533921 | 12.8802847 | 1550.120 |
| -2.8782635 | -1.4969069 | 1519.303 |
| 95.1028680 | 15.3441908 | 1528.015 |
| 53.0508449 | 67.9876337 | 1532.086 |
| 75.5722174 | 3.9519753 | 1502.610 |
| 92.1211206 | 17.0929100 | 1532.950 |
| 66.5940911 | 61.7554280 | 1547.222 |
| 68.3343318 | 24.4713866 | 1617.225 |
| 89.6319729 | 92.8701326 | 1535.141 |
| 100.8202775 | 87.0977263 | 1507.188 |
| 21.3615679 | 33.1135690 | 1628.742 |
| 88.6541885 | 39.6670205 | 1526.841 |
| 108.5622098 | 74.8217273 | 1485.196 |
| 56.0671732 | 2.7101703 | 1484.626 |
| 69.1604764 | 34.1565359 | 1550.667 |
| 63.3533009 | 91.0166136 | 1529.457 |
| -1.9578196 | 98.7146702 | 1502.354 |
| 81.3425835 | 21.5284310 | 1589.093 |
| 44.2450981 | 57.6431321 | 1514.573 |
| 88.3049812 | -8.4029509 | 1487.451 |
| 98.2021578 | 45.7207011 | 1481.716 |
| 11.6492440 | 64.9348700 | 1507.737 |
| 5.7902776 | 24.8026485 | 1599.414 |
| 82.7875774 | 59.4499204 | 1592.806 |
| 74.1220887 | 2.9185357 | 1502.794 |
| 32.0810399 | 57.3191678 | 1507.491 |
| 19.2687937 | -0.0460234 | 1525.650 |
| 24.6226843 | 23.3138199 | 1599.437 |
| 66.6141810 | 17.6189217 | 1577.792 |
| 35.0466244 | 5.9158673 | 1538.361 |
| 61.7549079 | 39.3237693 | 1528.888 |
| 89.1798913 | 57.3417052 | 1548.715 |
| 68.3177293 | 3.7032293 | 1525.710 |
| 30.9856599 | 101.5114146 | 1496.134 |
| -5.1436744 | 21.6430878 | 1556.844 |
| -6.5555396 | 62.7678419 | 1513.506 |
| 9.8439612 | 49.0257246 | 1556.732 |
| 31.3041346 | 79.5129649 | 1504.517 |
| 2.0242539 | 11.9979869 | 1539.037 |
| 84.7063420 | 13.1007343 | 1557.251 |
| 29.3631448 | -8.1769586 | 1503.023 |
| 67.0596468 | 54.3304944 | 1524.949 |
| 93.6544835 | 20.8858303 | 1522.962 |
| 95.7701951 | 109.8266496 | 1501.593 |
| 20.8136802 | 9.6575482 | 1552.333 |
| 89.8029411 | -3.7470808 | 1498.926 |
| 89.6190648 | 104.9048206 | 1507.858 |
| 18.5791090 | 94.9981342 | 1514.110 |
| 68.5679271 | 26.6114738 | 1588.468 |
| -2.7874155 | 52.1470068 | 1542.644 |
| 45.4848944 | -1.7112577 | 1523.542 |
| 55.9512635 | 7.8452460 | 1506.643 |
| 48.4728493 | 103.3569716 | 1502.020 |
| 15.6096432 | 72.9310911 | 1503.216 |
| 77.1553687 | 32.2952721 | 1586.116 |
| 101.7259463 | 90.9876603 | 1532.920 |
| 91.0791836 | 50.7894269 | 1493.142 |
| 68.2606661 | 82.4318992 | 1604.794 |
| 6.7598813 | 96.3158633 | 1492.373 |
| 107.6118367 | -0.7938033 | 1497.330 |
| 36.7417797 | 67.2738524 | 1499.485 |
| 11.3383602 | 11.4739348 | 1546.386 |
| 56.3474607 | 61.0645675 | 1523.533 |
| 49.9432404 | 44.6066225 | 1512.207 |
| 12.7349585 | 100.0015666 | 1500.048 |
| 34.2485292 | 90.7136576 | 1509.930 |
| -1.1026503 | 29.0980841 | 1590.837 |
| 88.5980615 | 23.1009244 | 1552.859 |
| 93.1525725 | 13.3031028 | 1534.580 |
| 92.8625568 | 90.9390584 | 1524.078 |
| 58.3477348 | 1.2338420 | 1519.394 |
| 50.7536376 | 95.3763274 | 1498.563 |
| 72.8735053 | 21.9861983 | 1617.793 |
| 89.6411641 | 69.1690246 | 1581.261 |
| 25.1111707 | 56.3215340 | 1496.755 |
| 1.5955143 | 86.9238931 | 1492.788 |
| 65.3676301 | 32.1324569 | 1570.713 |
| 61.0056921 | 5.5268771 | 1504.325 |
| 80.1208224 | 6.9373866 | 1519.488 |
| 13.9519384 | -3.6889649 | 1519.616 |
| 91.4471487 | 90.8848947 | 1531.403 |
| 35.9693894 | -0.9320418 | 1509.471 |
| 91.3207917 | 96.3104715 | 1513.320 |
| 51.3916477 | 43.7737725 | 1498.221 |
| 101.1608726 | -2.4532839 | 1502.029 |
| 18.7078511 | 64.4911674 | 1508.166 |
| 50.1847349 | 97.4818228 | 1508.136 |
| -6.5284906 | 65.4102213 | 1495.612 |
| 22.9682395 | 78.3188111 | 1504.898 |
| -1.6257895 | 22.5323168 | 1557.789 |
| 37.6809647 | 86.0740078 | 1505.012 |
| 44.5150829 | -8.0919361 | 1492.447 |
| 98.6483570 | 37.1308551 | 1508.043 |
| -3.8899612 | 63.9407826 | 1510.358 |
| 76.6220504 | 50.6024568 | 1542.991 |
| 34.5582256 | 6.0812392 | 1536.248 |
| 61.5165984 | 100.7611068 | 1505.314 |
| 1.4280745 | -4.7681216 | 1501.668 |
| 9.9571036 | 108.8500554 | 1515.744 |
| 58.1042503 | 53.5467778 | 1530.612 |
| 97.6922949 | 21.7455534 | 1522.981 |
| 57.5637215 | 17.2815492 | 1557.595 |
| 73.2940064 | 86.8986099 | 1636.235 |
| 42.2275999 | 108.1269370 | 1513.672 |
| 79.2007804 | 26.6891237 | 1589.146 |
| 25.0378917 | -7.3017076 | 1501.583 |
| 104.9990713 | 30.7436422 | 1517.994 |
| 60.1676537 | 38.1424053 | 1528.900 |
| 28.1064449 | 103.5866842 | 1501.682 |
| 46.2679811 | 45.6117505 | 1503.142 |
| 0.1100654 | 80.1602309 | 1481.154 |
| 9.7661822 | 20.6649067 | 1583.420 |
| 15.3392095 | 8.7165088 | 1546.996 |
| 86.2068609 | 24.2161395 | 1568.255 |
| 95.9371390 | 87.4804205 | 1505.517 |
| 71.5423485 | 75.8261241 | 1665.263 |
| 85.0595993 | 1.4235339 | 1522.612 |
| 62.7279188 | 35.0926339 | 1524.308 |
| 0.3784555 | 84.3519681 | 1500.357 |
| 109.1100813 | -6.7318973 | 1504.116 |
| 10.1067326 | 8.1089910 | 1523.197 |
| 14.6993165 | 74.4807939 | 1505.789 |
| 19.4556634 | 5.5710324 | 1528.607 |
| 53.8164148 | 40.7650830 | 1529.546 |
| 2.2608089 | -5.2131376 | 1515.577 |
| 69.6403102 | 95.0177066 | 1558.416 |
| 9.2092220 | 49.7994056 | 1572.438 |
| 67.8009084 | 44.0925053 | 1514.321 |
| 84.3480678 | 103.9824313 | 1536.524 |
| 61.4812275 | 62.2715877 | 1520.285 |
| 31.4714729 | -8.9967737 | 1488.522 |
| 48.6934323 | 10.0242483 | 1527.049 |
| 35.9006849 | 109.7142090 | 1502.821 |
| 93.0243156 | -7.7038707 | 1499.589 |
| 82.1181352 | 29.1630460 | 1561.115 |
| 55.5579532 | 34.3340920 | 1528.703 |
| 50.3538442 | 61.4322998 | 1517.249 |
| 65.3020222 | 6.8285881 | 1509.494 |
| 81.8677705 | 72.0445888 | 1734.414 |
| 72.2282170 | 95.8624118 | 1550.477 |
| 97.0193113 | 15.6787341 | 1499.625 |
| 41.4619855 | 24.8536569 | 1563.564 |
| 55.1013032 | 98.4986717 | 1535.853 |
| 40.0393734 | 13.6109056 | 1530.488 |
| -4.0411751 | -4.9687549 | 1517.070 |
| 26.2773770 | 47.8202663 | 1556.910 |
| 81.9814880 | 66.8288662 | 1680.413 |
| 24.8396541 | 77.3473062 | 1506.856 |
| 78.8747559 | -3.0755225 | 1484.717 |
| 26.6468710 | 67.3051887 | 1496.684 |
| -1.6635758 | 91.2182948 | 1504.817 |
| 92.6802950 | 96.2200029 | 1513.405 |
| -2.7809562 | 13.8046413 | 1512.486 |
| 39.1102889 | 2.5987683 | 1502.162 |
| 45.8712546 | 48.8523444 | 1499.620 |
| 89.2561780 | 20.1589514 | 1569.780 |
| -7.0787753 | 73.4842318 | 1504.415 |
| 57.8336239 | 78.1569571 | 1525.848 |
| 74.3947873 | 94.5381534 | 1551.992 |
| 9.4485060 | 94.2673116 | 1496.655 |
| 38.3881533 | 55.6955801 | 1504.315 |
| -6.6127293 | 40.7847500 | 1559.043 |
| 73.7611126 | 103.1028591 | 1527.117 |
| 17.9282801 | 5.0622493 | 1536.839 |
| 53.3900396 | 20.3366755 | 1574.441 |
| 45.4418399 | 73.3952099 | 1509.247 |
| 12.7222182 | 43.8270357 | 1616.504 |
| 77.6123913 | 64.6587233 | 1646.669 |
| 67.8873437 | 81.9926225 | 1582.662 |
| 37.6814209 | 38.7582247 | 1542.598 |
| 32.4833796 | 57.4068465 | 1516.853 |
| 10.6336045 | 41.3802584 | 1618.410 |
| 78.1553770 | 40.5129908 | 1527.939 |
| 103.1758267 | 75.3131802 | 1492.229 |
| -0.0172123 | 93.9257011 | 1509.484 |
| 41.0459263 | 99.9712554 | 1523.821 |
| 109.0236550 | 96.4530513 | 1507.679 |
| 21.8874990 | 101.9317074 | 1497.830 |
| 81.7513145 | -2.9231695 | 1500.372 |
| 99.3454042 | 88.8938747 | 1501.742 |
| 94.9134218 | 57.8260363 | 1527.223 |
| 84.0662564 | 91.0539421 | 1610.124 |
| 109.9569393 | 69.4030961 | 1505.429 |
| 7.6394863 | 83.2413842 | 1507.087 |
| 34.0495309 | 100.4868590 | 1486.270 |
| 28.4740328 | -9.7608886 | 1504.520 |
| 88.9875232 | -6.1349549 | 1504.028 |
| 43.4492684 | 23.5005026 | 1550.172 |
| 96.5806337 | 25.7377867 | 1543.924 |
| 21.5609224 | 86.1965763 | 1496.312 |
| 87.7217712 | 74.6415199 | 1643.099 |
| 15.9316886 | 7.5897939 | 1547.601 |
| 20.5556190 | 24.5636674 | 1609.879 |
| 80.5813421 | 4.1914909 | 1514.706 |
| 88.6072139 | 46.5673720 | 1496.160 |
| 20.1747839 | 102.8409574 | 1491.287 |
| -5.0209855 | 56.1619027 | 1500.766 |
| 90.0875346 | 69.9890892 | 1577.988 |
| 27.6607426 | 23.9461600 | 1577.866 |
| 6.9553972 | 65.0182821 | 1516.048 |
| 34.1058517 | 92.7306515 | 1488.535 |
| 89.2058332 | 11.1390503 | 1520.170 |
| 61.5390024 | 90.9235002 | 1517.024 |
| 59.4918529 | 36.4783816 | 1538.244 |
| 31.9344792 | 16.2839703 | 1560.817 |
| 4.7098450 | 49.7922594 | 1562.877 |
| 53.3536919 | -1.2171213 | 1496.778 |
| 42.9603342 | 100.3160730 | 1504.924 |
| 16.9773379 | -9.3945445 | 1517.950 |
| 63.2456879 | 49.9425281 | 1524.966 |
| 6.0022059 | 100.6571232 | 1503.972 |
| 33.3386119 | 5.3754525 | 1514.375 |
| 97.9236387 | 94.2334175 | 1513.962 |
| 49.4804126 | 35.4222640 | 1526.181 |
| 52.8546455 | 27.6743391 | 1571.022 |
| -5.0341750 | 68.9563156 | 1496.880 |
| 95.7745047 | 82.8733316 | 1502.377 |
| 82.7481690 | 107.2289593 | 1503.532 |
| 42.5722881 | 55.3928997 | 1514.392 |
| 51.5053520 | -4.3927146 | 1483.801 |
| 54.6841293 | 86.4875838 | 1534.524 |
| 70.0131909 | 108.8814492 | 1511.883 |
| 76.0295058 | 29.9439572 | 1587.988 |
| 85.6484853 | 26.1797903 | 1551.291 |
| 97.6697392 | 59.9019719 | 1523.718 |
| 20.8957098 | 0.5052321 | 1519.248 |
| 88.4729862 | 66.2547608 | 1590.906 |
| 20.1538164 | 77.0560286 | 1511.153 |
| 44.8187165 | 51.1533860 | 1518.425 |
| 68.0309262 | -3.8631841 | 1498.178 |
| 58.6289978 | 82.8772579 | 1525.467 |
| 25.7139866 | -3.9269104 | 1522.011 |
| 21.3003076 | 65.2011733 | 1497.728 |
| 62.4479316 | 79.3138901 | 1543.781 |
| 59.2408683 | 79.0786884 | 1548.860 |
| 63.8441639 | 101.3626718 | 1527.701 |
| 39.5953519 | 5.5493366 | 1525.488 |
| 58.4931732 | 23.2969240 | 1605.282 |
| 60.1602868 | 1.9332860 | 1503.454 |
| 1.9114440 | -5.7792464 | 1499.739 |
| 96.2715323 | 46.9828447 | 1494.714 |
| 39.0903538 | 104.7204956 | 1488.624 |
| 101.6557304 | 46.8468954 | 1492.734 |
| 105.0268924 | 37.1232838 | 1502.136 |
| 71.3148199 | 53.0517592 | 1534.068 |
| 24.2461192 | 61.9676012 | 1493.115 |
| 98.1935393 | 82.2192951 | 1544.857 |
| 39.4108665 | 22.5839583 | 1551.701 |
| 3.8308710 | -9.3380757 | 1513.457 |
| 14.3746430 | 27.5417999 | 1624.842 |
| 23.9683261 | 31.4856350 | 1597.986 |
| 70.1198113 | 109.7017905 | 1503.093 |
| 70.8996247 | 87.4279395 | 1612.542 |
| 96.7478936 | 11.5451906 | 1509.057 |
| 98.3181060 | 100.3779429 | 1513.373 |
| 42.5515716 | 23.9775342 | 1559.335 |
| 96.5981999 | 95.1618410 | 1516.647 |
| 41.9065579 | 16.1645199 | 1546.356 |
| 90.8175866 | 94.2590608 | 1538.461 |
| 33.8620463 | 107.5285164 | 1501.119 |
| 6.4959357 | 61.2849425 | 1508.263 |
| 92.4394683 | 77.7916470 | 1575.389 |
| 25.2709471 | 106.8835163 | 1491.833 |
| 20.8148430 | 27.7309045 | 1607.803 |
| 59.1652409 | 54.7402013 | 1518.248 |
| 26.7529195 | 12.0087846 | 1549.326 |
| 17.4847575 | 60.2774071 | 1508.729 |
| 95.9102697 | 71.3994509 | 1516.467 |
| 6.6121522 | 16.3705873 | 1573.759 |
| 36.5492713 | 100.5967313 | 1504.574 |
| 65.4099205 | 37.2396430 | 1546.641 |
| -4.3319132 | 78.5934965 | 1510.538 |
| 22.2900263 | 21.6663253 | 1598.827 |
| 76.1188283 | 95.5553568 | 1582.666 |
| 62.2883352 | 22.4480367 | 1588.908 |
| 38.0815533 | 64.5434269 | 1499.889 |
| 106.6215866 | 23.3670923 | 1513.965 |
| 94.0356429 | 79.3797497 | 1545.206 |
| 66.5149667 | 63.9329539 | 1556.575 |
| 17.1773796 | 25.5647566 | 1611.464 |
| 77.2371694 | 37.3248484 | 1527.759 |
| 75.1846810 | 92.0785114 | 1597.508 |
| 106.7595465 | 64.9866935 | 1498.387 |
| 52.0815477 | 18.3104257 | 1569.869 |
| 76.4454821 | 107.6737522 | 1528.958 |
| 25.2439926 | 73.0578172 | 1507.565 |
| 22.0903061 | 102.0663570 | 1491.416 |
| 9.6399116 | 14.8480610 | 1564.207 |
| 44.1345382 | 90.2871030 | 1500.925 |
| -3.0260081 | 43.1150788 | 1582.300 |
| 68.8843565 | 45.6053514 | 1517.184 |
| 15.9653080 | 75.8507296 | 1492.573 |
| 65.5993460 | -7.6848511 | 1491.313 |
| 26.7613682 | 16.0195621 | 1553.978 |
| 83.1410134 | -4.2698598 | 1501.872 |
| 21.2695025 | 37.2339109 | 1630.816 |
| 70.8444498 | 103.2475674 | 1516.087 |
| 103.9555270 | 65.6625380 | 1515.712 |
| 99.1280174 | 15.2649251 | 1505.829 |
| 2.7211394 | 79.0324802 | 1492.696 |
| 91.0219634 | 45.3216888 | 1500.330 |
| 22.3289261 | 85.5976951 | 1496.914 |
| -1.9080134 | 81.3258636 | 1488.583 |
| -5.7198598 | 72.1514194 | 1500.639 |
| 24.6421009 | 42.2301986 | 1590.441 |
| 30.2130708 | 53.5880265 | 1502.603 |
| 21.0951454 | 15.7748851 | 1562.054 |
| 83.8048238 | 53.9080075 | 1547.451 |
| 35.4662511 | -5.5924372 | 1509.002 |
| 84.4450979 | 31.4639356 | 1568.735 |
| 99.9266260 | 67.1487860 | 1518.432 |
| 46.1240466 | 15.5221679 | 1537.757 |
| 104.7312447 | 12.7352638 | 1515.059 |
| 19.2201695 | 69.6486299 | 1483.650 |
| 96.3367378 | 51.8959265 | 1510.629 |
| 102.1819532 | 50.7179467 | 1495.384 |
| 88.9371270 | 83.3708389 | 1606.308 |
| 46.7660689 | 69.0966059 | 1527.524 |
| 15.9327528 | 92.6088435 | 1489.406 |
| 77.0372848 | 32.9504997 | 1557.608 |
| 58.5118181 | 77.3999484 | 1541.762 |
| 56.9811699 | 53.7657752 | 1500.824 |
| 64.9082567 | -3.2443351 | 1505.016 |
| 38.2952744 | 94.5939577 | 1503.812 |
| 70.8722402 | 29.5239710 | 1598.947 |
| 72.0497879 | 89.0016703 | 1617.575 |
| 101.9334932 | 37.5686142 | 1495.838 |
| 109.8826672 | -2.9751045 | 1515.478 |
| 57.1366042 | 8.0383951 | 1517.843 |
| 12.3675564 | 84.2546153 | 1505.852 |
| 61.6454857 | 75.0934569 | 1543.400 |
| 49.7533377 | 57.2697188 | 1511.848 |
| 71.3972383 | 86.9312615 | 1603.218 |
| 45.5282608 | 65.5898775 | 1518.721 |
| 29.3858278 | 82.2726451 | 1506.451 |
| 22.8850113 | 64.8671459 | 1498.133 |
| 45.0699241 | 24.3136256 | 1569.493 |
| 65.9014948 | 5.8046858 | 1515.366 |
| 65.9690510 | 108.6768581 | 1525.685 |
| 10.7519959 | 86.2096313 | 1515.906 |
| 39.5647257 | 105.8044188 | 1501.214 |
| 20.8036622 | 47.6712879 | 1575.772 |
| 41.9689670 | 9.5280402 | 1520.868 |
| 79.5410652 | 41.5403248 | 1506.560 |
| 26.7119340 | 62.2381561 | 1491.121 |
| 87.0795461 | 107.8822616 | 1508.976 |
| 93.8081755 | 19.4210003 | 1555.853 |
| 91.5327959 | 29.0430997 | 1538.999 |
| 2.7030215 | 18.1596115 | 1567.092 |
| 100.7991998 | 64.0284909 | 1494.948 |
| 39.8600721 | 31.5285579 | 1542.262 |
| 44.1923083 | 55.2737889 | 1502.390 |
| 98.2832157 | 63.1467524 | 1517.818 |
| -3.5277962 | 87.7756079 | 1495.102 |
| 41.0449207 | -2.0839615 | 1503.252 |
| 0.6568744 | -4.1652459 | 1498.135 |
| 103.1660420 | 10.9878643 | 1513.850 |
| 60.3913352 | 62.6394496 | 1521.664 |
| 49.3752225 | 28.5603801 | 1532.891 |
| 81.8617240 | 30.3232149 | 1565.344 |
| 11.6491422 | 85.2267429 | 1500.841 |
| 5.1677096 | 43.6305394 | 1595.774 |
| -3.8947262 | 106.8776215 | 1509.734 |
| 66.4253047 | 105.6924791 | 1510.245 |
| 56.4852030 | 20.1350372 | 1581.431 |
| 77.0702226 | 51.6014682 | 1533.997 |
| -7.2689226 | 88.2718319 | 1511.910 |
| 71.2353067 | -5.4746269 | 1501.750 |
| 31.1148323 | 78.4540884 | 1524.668 |
| 102.6538761 | 11.7575106 | 1511.470 |
| 109.5170704 | 51.2340740 | 1499.621 |
| 42.0342141 | 33.1769020 | 1531.316 |
| 72.8749498 | -6.4039993 | 1497.093 |
| 99.1043105 | 68.9134452 | 1511.582 |
| 91.2494428 | 81.2928025 | 1556.601 |
| 104.7505704 | 105.2019218 | 1498.370 |
| 39.5326475 | 63.9007516 | 1506.227 |
| -9.0300590 | 20.7906424 | 1514.349 |
| 39.6397507 | -7.8658967 | 1488.148 |
| 108.9315627 | 74.0586210 | 1526.546 |
| 54.6685656 | -9.5554745 | 1478.984 |
| 48.0970196 | 2.8355904 | 1500.376 |
| 28.1432877 | 67.9735098 | 1498.586 |
| 22.2442323 | -5.0874047 | 1532.196 |
| 52.6888546 | 21.3399248 | 1556.263 |
| 74.2323138 | 94.1292340 | 1577.932 |
| -4.5672527 | 78.0267880 | 1510.912 |
| 41.2776241 | 17.1847075 | 1546.162 |
| 106.8536079 | 64.6798782 | 1500.722 |
| 94.3077619 | 103.2697906 | 1513.920 |
| 53.5451763 | 40.5375573 | 1517.424 |
| 95.4244042 | -7.5664026 | 1505.027 |
| -4.5798137 | -0.3586452 | 1508.937 |
| 89.6624644 | 99.3101567 | 1518.012 |
| 32.9560425 | 80.4317858 | 1512.044 |
| -1.7280515 | 19.3834340 | 1536.298 |
| 104.1473564 | 92.3759427 | 1500.146 |
| 13.8418379 | -9.2615675 | 1520.128 |
| 10.3614003 | -6.9285953 | 1524.814 |
| 23.2121472 | 102.1444907 | 1489.344 |
| 66.9251651 | 79.4080454 | 1569.973 |
| 88.1680190 | 46.9285504 | 1500.956 |
| 14.9520326 | 13.5341534 | 1558.593 |
| 39.0940108 | 76.8594967 | 1495.935 |
| 15.9254863 | 9.6555981 | 1537.437 |
| 100.0806638 | 92.2636345 | 1513.651 |
| -2.1782721 | 86.6995902 | 1506.084 |
| 54.2987261 | 2.3928148 | 1508.893 |
| 92.0805643 | 20.0592599 | 1534.972 |
| 72.5566401 | 4.6406742 | 1501.272 |
| 72.6624334 | 32.9188668 | 1563.051 |
| 89.5553570 | 77.8546149 | 1613.075 |
| 74.8512809 | 72.6723929 | 1713.797 |
| 47.9141089 | 30.5750335 | 1534.930 |
| 66.0321908 | 40.4925415 | 1515.921 |
| 37.1322801 | 98.1926650 | 1499.387 |
| 50.8429242 | -7.1966780 | 1506.301 |
| -9.2976044 | 70.0467883 | 1496.489 |
| 34.1314216 | -1.0925443 | 1509.968 |
| 83.0657166 | 52.7460320 | 1541.175 |
| 60.6463539 | 87.1843896 | 1531.879 |
| 25.4753125 | -0.3565201 | 1537.802 |
| 63.9951391 | 95.0380162 | 1537.010 |
| 81.1837876 | 3.9920236 | 1493.454 |
| 74.3581717 | 43.8768766 | 1501.473 |
| -3.1426746 | -0.0874256 | 1510.891 |
| 53.4858685 | 12.3895703 | 1527.188 |
| 63.1320570 | 32.2232916 | 1565.811 |
| 93.6906996 | 80.6107477 | 1547.330 |
| 41.8486828 | 89.5659122 | 1508.680 |
| 29.3908050 | 25.6783325 | 1591.965 |
| 97.1850477 | 50.6272007 | 1521.476 |
| 28.1837130 | -6.1465343 | 1520.011 |
| 37.7350850 | 42.0853067 | 1532.732 |
| 56.2257791 | -0.7359922 | 1513.249 |
| 52.2089753 | 92.6217944 | 1517.203 |
| 2.5385528 | 38.1774603 | 1590.946 |
| 96.4358540 | 42.6515159 | 1493.466 |
| 37.5199013 | 90.4118859 | 1508.448 |
| 55.8494664 | 11.3617214 | 1526.087 |
| 90.0510739 | 108.5290342 | 1509.947 |
| 81.5646677 | 27.8533115 | 1563.971 |
| 63.2692140 | 104.8468013 | 1512.176 |
| 21.9061902 | 29.5260864 | 1624.406 |
| 48.7354555 | 26.8942822 | 1549.051 |
| 65.8188237 | 83.8292525 | 1564.853 |
| 56.4777516 | 30.5094723 | 1561.733 |
| 66.6976356 | 85.4519467 | 1564.198 |
| 14.6818828 | 16.8287699 | 1555.984 |
| 68.9388928 | 51.3027483 | 1526.373 |
| 42.1837569 | 99.9531719 | 1490.173 |
| 59.7795871 | 83.3325606 | 1543.893 |
| 38.7832108 | 29.9532540 | 1557.967 |
| 1.6569089 | 48.9196271 | 1564.926 |
| 9.1434492 | 34.0696026 | 1627.608 |
| 2.9571689 | 90.6819687 | 1508.770 |
| 66.5165026 | 109.5292699 | 1511.177 |
| 70.1532239 | 42.7150412 | 1524.134 |
| -8.3372651 | 14.7380619 | 1514.970 |
| 75.3393121 | 68.1212666 | 1657.606 |
| 69.4156640 | -8.4492454 | 1499.460 |
| 27.7118924 | -6.8979296 | 1515.345 |
| 56.2264132 | 35.3259819 | 1519.458 |
| 100.0337624 | 1.7270084 | 1499.681 |
| -4.1995380 | 40.3426506 | 1580.441 |
| 2.6269938 | 26.4221574 | 1572.903 |
| 60.1177011 | 49.1849265 | 1514.151 |
| 49.5006055 | 15.8215728 | 1533.524 |
| 61.9652308 | 23.4285500 | 1597.530 |
| 34.6930682 | -1.9026141 | 1490.962 |
| 87.3495553 | 41.9261293 | 1514.999 |
| 17.6895303 | 82.4151665 | 1497.746 |
| 55.9984183 | 49.8219867 | 1493.417 |
| 58.3716456 | 105.7419121 | 1507.352 |
| 86.4720964 | 48.1162036 | 1523.178 |
| 77.7358201 | 17.3299498 | 1581.757 |
| 75.3203586 | 74.9087513 | 1721.552 |
| 82.3082523 | 42.0287462 | 1514.646 |
| 42.4700805 | 20.1111487 | 1541.480 |
| 73.4870868 | 1.7718971 | 1501.789 |
| 44.7152671 | 101.9405057 | 1506.463 |
| 0.6859109 | 27.2353932 | 1597.665 |
| 12.5278158 | 79.6288005 | 1525.128 |
| 31.5353685 | 24.1613830 | 1574.484 |
| 73.1009148 | 104.1031841 | 1528.343 |
| 7.6972773 | 54.0614536 | 1528.846 |
| 1.9605939 | 105.2482462 | 1511.436 |
| 84.4848087 | -8.5314229 | 1494.283 |
| 58.5529146 | 107.9059088 | 1501.270 |
| 3.9025865 | -3.1788485 | 1505.593 |
| 51.2491019 | 49.5866430 | 1501.208 |
| 54.3456158 | 18.5347186 | 1555.640 |
| 107.6847397 | 97.5226675 | 1500.701 |
| 24.2371608 | 45.0674648 | 1564.559 |
| 2.0416838 | 23.1091635 | 1552.594 |
| 58.9533605 | 17.0159400 | 1572.365 |
| 37.0426840 | 81.0124443 | 1520.064 |
| 41.0724442 | 93.7449142 | 1509.576 |
| -6.3120728 | 2.5762643 | 1510.855 |
| 30.0663304 | 42.0985865 | 1545.147 |
| 60.5333997 | 97.5168953 | 1523.771 |
| -5.7777314 | 17.7990940 | 1512.944 |
| 58.3574528 | 26.6432727 | 1591.132 |
| 99.5527244 | 14.1997530 | 1510.882 |
| 5.4314999 | 59.6360724 | 1502.675 |
| 27.2049331 | 58.1436156 | 1521.830 |
| 70.7133825 | 72.2600724 | 1628.442 |
| 30.2074410 | 15.7865605 | 1561.034 |
| 103.4756336 | 73.3337831 | 1497.947 |
| 52.2111962 | 97.3944622 | 1513.381 |
| 53.4971582 | 62.7045896 | 1522.696 |
| 107.2800456 | -7.6093157 | 1513.462 |
| 36.0828773 | 1.0648308 | 1516.719 |
| 9.8959398 | 38.0345417 | 1631.626 |
| 28.1724798 | 17.5128259 | 1569.083 |
| 11.4562865 | -1.4290530 | 1522.087 |
| 19.9634146 | 67.6977307 | 1510.486 |
| 77.8391698 | 84.2881631 | 1695.827 |
| 56.8943508 | 74.8307769 | 1534.896 |
| 13.6876207 | 19.0344309 | 1568.703 |
| 5.6377469 | 1.5353368 | 1530.449 |
| 98.3904664 | 47.9330029 | 1507.807 |
| 21.1893364 | 35.5331583 | 1619.810 |
| 65.8549123 | 35.4084335 | 1536.643 |
| 26.9715169 | 16.1515069 | 1558.015 |
| 68.4559015 | 62.8591251 | 1573.437 |
| 29.6779532 | 97.7635954 | 1499.683 |
| 38.5344473 | 62.5781893 | 1497.989 |
| 80.7834582 | 88.5448509 | 1662.155 |
| 93.9590048 | 104.6989836 | 1510.376 |
| 18.7237450 | 101.6639340 | 1515.993 |
| 66.4147318 | 79.8895007 | 1573.536 |
| 39.8472847 | 91.3630961 | 1490.340 |
| 47.8162886 | 58.5914050 | 1498.235 |
| 44.7660164 | 28.6819585 | 1544.230 |
| 97.9984412 | 40.8660252 | 1521.193 |
| 82.3198781 | 64.4363866 | 1644.222 |
| 54.7836843 | 84.7816309 | 1530.893 |
| 6.7822165 | 78.8207196 | 1513.144 |
| 50.7918310 | 10.1679152 | 1524.965 |
| 93.3292193 | 33.5288333 | 1529.130 |
| 27.8532724 | 67.4226551 | 1512.085 |
| 23.2761412 | 14.0967528 | 1548.178 |
| 109.4592876 | 84.4723061 | 1511.526 |
| 94.0956487 | 83.2652183 | 1533.264 |
| 81.4859804 | 25.0429054 | 1587.515 |
| -4.2042275 | 51.4680932 | 1527.781 |
| 1.5873534 | 6.6489744 | 1513.738 |
| 81.1939345 | 27.5667786 | 1583.733 |
| 99.8942479 | 18.5305668 | 1505.506 |
| 36.0007531 | 13.2020122 | 1543.004 |
| 12.9720460 | 71.6244480 | 1503.957 |
| 106.9683536 | 25.0433458 | 1505.885 |
| 2.9125746 | -3.1916066 | 1515.022 |
| 4.3695793 | 1.3057712 | 1520.694 |
| 4.1776137 | -2.0047883 | 1507.063 |
| 55.6018076 | 93.3213405 | 1525.231 |
| 38.1397455 | 85.4323965 | 1507.045 |
| 107.1729795 | 3.4698220 | 1515.953 |
| -8.0515425 | 32.4983143 | 1561.571 |
| 54.3171866 | 85.9232182 | 1527.439 |
| 39.9986217 | 88.7987212 | 1492.890 |
| 8.1586290 | 100.2437655 | 1483.974 |
| 30.5534911 | 106.8646714 | 1492.227 |
| 11.6125095 | 0.5649165 | 1499.356 |
| 17.5639922 | 12.4533953 | 1551.861 |
| 74.0025834 | 19.4443144 | 1587.757 |
| 63.8565898 | 50.4976860 | 1507.068 |
| 62.1451755 | 73.6397534 | 1566.301 |
| 108.5706315 | 39.9792884 | 1504.908 |
| -8.6882303 | 35.6378750 | 1575.384 |
| 85.7286828 | 25.2717346 | 1575.418 |
| 51.2058432 | 88.0449695 | 1535.663 |
| 12.0071275 | 4.6950964 | 1522.949 |
| 65.0518122 | 90.6750049 | 1527.383 |
| 4.7636615 | 28.2655413 | 1591.190 |
| 97.1729092 | 52.6163154 | 1519.422 |
| 46.5810002 | 22.0714145 | 1554.407 |
| -7.8478602 | 66.6177464 | 1493.155 |
| 49.6658752 | 51.2813031 | 1516.148 |
| 9.6442217 | 61.1968264 | 1516.189 |
| 0.2572076 | 70.8486255 | 1498.373 |
| 34.3806700 | 101.7533175 | 1500.140 |
| 102.8069192 | 60.0474518 | 1514.141 |
| 5.0835445 | 45.1251017 | 1581.161 |
| 104.9616358 | 86.3437259 | 1501.398 |
| 29.2241076 | 70.3265714 | 1494.560 |
| 51.9916461 | 89.2788112 | 1520.312 |
| 98.6504682 | 4.9058045 | 1507.826 |
| 87.9540403 | 90.5393373 | 1566.810 |
| 35.6069909 | 89.1899773 | 1510.483 |
| 27.3490404 | 72.2583796 | 1496.364 |
| 44.4005316 | 54.6944445 | 1505.193 |
| 78.0673646 | 108.2493421 | 1521.683 |
| 86.0345060 | 52.2455135 | 1521.489 |
| 50.7727878 | 78.7952264 | 1517.351 |
| -6.3879282 | 28.1116613 | 1552.978 |
| 16.4694860 | 45.3933454 | 1607.031 |
| 14.2883132 | 36.7943439 | 1639.067 |
| 40.4163508 | -6.8347082 | 1493.901 |
| 59.4043694 | 107.4301601 | 1514.780 |
| 26.9125733 | 4.4445444 | 1505.026 |
| 93.9689351 | 4.9250892 | 1518.441 |
| 12.8351709 | 52.2356144 | 1540.502 |
| -4.5957246 | 68.7504948 | 1504.623 |
| 90.9387196 | 38.9355121 | 1510.016 |
| -1.6660295 | 69.8818550 | 1490.004 |
| 100.7527438 | 102.5080020 | 1495.332 |
| 105.9911643 | 37.0228226 | 1494.549 |
| 79.7167846 | 94.4482587 | 1574.811 |
| 25.5021479 | 95.8724883 | 1489.824 |
| 64.8679455 | 65.9269781 | 1552.497 |
| 34.6925810 | 107.7239428 | 1508.946 |
| 97.4436983 | 53.5945425 | 1499.629 |
| 34.9634855 | 14.2292793 | 1545.566 |
| 50.6543713 | 42.7736884 | 1528.498 |
| 78.6840622 | 54.9441999 | 1561.713 |
| 103.4864686 | 0.6000343 | 1493.714 |
| 67.8934639 | 18.1810981 | 1585.116 |
| 41.2633752 | 65.9674251 | 1499.151 |
| 54.3623894 | -3.3101534 | 1503.777 |
| 85.5901665 | -8.2802795 | 1508.463 |
| 105.2418536 | 59.6567587 | 1509.805 |
| 100.6187182 | -1.9953644 | 1500.282 |
| 56.6909134 | 31.9379110 | 1553.912 |
| 66.4071529 | 58.9636526 | 1549.707 |
| 101.4136140 | 42.7639144 | 1497.432 |
| 15.8244842 | 86.5912274 | 1504.906 |
| 25.3529834 | 23.0743360 | 1587.073 |
| 6.5443245 | 76.7723261 | 1516.060 |
| 21.3540380 | 100.5852605 | 1502.077 |
| 25.2019942 | 14.2011983 | 1580.727 |
| 69.1657058 | 93.3994231 | 1547.872 |
| -0.6532370 | 49.1092228 | 1532.503 |
| 5.7016030 | 108.0617277 | 1493.535 |
| 50.6036347 | 8.3754791 | 1516.558 |
| 36.7557311 | 17.2599553 | 1566.080 |
| 91.0657953 | 105.8925509 | 1493.553 |
| -5.3899947 | 94.7447042 | 1501.037 |
| 48.0334802 | 21.4244828 | 1549.708 |
| 50.4329633 | 103.7323134 | 1519.404 |
| -5.0505084 | 66.8113191 | 1505.004 |
| 69.7325202 | 40.9201219 | 1515.943 |
| -4.8887698 | 28.0602180 | 1557.802 |
| 46.4621114 | 61.2889474 | 1503.009 |
| 64.2938416 | 35.7760528 | 1561.403 |
| 77.1037695 | 42.3929282 | 1526.516 |
| 96.2017672 | 101.9963953 | 1519.433 |
| 0.4662857 | 89.0292998 | 1483.743 |
| 72.0708827 | 59.2056660 | 1588.294 |
| 109.5364165 | 91.5752343 | 1489.485 |
| 79.2936642 | 14.1492138 | 1560.668 |
| -5.3118066 | 34.5973874 | 1585.209 |
| 63.0517392 | 6.3511438 | 1492.916 |
| 17.3393103 | 44.7068286 | 1585.357 |
| 62.9541152 | 80.9099935 | 1552.035 |
| 87.2019699 | 63.7778058 | 1582.010 |
| 90.3052440 | 69.0934487 | 1574.500 |
| 25.2688734 | 76.1048284 | 1493.799 |
| 21.4297049 | 38.1346542 | 1623.359 |
| 3.5175382 | -1.4756003 | 1509.298 |
| 76.0590468 | 91.6870374 | 1600.662 |
| 92.9337300 | 105.7332276 | 1511.960 |
| 71.6638974 | 35.1597497 | 1553.390 |
| 10.7209173 | 28.1735718 | 1633.483 |
| 51.3673127 | 108.6976360 | 1502.855 |
| 71.1383390 | 13.3389689 | 1549.392 |
| 75.5817593 | 14.6318629 | 1541.855 |
| 30.8308143 | 36.5810263 | 1589.775 |
| 67.7793316 | 16.1500442 | 1585.947 |
| 106.0849460 | 98.0246542 | 1485.666 |
| 11.2598220 | 15.8606408 | 1565.633 |
| 23.1088913 | 76.0537079 | 1479.768 |
| 10.4403996 | -7.6325232 | 1517.407 |
| 2.7615352 | 104.7978228 | 1529.142 |
| 89.0893202 | -8.5448018 | 1500.948 |
| 98.6259430 | -7.0139099 | 1510.369 |
| 14.6402051 | 103.9925252 | 1488.091 |
| 54.7124719 | 26.8612291 | 1572.300 |
| 47.1933835 | 32.0922527 | 1525.702 |
| 85.3700834 | 103.6307418 | 1515.363 |
| 1.7721502 | 15.0269578 | 1563.127 |
| 12.5727107 | 8.6170472 | 1526.653 |
| 45.5713691 | 53.7627569 | 1508.066 |
| 44.2215382 | 67.0370509 | 1497.606 |
| 74.5813477 | 101.3129812 | 1528.128 |
| 30.0115687 | 11.3819921 | 1541.190 |
| 90.4531909 | 61.6848080 | 1554.888 |
| 20.0863997 | 46.7913048 | 1570.847 |
| 47.1320663 | 38.5187942 | 1531.695 |
| 21.5359278 | 2.8504028 | 1525.109 |
| 102.1136188 | 50.2343568 | 1501.812 |
| 48.5367774 | 12.0062433 | 1535.235 |
| -2.3180114 | 55.7169349 | 1506.025 |
| 50.6102194 | 88.1279888 | 1510.649 |
| 77.8522285 | 1.2872542 | 1504.101 |
| 48.3893854 | 104.5696926 | 1492.684 |
| 66.3933251 | 101.9618268 | 1509.666 |
| 85.9375359 | 28.1935550 | 1555.860 |
| 67.8682558 | 91.8400077 | 1543.947 |
| 102.2557618 | 50.6659800 | 1493.771 |
| 89.4257922 | 46.7036010 | 1519.883 |
| 61.1501857 | 38.6312221 | 1523.934 |
| 3.1689991 | 21.9017038 | 1571.345 |
| 54.2624584 | 75.0700301 | 1520.234 |
| -2.1329539 | 0.6990997 | 1517.794 |
| 14.5326373 | 14.6623044 | 1569.137 |
| 92.0923677 | 14.5583923 | 1540.265 |
| 50.9768518 | 11.4237147 | 1523.812 |
| 25.7224895 | 62.8373531 | 1499.869 |
| 9.0913331 | 91.3946745 | 1497.405 |
| 56.7395689 | 65.3898101 | 1526.841 |
| -6.9715546 | 108.3857356 | 1505.603 |
| 86.1255730 | 37.8573761 | 1513.962 |
| 7.6243359 | -5.6062683 | 1508.418 |
| 14.3570259 | 78.8911665 | 1492.319 |
| 78.8684186 | 86.3794514 | 1701.078 |
| 50.1203243 | 6.5330506 | 1525.042 |
| 78.2810920 | 101.1884724 | 1548.463 |
| 91.9315635 | 32.0111449 | 1514.694 |
| -1.8781652 | 45.9827446 | 1553.455 |
| 19.2320025 | 48.2957502 | 1577.783 |
| 52.5425185 | 5.6003292 | 1518.864 |
| 5.5162289 | 101.0176554 | 1511.460 |
| 68.3529132 | 96.7494079 | 1542.238 |
| 29.1369224 | 15.5334408 | 1544.346 |
| -3.6228530 | 94.7946062 | 1485.693 |
| 80.2698530 | -8.4745485 | 1514.719 |
| 63.7084272 | 92.6554637 | 1534.721 |
| 5.0221692 | 5.5771130 | 1521.852 |
| 23.1965916 | 22.3894481 | 1595.744 |
| 107.4721211 | 61.1590065 | 1511.029 |
| 91.2873364 | 85.5051300 | 1557.369 |
| 35.2431964 | 29.0336821 | 1547.582 |
| 67.6516760 | 10.0403570 | 1512.336 |
| 88.2066043 | 14.8927116 | 1537.623 |
| 69.6267306 | 84.3787444 | 1594.476 |
library(plotly)
library(sm)
longitude=seq(min(long),max(long),length=100)
latitude=seq(min(lat),max(lat),length=100)
RES=sm.regression(cbind(long,lat),alti,eval.points=cbind(longitude,latitude), eval.grid=TRUE,method='cv',display='none',nbins=0)$estimate
Est_altitude=matrix(RES,100,100)
DT=data.frame(longitude,latitude,Est_altitude)
plot_ly(DT,x=~latitude,y=~longitude,z = ~Est_altitude, type = "surface")#%>% layout(xaxis = list(autorange = "reversed"))Estimated.Altitude=sm.regression(cbind(long,lat),alti,eval.points=Position, eval.grid=FALSE,method='cv',display='none',nbins=0)$estimate
lo=Position[,1]
la=Position[,2]
DT2=data.frame(Time,lo,la,Estimated.Altitude)
plot_ly(DT,x=~latitude,y=~longitude,z = ~Est_altitude, type = "surface")%>%
add_trace(x=~la,y=~lo,z = ~Estimated.Altitude, type = "scatter3d",mode='markers')plot_ly(DT2,x=~Time,y=~Estimated.Altitude, type = 'scatter', mode = 'lines')#plot_ly(DT2,x=~la,y=~lo,z = ~RES, type = "scatter3d",trace='add')
#%>% layout(xaxis = list(autorange = "reversed"))3.2.2 Discrimination
load('Jussac_EN.RData')
library(knitr)
library(kableExtra)
DF=data.frame(Type,LCB,LMS,LBM,LP,LM,LAM)
kable(DF,'html',caption="Data given in Jussac.RData") %>%
kable_styling() %>%
scroll_box(width = "100%", height='7cm')| Type | LCB | LMS | LBM | LP | LM | LAM |
|---|---|---|---|---|---|---|
| Dog | 129 | 64 | 95 | 17.5 | 11.2 | 13.8 |
| Dog | 154 | 74 | 76 | 20.0 | 14.2 | 16.5 |
| Dog | 170 | 87 | 71 | 17.9 | 12.3 | 15.9 |
| Dog | 188 | 94 | 73 | 19.5 | 13.3 | 14.8 |
| Dog | 161 | 81 | 55 | 17.1 | 12.1 | 13.0 |
| Dog | 164 | 90 | 58 | 17.5 | 12.7 | 14.7 |
| Dog | 203 | 109 | 65 | 20.7 | 14.0 | 16.8 |
| Dog | 178 | 97 | 57 | 17.3 | 12.8 | 14.3 |
| Dog | 212 | 114 | 65 | 20.5 | 14.3 | 15.5 |
| Dog | 221 | 123 | 62 | 21.2 | 15.2 | 17.0 |
| Dog | 183 | 97 | 52 | 19.3 | 12.9 | 13.5 |
| Dog | 212 | 112 | 65 | 19.7 | 14.2 | 16.0 |
| Dog | 220 | 117 | 70 | 19.8 | 14.3 | 15.6 |
| Dog | 216 | 113 | 72 | 20.5 | 14.4 | 17.7 |
| Dog | 216 | 112 | 75 | 19.6 | 14.0 | 16.4 |
| Dog | 205 | 110 | 68 | 20.8 | 14.1 | 16.4 |
| Dog | 228 | 122 | 78 | 22.5 | 14.2 | 17.8 |
| Dog | 218 | 112 | 65 | 20.3 | 13.9 | 17.0 |
| Dog | 190 | 93 | 78 | 19.7 | 13.2 | 14.0 |
| Dog | 212 | 111 | 73 | 20.5 | 13.7 | 16.6 |
| Dog | 201 | 105 | 70 | 19.8 | 14.3 | 15.9 |
| Dog | 196 | 106 | 67 | 18.5 | 12.6 | 14.2 |
| Dog | 158 | 71 | 71 | 16.7 | 12.5 | 13.3 |
| Dog | 255 | 126 | 86 | 21.4 | 15.0 | 18.0 |
| Dog | 234 | 113 | 83 | 21.3 | 14.8 | 17.0 |
| Dog | 205 | 105 | 70 | 19.0 | 12.4 | 14.9 |
| Dog | 186 | 97 | 62 | 19.0 | 13.2 | 14.2 |
| Dog | 241 | 119 | 87 | 21.0 | 14.7 | 18.3 |
| Dog | 220 | 111 | 88 | 22.5 | 15.4 | 18.0 |
| Dog | 242 | 120 | 85 | 19.9 | 15.3 | 17.6 |
| Wolf | 199 | 105 | 73 | 23.4 | 15.0 | 19.1 |
| Wolf | 227 | 117 | 77 | 25.0 | 15.3 | 18.6 |
| Wolf | 228 | 122 | 82 | 24.7 | 15.0 | 18.5 |
| Wolf | 232 | 123 | 83 | 25.3 | 16.8 | 15.5 |
| Wolf | 231 | 121 | 78 | 23.5 | 16.5 | 19.6 |
| Wolf | 215 | 118 | 74 | 25.7 | 15.7 | 19.0 |
| Wolf | 184 | 100 | 69 | 23.3 | 15.8 | 19.7 |
| Wolf | 175 | 94 | 73 | 22.2 | 14.8 | 17.0 |
| Wolf | 239 | 124 | 77 | 25.0 | 16.8 | 27.0 |
| Wolf | 203 | 109 | 70 | 23.3 | 15.0 | 18.7 |
| Wolf | 226 | 118 | 72 | 26.0 | 16.0 | 19.4 |
| Wolf | 226 | 119 | 77 | 26.5 | 16.8 | 19.3 |
Quant=cbind(LCB,LMS,LBM,LP,LM,LAM)choice=combn(1:6,2)
NC=ncol(choice)
RES=list()
for (j in 1:NC)
{
x=Quant[,choice[,j]]
ju=Jussac[choice[,j]]
cat(colnames(x),'\n')
RE=Classif_NP(x,Type,ju)
cat(paste('Error:',RE$Err,' Auc:',RE$Auc,'\n'))
cat('Confusion Matrix\n')
print(RE$M_table)
cat('Prob0\n')
print(RE$Prob0)
cat('\n \n \n')
RES[[j]]=RE
}LCB LMS
Error: 0.30952380952381 Auc: 0.652777777777778
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 11 1
Prob0
Dog Wolf
[1,] 0.9431089 0.05689111
LCB LBM
Error: 0.261904761904762 Auc: 0.747222222222222
Confusion Matrix
Class
Y Dog Wolf
Dog 27 3
Wolf 8 4
Prob0
Dog Wolf
[1,] 0.6338983 0.3661017
LCB LP
Error: 0.0476190476190477 Auc: 0.988888888888889
Confusion Matrix
Class
Y Dog Wolf
Dog 29 1
Wolf 1 11
Prob0
Dog Wolf
[1,] 0.9630371 0.0369629
LCB LM
Error: 0.166666666666667 Auc: 0.933333333333333
Confusion Matrix
Class
Y Dog Wolf
Dog 27 3
Wolf 4 8
Prob0
Dog Wolf
[1,] 1.002257 -0.002257412
LCB LAM
Error: 0.0952380952380952 Auc: 0.847222222222222
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 2 10
Prob0
Dog Wolf
[1,] 0.9635584 0.03644162
LMS LBM
Error: 0.261904761904762 Auc: 0.802777777777778
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 9 3
Prob0
Dog Wolf
[1,] 0.6363852 0.3636148
LMS LP
Error: 0.0714285714285714 Auc: 0.994444444444444
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 1 11
Prob0
Dog Wolf
[1,] 0.8624162 0.1375838
LMS LM
Error: 0.0952380952380952 Auc: 0.941666666666667
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 2 10
Prob0
Dog Wolf
[1,] 0.7724708 0.2275292
LMS LAM
Error: 0.142857142857143 Auc: 0.891666666666667
Confusion Matrix
Class
Y Dog Wolf
Dog 27 3
Wolf 3 9
Prob0
Dog Wolf
[1,] 0.6542453 0.3457547
LBM LP
Error: 0.0714285714285714 Auc: 0.977777777777778
Confusion Matrix
Class
Y Dog Wolf
Dog 29 1
Wolf 2 10
Prob0
Dog Wolf
[1,] 0.8792215 0.1207785
LBM LM
Error: 0.0476190476190477 Auc: 0.972222222222222
Confusion Matrix
Class
Y Dog Wolf
Dog 30 0
Wolf 2 10
Prob0
Dog Wolf
[1,] 0.9682595 0.03174052
LBM LAM
Error: 0.119047619047619 Auc: 0.888888888888889
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 3 9
Prob0
Dog Wolf
[1,] 0.755392 0.244608
LP LM
Error: 0.0714285714285714 Auc: 0.969444444444444
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 1 11
Prob0
Dog Wolf
[1,] 0.9381193 0.06188067
LP LAM
Error: 0.0714285714285714 Auc: 0.961111111111111
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 1 11
Prob0
Dog Wolf
[1,] 0.9514529 0.04854711
LM LAM
Error: 0.0714285714285714 Auc: 0.919444444444444
Confusion Matrix
Class
Y Dog Wolf
Dog 28 2
Wolf 1 11
Prob0
Dog Wolf
[1,] 0.9988722 0.001127843
### LBM and LM
x=Quant[,choice[,11]]
ju=Jussac[choice[,11]]
Classif_NP(x,Type,ju,plot=T)
$Class
[1] Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog
[16] Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog Dog
[31] Wolf Wolf Dog Dog Wolf Wolf Wolf Wolf Wolf Wolf Wolf Wolf
Levels: Dog Wolf
$Prob
Dog Wolf
[1,] 1.006267278 -6.267278e-03
[2,] 0.837977693 1.620223e-01
[3,] 1.000017102 -1.710184e-05
[4,] 0.997925637 2.074363e-03
[5,] 1.000000003 -3.104118e-09
[6,] 0.999999970 3.018117e-08
[7,] 1.006374765 -6.374765e-03
[8,] 0.999999883 1.169261e-07
[9,] 0.941696974 5.830303e-02
[10,] 0.514238557 4.857614e-01
[11,] 1.000000012 -1.152800e-08
[12,] 0.975014745 2.498525e-02
[13,] 0.810642337 1.893577e-01
[14,] 0.676109358 3.238906e-01
[15,] 0.944273438 5.572656e-02
[16,] 0.976878679 2.312132e-02
[17,] 0.883728678 1.162713e-01
[18,] 1.009338134 -9.338134e-03
[19,] 1.035786657 -3.578666e-02
[20,] 0.996829363 3.170637e-03
[21,] 0.810642337 1.893577e-01
[22,] 1.000027833 -2.783275e-05
[23,] 0.999999293 7.065808e-07
[24,] 0.933654853 6.634515e-02
[25,] 0.761215898 2.387841e-01
[26,] 1.000010233 -1.023271e-05
[27,] 0.999925709 7.429136e-05
[28,] 1.440744196 -4.407442e-01
[29,] 1.351458374 -3.514584e-01
[30,] 0.708027304 2.919727e-01
[31,] 0.269828188 7.301718e-01
[32,] 0.309883664 6.901163e-01
[33,] 0.771335724 2.286643e-01
[34,] 0.816009272 1.839907e-01
[35,] 0.002536301 9.974637e-01
[36,] 0.029838852 9.701611e-01
[37,] 0.235212890 7.647871e-01
[38,] 0.426114402 5.738856e-01
[39,] -0.002457640 1.002458e+00
[40,] 0.427975936 5.720241e-01
[41,] 0.001168787 9.988312e-01
[42,] -0.002457640 1.002458e+00
$M_table
Class
Y Dog Wolf
Dog 30 0
Wolf 2 10
$Err
[1] 0.04761905
$Class0
[1] Dog
Levels: Dog Wolf
$Prob0
Dog Wolf
[1,] 0.9682595 0.03174052
$Auc
[1] 0.9722222Comparison with Linear Discriminant Analysis
Type.lda <- lda(Type~., data=DF)
LD1 <- Type.lda$scaling[,1] # vecteur 1
CV=lda(Type~., data=DF,CV=TRUE)
### Estimated classes LOOCV
Classe_est=CV$class
### Confusion matrix
table(Type,Classe_est) Classe_est
Type Dog Wolf
Dog 29 1
Wolf 0 12### Error rate on learning sample
mean(Classe_est!=Type)[1] 0.02380952### ROC Curve
RO=ROC(Type=='Wolf',CV$posterior[,2],plot=TRUE)
### Class0
predict(Type.lda,Jussac)$class[1] Dog
Levels: Dog Wolf### Prob0
predict(Type.lda,Jussac)$posterior[1][1] 0.99992453.2 Functional case
3.2.1 Regression
fdaplot=function(X,ylim=NULL)
{
matplot(t(X),type='l',ylim=ylim)
}
library("fda")
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.2638plot(R0$Estimated.values,Y)
abline(0,1)
R1=funopare.kernel.cv(Y,CURVES, CURVES,1,10,c(0,1))
R1$Mse[1] 44.27052plot(R1$Estimated.values,Y)
abline(0,1)
R2=funopare.kernel.cv(Y,CURVES, CURVES,2,10,c(0,1))
R2$Mse[1] 5.784178plot(R2$Estimated.values,Y)
abline(0,1)
R3=funopare.kernel.cv(Y,CURVES, CURVES,3,10,c(0,1))
R3$Mse[1] 9.118275plot(R3$Estimated.values,Y)
abline(0,1)
Rpca=funopare.kernel.cv(Y,CURVES, CURVES,semimetric='pca',10)
Rpca$Mse[1] 122.7028plot(Rpca$Estimated.values,Y)
abline(0,1)
3.2.2 Discrimination
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)
}
#######
Classif_NP_fun=function(X,Y,X0=NA,plot=FALSE,semimetric,...)
{
X=as.matrix(X)
X0=as.matrix(X0)
n=length(Y)
V=sort(unique(Y))
n_V=length(V)
Prob=matrix(NA,n,n_V)
colnames(Prob)=V
Class=rep(NA,n)
test_sample=FALSE
if (!is.na(max(X0)) & (any(dim(X)!=dim(X0))))
{
test_sample=TRUE
}
if (!is.na(max(X0)) & (all(dim(X)==dim(X0))))
{
if (any(X!=X0)){test_sample=TRUE}
}
if (test_sample)
{
P0=matrix(NA,nrow(X0),n_V)
Class0=rep(NA,n)
}
for (v in 1:n_V)
{
z=as.numeric(Y==V[v])
Prob[,v]=funopare.kernel.cv(z,X,X,semimetric=semimetric,...)$Estimated.values
if (test_sample) {P0[,v]=funopare.kernel.cv(z,X,X0,semimetric=semimetric,...)$Predicted.values}
}
if (n_V==2) {Roc=ROC(Y==V[2],Prob[,2],plot)}
Class=V[apply(Prob,1,which.max)]
V_est=sort(unique(Class))
if (length(V_est)==n_V){M_table=table(Y,Class)}
else {
M_table=matrix(0,n_V,n_V)
M_table0=table(Y,Class)
for (j in 1:length(V_est)) {M_table[,which(V==V_est[j])]=M_table0[,j]}
}
Err=1-(sum(diag(M_table))/sum(M_table))
if (test_sample) {Class0=V[apply(P0,1,which.max)]}
if (test_sample) {return(list(Class=Class, Prob=Prob, M_table=M_table, Err=Err, Class0=Class0,Prob0=P0,Auc=ifelse(n_V==2,Roc$AUC,NA)))}
else {return(list(Class=Class, Prob=Prob, M_table=M_table, Err=Err,Auc=ifelse(n_V==2,Roc$AUC,NA)))}
}
learn=1:1000
RE=Classif_NP_fun(CURVES[learn,],PHONEME[learn],CURVES[-learn,],semimetric = 'mplsr',q=4)
### On learning sample
RE$M_table Class
Y AA AO DCL IY SH
AA 153 47 0 0 0
AO 26 173 1 0 0
DCL 0 0 199 1 0
IY 0 0 1 199 0
SH 0 0 0 0 200RE$Err[1] 0.076### On testing sample
table(PHONEME[-learn],RE$Class0)
AA AO DCL IY SH
AA 152 48 0 0 0
AO 29 171 0 0 0
DCL 1 0 199 0 0
IY 0 0 0 200 0
SH 0 0 0 0 200mean(PHONEME[-learn]!=RE$Class0)[1] 0.078