Guess, Poisson, Max-Log-Likelihood

This commit is contained in:
Noah L. Schrick 2022-04-11 16:21:34 -05:00
parent f02246966c
commit f4b9a9f7a6

View File

@ -7,22 +7,47 @@ library(igraph)
library(igraphdata) library(igraphdata)
data(yeast) data(yeast)
g <- yeast g <- yeast
g.netname <- "Yeast"
################# Set up Work #################
g.vec <- degree(g)
g.hist <- hist(g.vec, freq=FALSE, main=paste("Histogram of the", g.netname,
" Network"))
legend("topright", c("Guess", "Poisson", "Least-Squares Fit",
"Max Log-Likelihood"), lty=c(1,2,3,4), col=c("#40B0A6",
"#006CD1", "#E66100", "#D35FB7"))
g.mean <- mean(g.vec)
g.seq <- 0:max(g.vec) # x-axis
################# Guessing Alpha #################
alpha.guess <- 1.5
lines(g.seq, g.seq^(-alpha.guess), col="#40B0A6", lty=1, lwd=3)
################# Poisson #################
g.pois <- dpois(g.seq, g.mean, log=F)
lines(g.seq, g.pois, col="#006CD1", lty=2, lwd=3)
################# Linear model: Least-Squares Fit ################# ################# Linear model: Least-Squares Fit #################
g.hist <- hist(degree(g), freq=FALSE) g.breaks <- g.hist$breaks[-c(1)] # remove 0
g.seq <- 0:max(degree(g)) # x-axis
g.breaks <- g.hist$breaks[-c(1,2)] # remove 0 and low degrees
g.probs <- g.hist$density[-1] # make lengths match g.probs <- g.hist$density[-1] # make lengths match
# Need to clean up probabilities that are 0 # Need to clean up probabilities that are 0
nz.probs.mask <- g.probs!=0 nz.probs.mask <- g.probs!=0
g.breaks.clean <- g.breaks[nz.probs.mask] g.breaks.clean <- g.breaks[nz.probs.mask]
g.probs.clean <- g.breaks[nz.probs.mask] g.probs.clean <- g.probs[nz.probs.mask]
plot(log(g.breaks.clean), log(g.probs.clean)) #plot(log(g.breaks.clean), log(g.probs.clean))
g.fit <- lm(log(g.probs.clean)~log(g.breaks.clean)) g.fit <- lm(log(g.probs.clean)~log(g.breaks.clean))
summary(g.fit) summary(g.fit)
coef(g.fit)[2] alpha.LM <- coef(g.fit)[2]
lines(g.seq, g.seq^(-alpha.LM), col="#E66100", lty=3, lwd=3)
################# Max-Log-Likelihood ################# ################# Max-Log-Likelihood #################
n <- length(g.breaks.clean)
kmin <- g.breaks.clean[1]
alpha.ML <- 1 + n/sum(log(g.breaks.clean/kmin))
alpha.ML
lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4, lwd=3)