# NOT RUN {
# ---------------------------------------------------- #
# FITTING THE MODEL
# ---------------------------------------------------- #
# Load the 'city' data and wrap as 'riemobj'
data(cities)
locations = cities$cartesian
embed2 = array(0,c(60,2))
for (i in 1:60){
embed2[i,] = sphere.xyz2geo(locations[i,])
}
# Fit the model with different numbers of clusters
k2 = moSL(locations, k=2)
k3 = moSL(locations, k=3)
k4 = moSL(locations, k=4)
# Visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3))
plot(embed2, col=k2$cluster, pch=19, main="K=2")
plot(embed2, col=k3$cluster, pch=19, main="K=3")
plot(embed2, col=k4$cluster, pch=19, main="K=4")
par(opar)
# ---------------------------------------------------- #
# USE S3 METHODS
# ---------------------------------------------------- #
# Use the same 'locations' data as new data
# (1) log-likelihood
newloglkd = round(loglkd(k3, locations), 5)
fitloglkd = round(k3$loglkd, 5)
print(paste0("Log-likelihood for K=3 fitted : ", fitloglkd))
print(paste0("Log-likelihood for K=3 predicted : ", newloglkd))
# (2) label
newlabel = label(k3, locations)
# (3) density
newdensity = density(k3, locations)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab