pp <- numeric(13)
for (i in 1:13) {
n <- i
nn <- 0:(n-1)
x <- (12-nn)/12
pp[i] <- prod(x)
}
round(data.frame(pn=pp, pnc=1-pp), 4)
# ue2.9
f <- function(N, k) {
NN <- 0:(N-k)
x <- (-1)^NN/factorial(NN)
sum(x)/factorial(k)
}
N <- 8
round(sapply(0:8, f, N=8), 4)
factorial(8)*sapply(0:8, f, N=8)
View(f)
# ue2.8
pp <- numeric(13)
for (i in 1:13) {
n <- i
nn <- 0:(n-1)
x <- (12-nn)/12
pp[i] <- prod(x)
}
round(data.frame(pn=pp, pnc=1-pp), 4)
help
help(prod)
# ue2.8
pp <- numeric(13)
for (i in 1:13) {
n <- i
nn <- 0:(n-1)
x <- (12-nn)/12
pp[i] <- prod(x)
}
round(data.frame(pn=pp, pnc=1-pp), 4)
f <- function(N, k) {
NN <- 0:(N-k)
x <- (-1)^NN/factorial(NN)
sum(x)/factorial(k)
}
N <- 8
round(sapply(0:8, f, N=8), 4)
factorial(8)*sapply(0:8, f, N=8)
N <- 10000
DK <- numeric(N)
for (i in 1:N) {
x <- sort(runif(2))
y <- sort(c(x[1],x[2]-x[1],1-x[2]))
DK[i] <- ifelse(y[1]+y[2] > y[3], 1, 0)
}
DKum <- cumsum(DK)
plot(1:N, DKum/(1:N), type="l", lty=1, lwd=2, log="x", xlab="n", col="blue3", ylab=expression(paste(h[n],"(Dreieck konstruierbar)")))
abline(h=1/4, lty=2)
# ue2.23
p <- seq(0, 1, by=0.001)
p.intakt <- 2*p^2+2*p^3-5*p^4+2*p^5
plot(p, p.intakt, type="l", lty=1, lwd=2, asp=1, xlab="p", ylab="P( System intakt )", axes=F)
lines(c(0,1), c(0,1), lty=2)
axis(1, pos=0)
axis(2, pos=0)
axis(3, pos=1)
axis(4, pos=1)
x <- 1:6
p.max <- (2*x-1)/36
Fmax <- cumsum(p.max)
p.min <- (13-2*x)/36
Fmin <- cumsum(p.min)
plot(c(0.5,x,6.5), c(0,Fmax,1), type="s", lty=1, lwd=3, xlab="x", ylab="F(x)", col="grey50", main=expression(paste("Verteilungsfunktion für ", X[min], " und ", X[max])))
lines(c(0.5,x,6.5), c(0,Fmin,1), type="s", lty=1, lwd=3, col="grey50")
abline(h=c(0,1), lty=2)
abline(v=c(1,6), lty=2)
text(locator(1), "Max", cex=1.5)
# ue3.2
N <- 20; M <- 10
p <- M/(N+M)
x <- 0:15
Fx <- 1-(1-p)^x
plot(x, Fx, type="s", lty=1, lwd=3, xlab="x", ylab="F(x)", col="grey50", main=paste("Verteilungsfunktion (", "N = ", N, ", M = ", M, ")"))
abline(h=c(0,1), lty=2)
help(replicate)
help(sapply)
v2 <- replicate(100, erm <- sample_gnp(n=200, p=0.02))
v2 <- replicate(100, sample_gnp(n=200, p=0.02))
library("igraph")
v2 <- replicate(100, sample_gnp(n=200, p=0.02))
v2
degree(v2)
v1 <- sapply(1:100, function(n)source("my_script_file.R")$value)
sapply(v2, mean(degree(erm)))
f1 <- function(x) {
mean(degree(x))
}
sapply(v2, f1)
f1 <- function(x) {
#mean(
degree(x)
# )
}
sapply(v2, f1)
f1 <- function(x) {
#mean(
#degree(
x
#)
# )
}
sapply(v2, f1)
f1 <- function(x) {
mean(
#degree(
x
#)
)
}
sapply(v2, f1)
f1 <- function(x) {
mean(x)
#degree(
#x
#)
#)
}
sapply(v2, f1)
f1 <- function(x) {
mean(degree(x))
}
v3 <- sapply(1:100, sample_gnp(n=200, p=0.02))
v3 <- rapply(1:100, sample_gnp(n=200, p=0.02))
v3 <- lapply(1:100, sample_gnp(n=200, p=0.02))
v3 <- apply(1:100, sample_gnp(n=200, p=0.02))
x <- list()
for (i in seq_along(1:100)) {
erm <- sample_gnp(n=200, p=0.02)
list.append(x, mean(degree(erm)))
}
x <- list()
for (i in seq_along(1:100)) {
erm <- sample_gnp(n=200, p=0.02)
c(x, mean(degree(erm)))
}
x
for (i in seq_along(1:100)) {
i
erm <- sample_gnp(n=200, p=0.02)
c(x, mean(degree(erm)))
}
for (i in seq_along(1:100)) {
erm <- sample_gnp(n=200, p=0.02)
c(x, mean(degree(erm)))
}
for(i in seq_along(1:100)){
erm <- sample_gnp(n=200, p=0.02)
c(x, mean(degree(erm)))
}
sapply(v2, function(x) {
mean(degree(x))
})
x <- list()
for(i in seq_along(1:100)){
erm <- sample_gnp(n=200, p=0.02)
c(x, mean(degree(erm)))
}
for(i in 1:100){
erm <- sample_gnp(n=200, p=0.02)
c(x, mean(degree(erm)))
}
c(x, 1)
x = c(x, 1)
x <- list()
for(i in seq_along(1:100){
erm <- sample_gnp(n=200, p=0.02)
x= c(x, mean(degree(erm)))
}
#TODO: How close is this mean to the expected average degree?
########## 2.3 ##########
par(mfrow=c(1,2))
#p=0
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
########## 2.4 ##########
#p=0.01
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0.01)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
#p=0.1
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0.1)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
#p=0.6
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0.6)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
#While the average degree doesn't change we can see that the distribution and the graph do. With a higher p we get more connections to non-neighbor nodes and the degree distribution gets broader.
par(mfrow=c(1,1))
########## 2.5 ##########
#Generate network
ba <- barabasi.game(200, 1, 2, directed = F)
#Average degree
mean(degree(ba))
#Plot degree distribution
plot(degree.distribution(ba), xlab = "degree")
#The average degree is still around 4, the degree distribution has a high peak at 4 but is spread out to 27
# > prop.table(table(degree(ba)))
#
# 2     3     4     5     6     7     8     9    10    11    12    13    15    18    19    22    27
# 0.480 0.190 0.090 0.070 0.045 0.020 0.030 0.015 0.005 0.010 0.010 0.005 0.005 0.005 0.010 0.005 0.005
########## 2.6 ##########
#Generate network
ba <- barabasi.game(2000, 1, 2, directed = F)
#Average degree
mean(degree(ba))
#Plot the network
V(ba)$size <- degree(ba, v = V(ba))
plot(ba, vertex.label=NA)
#Plot the cumulative degree distribution
par(mfrow=c(1,2))
plot(degree.distribution(ba, cumulative = T), xlab = "degree k", ylab = "Normal scale")
plot(degree.distribution(ba, cumulative = T), log="xy", xlab = "degree k", ylab = "Log-scale")
for(i in seq_along(1:100){
erm <- sample_gnp(n=200, p=0.02)
x= c(x, mean(degree(erm)))
}
#TODO: How close is this mean to the expected average degree?
########## 2.3 ##########
par(mfrow=c(1,2))
#p=0
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
########## 2.4 ##########
#p=0.01
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0.01)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
#p=0.1
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0.1)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
#p=0.6
#Generate the network
ws <- watts.strogatz.game(1, 200, 2, 0.6)
#Plot the network
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
#Average degree
mean(degree(ws))
#Plot degree distribution
plot(degree.distribution(ws), xlab = "degree")
#While the average degree doesn't change we can see that the distribution and the graph do. With a higher p we get more connections to non-neighbor nodes and the degree distribution gets broader.
par(mfrow=c(1,1))
########## 2.5 ##########
#Generate network
ba <- barabasi.game(200, 1, 2, directed = F)
#Average degree
mean(degree(ba))
#Plot degree distribution
plot(degree.distribution(ba), xlab = "degree")
#The average degree is still around 4, the degree distribution has a high peak at 4 but is spread out to 27
# > prop.table(table(degree(ba)))
#
# 2     3     4     5     6     7     8     9    10    11    12    13    15    18    19    22    27
# 0.480 0.190 0.090 0.070 0.045 0.020 0.030 0.015 0.005 0.010 0.010 0.005 0.005 0.005 0.010 0.005 0.005
########## 2.6 ##########
#Generate network
ba <- barabasi.game(2000, 1, 2, directed = F)
#Average degree
mean(degree(ba))
#Plot the network
V(ba)$size <- degree(ba, v = V(ba))
plot(ba, vertex.label=NA)
#Plot the cumulative degree distribution
par(mfrow=c(1,2))
plot(degree.distribution(ba, cumulative = T), xlab = "degree k", ylab = "Normal scale")
plot(degree.distribution(ba, cumulative = T), log="xy", xlab = "degree k", ylab = "Log-scale")
for(i in seq_along(1:100)){
erm <- sample_gnp(n=200, p=0.02)
x= c(x, mean(degree(erm)))
}
library("igraph")
x <- list()
for(i in seq_along(1:100)){
erm <- sample_gnp(n=200, p=0.02)
x= c(x, mean(degree(erm)))
}
mean(x)
average(x)
x <- c(x, mean(degree(erm)))
x <- list()
for(i in seq_along(1:100)){
erm <- sample_gnp(n=200, p=0.02)
x <- c(x, mean(degree(erm)))
}
x
choose(200, 2)*0.02
choose(200, 2)
choose(200, 2)*0.02
choose(200, 2)*2
choose(200, 0.02)*0.02
choose(200, 2)*0.02
(200-1)*0.02
erm <- function() {
mean(degree(sample_gnp(n=200, p=0.02)))
}
results <- ldply(1:100, function(i) data.frame(iteration = i, means = erm() ))
mean(results$means)
library(plyr)
erm <- function() {
mean(degree(sample_gnp(n=200, p=0.02)))
}
results <- ldply(1:100, function(i) data.frame(iteration = i, means = erm() ))
mean(results$means)
sd(results$means)
library("plyr")
results <- ldply(1:100, function(i) data.frame(iteration = i, means = {
mean(degree(sample_gnp(n=200, p=0.02)))
} ))
mean(results$means)
sd(results$means)
results <- ldply(1:100, function() data.frame(iteration = i, means = {
mean(degree(sample_gnp(n=200, p=0.02)))
}))
library("igraph")
library("plyr")
par(mfrow=c(2,2))
erm <- sample_gnm(n=200, m=400)
mean(degree(erm))
prop.table(table(degree(erm)))
V(erm)$color <- degree(erm, v = V(erm))
plot(erm, vertex.size=3, vertex.label=NA)
plot(degree.distribution(erm), xlab = "degree k", ylab = "Fraction of nodes with degree k")
erm2 <- sample_gnm(n=200, m=400)
mean(degree(erm2))
prop.table(table(degree(erm2)))
V(erm2)$color <- degree(erm2, v = V(erm))
plot(erm2, vertex.size=3, vertex.label=NA)
plot(degree.distribution(erm2), xlab = "degree k", ylab = "Fraction of nodes with degree k")
par(mfrow=c(1,1))
erm <- sample_gnp(n=200, p=0.02)
mean(degree(erm))
x <- ldply(1:100, function(i) data.frame(iteration = i, means = {
mean(degree(sample_gnp(n=200, p=0.02)))
}))
mean(x$means)
sd(x$means)
(200-1)*0.02
par(mfrow=c(1,2))
ws <- watts.strogatz.game(1, 200, 2, 0)
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
mean(degree(ws))
plot(degree.distribution(ws), xlab = "degree")
ws <- watts.strogatz.game(1, 200, 2, 0.01)
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
mean(degree(ws))
plot(degree.distribution(ws), xlab = "degree")
ws <- watts.strogatz.game(1, 200, 2, 0.1)
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
mean(degree(ws))
plot(degree.distribution(ws), xlab = "degree")
ws <- watts.strogatz.game(1, 200, 2, 0.6)
V(ws)$color <- degree(ws, v = V(ws))
plot(ws, layout=layout.circle, vertex.label=NA, vertex.size=3)
mean(degree(ws))
plot(degree.distribution(ws), xlab = "degree")
par(mfrow=c(1,1))
ba <- barabasi.game(200, 1, 2, directed = F)
mean(degree(ba))
plot(degree.distribution(ba), xlab = "degree")
ba <- barabasi.game(2000, 1, 2, directed = F)
mean(degree(ba))
V(ba)$size <- degree(ba, v = V(ba))
plot(ba, vertex.label=NA)
par(mfrow=c(1,2))
plot(degree.distribution(ba, cumulative = T), xlab = "degree k", ylab = "Normal scale")
plot(degree.distribution(ba, cumulative = T), log="xy", xlab = "degree k", ylab = "Log-scale")
setwd("C:/Users/Peter/Uni/Eco")
########## 4.6 ##########
wagedata = read.table("wagedata.csv")
log_edu = lm(log(wage) ~ educ + exper, data = wagedata)
sq_log_edu = lm(log(wage) ~ educ + exper + I(exper*exper), data = wagedata)
#a)
summary(log_edu) #R-squared:  0.2493, Adjusted R-squared:  0.2465
summary(sq_log_edu) #R-squared:  0.3003, Adjusted R-squared:  0.2963
#b)
b_exp = log_edu$coefficients["exper"]
b_exp2 = sq_log_edu$coefficients["exper"]
exper = wagedata$exper
plot(exper, b_exp*exper+b_exp2*exper*exper)
sq_log_edu$coefficients
sq_log_edu$coefficients["I(exper * exper)"]
#b)
b_exp = sq_log_edu$coefficients["exper"]
b_exp2 = sq_log_edu$coefficients["I(exper * exper)"]
exper = wagedata$exper
plot(exper, b_exp*exper+b_exp2*exper*exper)
plot(b_exp*exper+b_exp2*exper*exper, exper)
help(plot)
########## 4.7 ##########
bwght = read.table("bwght.csv", header=TRUE, sep = ",")
#a)
weight = lm(bwght ~ cigs + faminc, data = bwght)
summary(weight)
#b)
bwght2 = bwght
bwght2$bwght = bwght$bwght * 453.59237
weight2 = lm(bwght ~ cigs + faminc, data = bwght2)
summary(weight2)
#c)
weight_nI = lm(bwght ~ cigs + faminc + 0, data = bwght)
summary(weight_nI) #0.7237
1 - sum(weight_nI$residuals) / sum(bwght$bwght) #0.7402849
sum(weight_nI$fitted.values) / sum(bwght$bwght) #0.7402849
#b)
ounce_to_gramm = 28.3495
bwght2 = bwght
bwght2$bwght = bwght$bwght * ounce_to_gramm
weight2 = lm(bwght ~ cigs + faminc, data = bwght2)
summary(weight2)
#R² still is 0.0298, all other betas are also (nearly) equal after conversion
3316.1581 * ounce_to_gramm #Intercept: 53058.57 converted, now 53058.57
-13.1374 * ounce_to_gramm #cigs: -210.1992 converted, now -210.20
summary(weight)
#R² still is 0.0298, all other betas are also (nearly) equal after conversion
116.97413 * ounce_to_gramm #Intercept: 53058.57 converted, now 94011.42
-0.46341 * ounce_to_gramm #cigs: -210.1992 converted, now -372.4387
0.09276 * ounce_to_gramm #faminc: 42.07523 converted, now 42.08
summary(weight)
summary(weight2)
summary(weight_nI) #0.7237
log_edu = lm(log(wage) ~ educ + exper, data = wagedata)
sq_log_edu = lm(log(wage) ~ educ + exper + I(exper*exper), data = wagedata)
#a)
summary(log_edu) #R-squared:  0.2493, Adjusted R-squared:  0.2465
summary(sq_log_edu) #R-squared:  0.3003, Adjusted R-squared:  0.2963
#b)
b_exp = sq_log_edu$coefficients["exper"]
b_exp2 = sq_log_edu$coefficients["I(exper * exper)"]
exper = wagedata$exper
plot(b_exp*exper+b_exp2*exper*exper, exper)
plot(exper, b_exp*exper+b_exp2*exper*exper)
R
summary(weight)
1 - sum(weight_nI$residuals^2) / sum(bwght$bwght^2) #0.7402849
sum(weight_nI$fitted.values^2) / sum(bwght$bwght^2) #0.7402849
