Nothing
###########################################################################
# Copyright 2009 Michael Schweinberger #
# #
# This file is part of hergm. #
# #
# hergm is free software: you can redistribute it and/or modify #
# it under the terms of the GNU General Public License as published by #
# the Free Software Foundation, either version 3 of the License, or #
# (at your option) any later version. #
# #
# hergm is distributed in the hope that it will be useful, #
# but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
# along with hergm. If not, see <http://www.gnu.org/licenses/>. #
# #
###########################################################################
hergm.loss <- function(n_nodes, z, p)
# Loss function (Stephens, 2000, discussion): minus log joint classification probability
# input: number of nodes, categories of nodes, classification probabilities
# output: loss
{
loss <- 0
for (i in 1:n_nodes) # Node
{
k <- z[i] # Category of node
loss <- loss - log(p[i,k]) # Loss associated with classification of node to category
}
loss
}
hergm.permute <- function(n_nodes, z, nu)
# Permute labels of category indicators z in accordance with permutation nu
# input: number of nodes, categories of nodes, permutation
# output: permuted categories of nodes
{
nu_z <- vector(mode = "numeric", length = n_nodes)
for (i in 1:n_nodes)
{
k <- z[i] # Present value
nu_z[i] <- nu[k] # Permuted value
}
nu_z
}
hergm.permute_indicator <- function(n_nodes, n_categories, n_sample, indicator, permutations)
{
nu_indicator <- matrix(data = 0, nrow = n_sample, ncol = n_nodes) # Classification probabilities
z <- vector(mode = "numeric", length = n_nodes)
nu <- vector(mode = "numeric", length = n_categories)
nu_z <- vector(mode = "numeric", length = n_nodes)
for (i in 1:n_sample)
{
z <- indicator[i,]
nu <- permutations[i,] # Set nu to minimizing permutation
nu_indicator[i,] <- hergm.permute(n_nodes,z,nu) # Permute category indicators z in accordance with permutation nu
}
nu_indicator
}
hergm.sa <- function(n_nodes, n_categories, n_permutations, permutations, z, p, last_minimizer)
# Simulated Annealing
# input: number of nodes, number categories, number of permutations, MCMC sample point of categories of nodes, classification probabilities, last minimizer
# output: sample minimum, minimizer
{
nu <- permutations[last_minimizer,] # Set nu to i-th permutation
nu_z <- hergm.permute(n_nodes,z,nu) # Permute z in accordance with nu
last_loss <- hergm.loss(n_nodes,nu_z,p) # Compute loss under nu_z and p
last <- last_minimizer # Permutation
minimum <- last_loss # Minimum of loss function
minimizer <- last_minimizer # Minimizer of loss function
t <- 2 # Temperature
factor <- 0.9 # Factor by which temperature is decreased
n_iterations <- n_categories # Number of iterations
while (t > 0.1) # Number of batches with constant temperature
{
t <- factor * t # Decrease temperature
for (iteration in 1:n_iterations) # Number of iterations per batch with constant temperature
{
# Candidate: walk on "ring"
u <- runif(1, min = -10, max = 10) # Uniform(-10,10)
i <- last + round(u)
if (i < 1) i <- i + n_permutations
else if (i > n_permutations) i <- i - n_permutations
nu <- permutations[i,] # Set nu to i-th permutation
nu_z <- hergm.permute(n_nodes,z,nu) # Permute z in accordance with nu
loss <- hergm.loss(n_nodes,nu_z,p) # Compute loss under nu_z and p
# To accept or not
if (is.infinite(last_loss)) last_loss <- 1e5
if (is.infinite(loss)) loss <- 1e5
d <- loss - last_loss
if (d < 0) accept <- 1
else
{
pr <- exp(- d / t)
if (u < pr) accept <- 1
else accept <- 0
}
# If accepted, store
if (accept == 1)
{
if (loss < minimum) # Store sample minimum and minimizer
{
minimum <- loss
minimizer <- i
}
last_loss <- loss # Store loss
last <- i # Store permutation
}
}
}
SA <- list()
SA$minimum <- minimum
SA$minimizer <- minimizer
SA
}
hergm.step_1 <- function(n_nodes, n_categories, n_sample, nu_indicator)
# Step 1 of minimization algorithm
# input: number of nodes, number of categories, MCMC sample size, MCMC sample of categories of nodes
# output: classification probabilities which minimize loss function
{
p <- matrix(data = 0, nrow = n_nodes, ncol = n_categories) # Classification probabilities
for (i in 1:n_nodes)
{
for (k in 1:n_categories)
{
p[i,k] <- sum(nu_indicator[,i] == k) / n_sample # Sample proportion of times node i is member of category k
}
}
p
}
hergm.step_2 <- function(n_sample, n_nodes, n_categories, n_permutations, permutations, indicator, p, last_minimizers)
# Step 2 of minimization algorithm
# input: MCMC sample size, number of nodes, number categories, number of permutations, permutations, MCMC sample of categories of nodes, classification probabilities
# output: permutations which minimize loss function
{
min_permutations <- matrix(data = 0, nrow = n_sample, ncol = n_categories)
z <- vector(mode = "numeric", length = n_nodes)
nu <- vector(mode = "numeric", length = n_categories)
present_loss <- 0
for (h in 1:n_sample)
{
z <- indicator[h,]
if (n_categories <= 5) # If number of categories small: complete enumeration
{
# cat("Relabeling algorithm with complete enumeration")
min_loss <- Inf
for (i in 1:n_permutations)
{
nu <- permutations[i,] # Set nu to i-th permutation
nu_z <- hergm.permute(n_nodes,z,nu) # Permute z in accordance with nu
loss <- hergm.loss(n_nodes,nu_z,p) # Compute loss under nu_z and p
if (loss < min_loss)
{
min_loss <- loss
minimizer <- i
}
}
}
else # If number of categories large: simulated annealing
{
# cat("\nRelabeling algorithm with simulated annealing")
last_minimizer <- last_minimizers[h]
m <- hergm.sa(n_nodes,n_categories,n_permutations,permutations,z,p,last_minimizer)
min_loss <- m$minimum
minimizer <- m$minimizer
last_minimizers[h] <- minimizer
}
present_loss <- present_loss + min_loss # Value of loss function at the end of step 2
min_permutations[h,] <- permutations[minimizer,] # Store minimizing permutation for h-th MCMC sample point
}
step_2 <- list()
step_2$present_loss <- present_loss
step_2$min_permutations <- as.matrix(min_permutations)
step_2
}
hergm.min_loss_1 <- function(n_categories, indicator, stop_criterion, verbose)
# Minimize posterior expected loss
# input: number of categories, convergence criterion, stop criterion
# output: minimizing permutations, classification probabilities, relabeled MCMC output of categories of nodes
{
if (verbose >= 0)
{
cat("\nRelabeling MCMC sample")
cat("\n----------------------")
}
n_permutations <- factorial(n_categories)
permutations <- hergm.permutation.wrapper(n_categories) # Generate possible permutations of the category labels
if (min(indicator) == 0) # Category indicators: if category labels are 0..n_categories-1, then translate by 1 to obtain category labels 1..n_categories
{
indicator <- indicator + 1
# cat("\nInputted category indicators translated by 1")
}
n_sample <- nrow(indicator) # MCMC sample size
n_nodes <- ncol(indicator) # Number of nodes
p <- matrix(data = 0, nrow = n_nodes, ncol = n_categories) # Classification probabilities
min_permutations <- matrix(data = 0, nrow = n_sample, ncol = n_categories)
for (i in 1:n_sample)
{
k <- round(runif(n = 1, min = 1, max = n_permutations)) # Sample permutation
min_permutations[i,] <- permutations[k,]
}
min_min_loss <- Inf
min_min_permutations <- min_permutations
last_minimizers <- vector(mode = "numeric", length = n_sample) # Minimizers
last_minimizers[] <- 1
last_loss <- Inf
change <- Inf
i <- 0
while (i < stop_criterion) # Stop criterion
{
i <- i + 1
if (verbose >= 0) cat("\n\nIteration ",i)
nu_indicator <- hergm.permute_indicator(n_nodes,n_categories,n_sample,indicator,min_permutations)
p <- hergm.step_1(n_nodes,n_categories,n_sample,nu_indicator) # Minimize loss with respect to classification probabilites
step_2 <- hergm.step_2(n_sample,n_nodes,n_categories,n_permutations,permutations,indicator,p,last_minimizers) # Minimize loss with respect to permutations of category labels
min_permutations <- step_2$min_permutations
present_loss <- step_2$present_loss
if (present_loss < min_min_loss)
{
min_min_loss <- present_loss
min_min_permutations <- min_permutations
}
if (verbose >= 0) cat("\nLoss:",present_loss)
change <- present_loss - last_loss
if (verbose >= 0) cat("\nChange",change)
last_loss <- present_loss
if (abs(change) <= 1) i <- stop_criterion
}
min_permutations <- min_min_permutations
indicator <- hergm.permute_indicator(n_nodes,n_categories,n_sample,indicator,min_permutations)
if (verbose >= 0) cat("\n\nMCMC sample relabeled\n")
minimizer <- list()
minimizer$loss <- last_loss
minimizer$min_permutations <- min_permutations
minimizer$p <- p
minimizer$indicator <- indicator
minimizer
}
hergm.min_loss_2 <- function(n_categories, indicator, verbose)
# Minimize posterior expected loss
# input: number of categories, convergence criterion, stop criterion
# output: minimizing permutations, classification probabilities, relabeled MCMC output of categories of nodes
{
if (verbose >= 0)
{
cat("\nRelabeling MCMC sample...")
cat("\n-------------------------")
}
n_permutations <- factorial(n_categories)
permutations <- hergm.permutation.wrapper(n_categories) # Generate possible permutations of the category labels
if (min(indicator) == 0) # Category indicators: if category labels are 0..n_categories-1, then translate by 1 to obtain category labels 1..n_categories
{
indicator <- indicator + 1
# cat("\nInputted category indicators translated by 1")
}
n_sample <- nrow(indicator) # MCMC sample size
n_nodes <- ncol(indicator) # Number of nodes
p <- matrix(data = 0, nrow = n_nodes, ncol = n_categories) # Classification probabilities
min_permutations <- matrix(data = 0, nrow = n_sample, ncol = n_categories)
#print("Indicators:")
#print(indicator)
for (s in 1:n_sample)
{
n_unique <- 0
for (i in 1:n_nodes)
{
k <- indicator[s,i] # Label of block of node i
if (min_permutations[s,k] == 0) # Check whether block label has not shown up. If block label has not shown up, increment number of unique block labels and relabel block
{
n_unique <- n_unique + 1
min_permutations[s,k] <- n_unique
indicator[s,i] <- min_permutations[s,k]
}
else indicator[s,i] <- min_permutations[s,k] # Otherwise, relabel block by retrieving permuted block label from min_permutations
}
}
#print("Relabeled indicators:")
#print(indicator)
for (i in 1:n_nodes)
{
z <- indicator[,i] # Column i is vector of relabeled block memberships of node i
for (k in 1:n_categories)
{
p[i,k] <- sum(z==k) / n_sample # Estimated posterior classification probabilities based on relabeled MCMC sample
}
}
#print("Estimated posterior classification probabilities based on relabeled MCMC sample:")
#print(p)
if (verbose >= 0) cat("\n\nMCMC sample relabeled\n")
minimizer <- list()
minimizer$loss <- NULL
minimizer$min_permutations <- min_permutations
minimizer$p <- p
minimizer$indicator <- indicator
minimizer
}
hergm.permute_mcmc <- function(mcmc, n_categories, min_permutations)
# Relabel MCMC sample
# input: name of file containing MCMC sample, minimizing permutations
# output: MCMC sample relabeled
{
min_mcmc <- matrix(data = 0, nrow = nrow(mcmc), ncol = n_categories)
for (h in 1:nrow(mcmc))
{
for (i in 1:n_categories)
{
k <- min_permutations[h,i] # Present permutation: identity permutation
min_mcmc[h,k] <- mcmc[h,i] # Minimizing permutation
}
}
min_mcmc
}
hergm.pie <- function(filename, n_nodes, p)
# Make one pie chart for each node, representing its classification probabilities
# input: name of file to which output ought to be written, number of nodes, classification probabilities
# output: pie charts
{
for (i in 1:n_nodes)
{
file_i <- paste(filename, "_", i, ".pdf", sep = "")
pdf(file_i)
pie(p[i,], radius = 1, clockwise = TRUE, col = c("white", "orange", "brown", "red", "black"), main = "", labels = c(1:100), cex = 3.5)
dev.off()
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.