In [1]:
#Primero leemos los datos
dat <- read.table('clean_cat.dat', header = TRUE)
summary(dat)
In [2]:
#Guardamos los datos q no vamos a usar para claisificar y los sacamos de la tabla
names <- dat$Id
visual <- dat$visual
ra <- dat$ra
dec <- dat$dec
dat <- dat[,-c(1,2,3,14)]
In [3]:
library('astrolibR')
gal.coord <- glactc(ra,dec, j = 1,, degree = TRUE, year = 2000)
l <- gal.coord$gl
b <- gal.coord$gb
In [10]:
#Hagamos un clustering!
library('mclust') #Cargamos la libreria
mc <- Mclust(dat, G = 1:2) #Probamos con una mixtura de gaussianas con 2 componentes como maximo
In [11]:
#Hagamos un par de graficos para ver como quedo
plot(mc,what='classification')
In [12]:
#Veamos si alguna de los clusters que encontro el mclust se corresponde con las confirmadas visualmente o las no
# confirmadas visualmente.
class <- mc$classification
dat_aux <- data.frame(dat, visual, class)
dat1 <- subset(dat_aux, dat_aux$class == 1)
dat2 <- subset(dat_aux, dat_aux$class == 2)
print('Cluster 1')
table(dat1$visual)
print('Cluster 2')
table(dat2$visual)
In [15]:
#Probemos ahora con un par mas de componentes
mc <- Mclust(dat, G = 1:5) #Probamos con una mixtura de gaussianas con 5 componentes como maximo
mc$G
In [14]:
#Hagamos un par de graficos para ver como quedo
plot(mc,what='classification')
In [16]:
# Grafico del criterio de informacion bayesiano para decidir con que modelo quedarse
plot(mc,what='BIC')
Si bien el primer algoritmo de clustering encuentra 2 grupos diferentes, ninguno de estos separa las que no fueron confirmadas visualmente de las que si, lo que quiere decir que estas 2 clases son muy similares!
Cuando le permito probar hasta 5 componentes, encuentra 5 grupos, pero de vuelta ninguno de estos separa las que no fueron confirmadas visualmente de las que si. Por otro lado en el grafico de BIC se ve que la performance de los modelos no cambia mucho al aumentarle el numero de componentes, lo que significa que es medio ruidoso, que no hay una clara distinción entre estas 5 componentes y que probablemente no haya una distinción física real entre estos clusters.
In [6]:
#Primero leemos los datos
dat <- read.table('clean_cat2.dat', header = TRUE)
summary(dat)
In [7]:
#Guardamos los datos q no vamos a usar para claisificar y los sacamos de la tabla
names <- dat$Id
visual <- dat$visual
ra <- dat$ra
dec <- dat$dec
dat <- dat[,-c(1,2,3,18)]
In [12]:
Z <- vector(length = length(dat$Z)) #Creamos un vector
Z[which(is.na(dat$Z) == TRUE)] = 0
Z[which(is.na(dat$Z) == FALSE)] = 1
dat$Z <- Z #Reasignamos la magnitud Z
Z2 <- vector(length = length(dat$Z2)) #Creamos un vector
Z2[which(is.na(dat$Z2) == TRUE)] = 0
Z2[which(is.na(dat$Z2) == FALSE)] = 1
dat$Z2 <- Z2 #Reasignamos la magnitud Z2
Y <- vector(length = length(dat$Y)) #Creamos un vector
Y[which(is.na(dat$Y) == TRUE)] = 0
Y[which(is.na(dat$Y) == FALSE)] = 1
dat$Y <- Y #Reasignamos la magnitud Y
Y2 <- vector(length = length(dat$Y2)) #Creamos un vector
Y2[which(is.na(dat$Y2) == TRUE)] = 0
Y2[which(is.na(dat$Y2) == FALSE)] = 1
dat$Y2 <- Y2 #Reasignamos la magnitud Y2
In [13]:
#Veamos como quedaron los datos
summary(dat)
In [19]:
#Hagamos un clustering!
library('mclust') #Cargamos la libreria
mc2 <- Mclust(dat, G = 1:2) #Probamos con una mixtura de gaussianas con 2 componentes como maximo
mc5 <- Mclust(dat, G = 1:5) #Probamos con una mixtura de gaussianas con 5 componentes como maximo
In [20]:
#Hagamos un par de graficos para ver como quedo
mc2$G
mc5$G
Si le agrego informacion sobre las magnitudes z e y (1 si las tiene, 0 sino) me encuentra un solo cluster...
In [4]:
pos <- data.frame(l, b)
plot(pos,pch=20)
In [5]:
#Cortamos por tile por que sino el mclust va a separar los 2 tiles simplemente
tile1 <- subset(pos, pos$l > 300)
tile2 <- subset(pos, pos$l < 300)
length(tile1$l)
length(tile2$l)
In [6]:
library('mclust')
mc1 <- Mclust(tile1, G = 5:15)
mc2 <- Mclust(tile2, G = 5:15)
mc1$G
mc2$G
In [7]:
plot(mc1,what='classification')
dev.copy(pdf,'cumulo.pdf')
dev.off()
In [13]:
class <- mc1$classification
tile1 <- data.frame(tile1,class)
plot(tile1$l,tile1$b,pch=20,xlab='l',ylab='b', cex=0.1)
#for(i in 1:mc1$G){
sub <- subset(tile1, tile1$class == 5)
points(sub$l, sub$b, pch = 1, col = 1)
# }
In [25]:
dat <- data.frame(dat, l, b)
tile1 <- subset(dat, dat$l > 300)
tile2 <- subset(dat, dat$l < 300)
length(tile1$l)
length(tile2$l)
tile1 <- data.frame(tile1, class)
cum <- subset(tile1, tile1$class == 5)
plot(tile1$J2, (tile1$J2-tile1$Ks2),pch=20)
points(cum$J2, (cum$J2-cum$Ks2), pch = 20, col = 'red')
hist((cum$J2-cum$Ks2))
In [20]:
mc1$parameters$variance
In [21]:
plot(mc2,what='classification')
In [22]:
mc2$parameters$variance
In [36]:
library('cluster')
In [85]:
gap1 <- clusGap(tile1, K.max = 10, FUNcluster = mclust.gap)
In [86]:
gap1
plot(gap1)
In [83]:
gap2 <- clusGap(tile2, K.max = 10, FUNcluster = mclust.gap)
In [84]:
gap2
plot(gap2)
In [81]:
mclust.gap <- function(data, k){
class <- vector(length = nrow(data))
class <- as.vector(Mclust(data, G = k)$classification)
class <- list('cluster' = class)
return(class)
}