Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
Scarica il documento per vederlo tutto.
vuoi
o PayPal
tutte le volte che vuoi
Formattazione del testo con tag HTML
PBcito.lda <- lda(group ~ Infg, data=cito)cito.lda
# Supponiamo di avere un nuovo valore di Infg ("0"), calcoliamo la probabilità aposteriori e la classificazione (implicitamente fissiamo una threshold)
x <- data.frame(Infg = 0)
# Il comando predict() restituisce la classe con posterior maggiore e le probabilità aposteriori per le classi
predict(cito.lda, x)$class
predict(cito.lda, x)$posterior
# Plottiamo
x11()
plot(cito$Infg[A], rep(0, length(A)), pch=16, col='blue', ylim=c(0,1),xlab='x', ylab='estimated posterior', main="LDA", xlim=range(cito$Infg))
points(cito$Infg[B], rep(0, length(B)), pch=16, col='red')
abline(v=0, col='grey')
points(c(0,0),c(predict(cito.lda, data.frame(Infg = 0))$posterior),col=c('blue', 'red'), pch='*', cex=2.5)
# memorizzo la posterior su una griglia
x <- data.frame(Infg=seq(-10, 35, 0.5))
cito.LDA.A <- predict(cito.lda, x)$posterior[,1]
posterior per il gruppo Acito.LDA.B <- predict(cito.lda, x)$posterior[,2] # posterior per il gruppo B
Stima delle classi:predict(cito.lda, x)$class
Aggiungiamo al plot anche la posteriorlines(x[,1], cito.LDA.A, type='l', col='blue', xlab='x', ylab='estimated posterior',main="LDA")
points(x[,1], cito.LDA.B, type='l', col='red')
abline(h = 0.5)
legend(-10, 0.9, legend=c('P(A|X=x)', 'P(B|X=x)'), fill=c('blue','red'), cex = 0.7)
CAMBIARE LA PRIOR: nella popolazione, la probabilità di essere gruppo A è 5%, B95%
N.B.: i gruppi sono in ordine alfabeticocito.lda.1 <- lda(group ~ Infg, data=cito, prior=c(0.05,0.95))
cito.lda.1x <- data.frame(Infg=seq(-10, 35, 0.5))
cito.LDA.A.1 <- predict(cito.lda.1, x)$posterior[,1] # posterior per Acito.LDA.B.1 <- predict(cito.lda.1, x)$posterior[,2] # posterior per Bx11()plot (x[,1], cito.LDA.A.1, type='l', col='blue',
m <- colMeans(iris2)
m1 <- colMeans(iris2[i1,])
m2 <- colMeans(iris2[i2,])
m3 <- colMeans(iris2[i3,])
S1 <- cov(iris2[i1,])
S2 <- cov(iris2[i2,])
S3 <- cov(iris2[i3,])
Sp <- ((n1-1)*S1+(n2-1)*S2+(n3-1)*S3)/(n-g)
# LDA altro modo
da.iris <- lda(iris2, species.name)
# stima da solo prior=c(0.99,0.005, 0.005))
lda.iris
# Visualizziamo in 3D le posterior per ogni gruppo
x <- seq(min(iris[,1]), max(iris[,1]), length=50)
y <- seq(min(iris[,2]), max(iris[,2]), length=50)
xy <- expand.grid(Sepal.Length=x, Sepal.Width=y)
# plotto punti su piano xy, su z valore densità gaussiana stimata per ogni gruppo e moltiplicata per prior
open3d()
points3d(iris2[i1,1], iris2[i1,2], 0, col='red', pch=15)
points3d(iris2[i2,1], iris2[i3,2], 0, col='green', pch=15)
points3d(iris2[i3,1], iris2[i2,2], 0, col='blue',
pch=15)surface3d(x,y,matrix(dmvnorm(xy, m1, Sp) / 3, 50), alpha=0.4, color='red') surface3d(x,y,matrix(dmvnorm(xy, m2, Sp) / 3, 50), alpha=0.4, color='green', add=T) surface3d(x,y,matrix(dmvnorm(xy, m3, Sp) / 3, 50), alpha=0.4, color='blue', add=T) box3d()axes3d() # classifico ogni punto con campana più alta # Ora plottiamo la posteriorden = matrix(dmvnorm(xy, m1, Sp) / 3) + matrix(dmvnorm(xy, m2, Sp) / 3) +matrix(dmvnorm(xy, m3, Sp) / 3) open3d() points3d(iris2[i1,1], iris2[i1,2], 0, col='red', pch=15) points3d(iris2[i2,1], iris2[i3,2], 0, col='green', pch=15) points3d(iris2[i3,1], iris2[i2,2], 0, col='blue', pch=15) surface3d(x,y,matrix(dmvnorm(xy, m1, Sp) / 3/den, 50), alpha=0.4, color='red') surface3d(x,y,matrix(dmvnorm(xy, m2, Sp) / 3/den, 50), alpha=0.4, color='green',add=T) surface3d(x,y,matrix(dmvnorm(xy, m3, Sp) / 3/den, 50), alpha=0.4, color='blue',add=T) box3d()axes3d() # Proiettiamo le linee
che delimitano i valori di classificazione sul piano xy
<code>plot(iris2, main='Iris Sepal', xlab='Sepal.Length', ylab='Sepal.Width', pch=20)
points(iris2[i1,], col='red', pch=20)
points(iris2[i2,], col='green', pch=20)
points(iris2[i3,], col='blue', pch=20)
legend(min(iris[,1]), max(iris[,2]), legend=levels(species.name),fill=c('red','green','blue'), cex=.7)
# Plottiamo le medie
points(lda.iris$means, pch=4,col=c('red','green','blue') , lwd=2, cex=1.5)
# Rifacciamo la griglia
x <- seq(min(iris[,1]), max(iris[,1]), length=200)
y <- seq(min(iris[,2]), max(iris[,2]), length=200)
xy <- expand.grid(Sepal.Length=x, Sepal.Width=y)
# Calcoliamo le posterior
z <- predict(lda.iris, xy)$post # data da P_i*f_i(x,y)
# confini delle aree di previsione
z1 <- z[,1] - pmax(z[,2], z[,3]) # P_1*f_1(x,y)-max{P_j*f_j(x,y)}
z2 <- z[,2] - pmax(z[,1], z[,3]) #
P_2*f_2(x,y)-max{P_j*f_j(x,y)}z3 <- z[,3] - pmax(z[,1], z[,2])
# P_3*f_3(x,y)-max{P_j*f_j(x,y)}which(z1>0)[1]z[3801,]
# Facciamo il plot delle curve di livello w ottengo piano diviso in tre regioni con
contorni rettilinei
contour(x, y, matrix(z1, 200), levels=0, drawlabels=F, add=T)
contour(x, y, matrix(z2, 200), levels=0, drawlabels=F, add=T)
contour(x, y, matrix(z3, 200), levels=0, drawlabels=F, add=T)
# Facciamo previsione ottimista sul "training set" con matrice di confusione
Lda.iris <- predict(lda.iris, iris2)
Lda.iris$classspecies.nametable(true.class=species.name, estimated.class=Lda.iris$class)
errors <- (Lda.iris$class != species.name)
errorssum(errors)length(species.name)lda_tr_error_rate <- sum(errors)/length(species.name)lda_tr_error_rate
# QDA (già verificato prima e unica assunzione con LDA)
qda.iris <- qda(iris2, species.name)
qda.iris
# Partizioniamo (esattamente con lo stesso metodo di prima) lo spazio usando QDA
x11()
plot(iris2, main='Iris Sepal',
xlab='Sepal.Length', ylab='Sepal.Width', pch=20) points(iris2[i1,], col='red', pch=20) points(iris2[i2,], col='green', pch=20) points(iris2[i3,], col='blue', pch=20) legend(min(iris[,1]), max(iris[,2]), legend=levels(species.name),fill=c('red','green','blue')) points(qda.iris$means, col=c('red','green','blue'), pch=4, lwd=2, cex=1.5) x <- seq(min(iris[,1]), max(iris[,1]), length=200) y <- seq(min(iris[,2]), max(iris[,2]), length=200) xy <- expand.grid(Sepal.Length=x, Sepal.Width=y) z <- predict(qda.iris, xy)$post z1 <- z[,1] - pmax(z[,2], z[,3]) z2 <- z[,2] - pmax(z[,1], z[,3]) z3 <- z[,3] - pmax(z[,1], z[,2]) contour(x, y, matrix(z1, 200), levels=0, drawlabels=F, add=T) contour(x, y, matrix(z2, 200), levels=0, drawlabels=F, add=T) contour(x, y, matrix(z3, 200), levels=0, drawlabels=F, add=T)```html
m2, S2) / 3) +matrix(dmvnorm(xy, m3, S3)/3)# Guardiamo di nuovo matrice di confusione e traning errorQda.iris <- predict(qda.iris, iris2)Qda.iris$classspecies.nametable(classe.vera=species.name, classe.allocata=Qda.iris$class)erroriq <- (Qda.iris$class != species.name)erroriqqda_tr_error_rate <- sum(erroriq)/length(species.name)qda_tr_error_ratelda_tr_error_rate# K-NEAREST NEIGHBOR classifierpacman::p_load(MASS,GGally,rgl,tidyverse, ggfortify, class, boot)library(datasets)cito <- read.table('Data/cytokines.txt', header=T)x <- data.frame(Infg=seq(-10, 35, 0.5))cito.knn <- knn(train = cito$Infg, test = x, cl = cito$group, k = 3, prob=T)cito.knn# Restituisce prima la classificazione del test set e poi, se sono state richieste, leprobabilità/frequenze - da estrarre con "attr"# N.B. i gruppi "A" e "B" sono in ordine alfabetico, e diventano rispettivamente classi"0" e "1"summary(cito.knn)aux =
```attributes(cito.knn)
#solo frequenze
aux$probcito.knn.class <- (cito.knn == 'B')
# Per avere la "probabilità a posteriori" di group B devo prendere "$prob" per quelli
# classificati come "B", per gli altri "1-$prob"
cito.knn.B <- ifelse(cito.knn.class==1,attributes(cito.knn)$prob,1 - attributes(cito.knn)$prob)
x11()
points(x[,1], cito.knn.B, type='l', col='red', lty=4)
abl