|
|
# Rootogramm
##########################################################################
data(buffalo,package="gss")
a <- hist(buffalo,10)
a1 <- density(buffalo)
n <- length(buffalo)
a11 <- dnorm(a1$x,mean=mean(buffalo),sd=sd(buffalo))
a1$y <- a11
# schaetze Histogramm
binw <- a1$x[2]-a1$x[1] # Breite der aequidistanten Intervalle
xhist <- a1$y*binw
xhist <- xhist/sum(xhist) # gesamte Flaeche 1
# rekonstruiertes Histogramm:
k <- length(a$breaks)
xrek <- rep(0,k+1)
for (i in 1:(k+1)){
if (i==1){xrek[i] <- sum(xhist[a1$x
if (i==(k+1)){xrek[i] <- sum(xhist[a1$x>=a$breaks[k]])}
if ((i>1) & (i<(k+1))){xrek[i] <- sum(xhist[(a1$x>=a$breaks[i-1]) & (a1$x
}
# plot sqrt histogram
pdf("histdicht1.pdf",width=7,height=6)
par(mar=c(4.5,4.5,1,1))
b1 <- a$mids[2]-a$mids[1]
xdens <- c(a$mids[1]-b1,a$mids,a$mids[length(a$mids)]+b1)
plot(0,0,xlim=c(xdens[1],xdens[length(xdens)]),ylim=c(0,max(sqrt(a$counts))),type="n",
xlab="Schneefall-Daten",ylab="Wurzel der Häufigkeiten",cex.lab=1.3)
lines(xdens,sqrt(xrek*n),lty=2)
for (i in 1:length(a$counts)){
x1 <- c(a$mids[i]-b1/2,a$mids[i]+b1/2,a$mids[i]+b1/2,a$mids[i]-b1/2,a$mids[i]-b1/2)
y1 <- c(0,0,sqrt(a$counts[i]),sqrt(a$counts[i]),0)
polygon(x1,y1)
}
dev.off()
|
|
|
# Rootogramm
##########################################################################
data(buffalo,package="gss")
a <- hist(buffalo,10)
a1 <- density(buffalo)
n <- length(buffalo)
a11 <- dnorm(a1$x,mean=mean(buffalo),sd=sd(buffalo))
a1$y <- a11
# schaetze Histogramm
binw <- a1$x[2]-a1$x[1] # Breite der aequidistanten Intervalle
xhist <- a1$y*binw
xhist <- xhist/sum(xhist) # gesamte Flaeche 1
# rekonstruiertes Histogramm:
k <- length(a$breaks)
xrek <- rep(0,k+1)
for (i in 1:(k+1)){
if (i==1){xrek[i] <- sum(xhist[a1$x
if (i==(k+1)){xrek[i] <- sum(xhist[a1$x>=a$breaks[k]])}
if ((i>1) & (i<(k+1))){xrek[i] <- sum(xhist[(a1$x>=a$breaks[i-1]) & (a1$x
}
# plot sqrt histogram reverse
pdf("histdicht2.pdf",width=7,height=6)
par(mar=c(4.5,4.5,1,1))
plot(0,0,xlim=c(xdens[1],xdens[length(xdens)]),ylim=c(max(sqrt(a$counts)),-4),type="n",
xlab="Schneefall-Daten",ylab="Wurzel der Häufigkeiten",cex.lab=1.3)
lines(xdens,sqrt(xrek*n),lty=2)
acount <- c(0,sqrt(a$counts),0)
xrek1 <- sqrt(xrek*n)
amids <- c(a$mids[1]-b1,a$mids,a$mids[length(a$mids)]+b1)
for (i in 1:length(acount)){
x1 <- c(amids[i]-b1/2,amids[i]+b1/2,amids[i]+b1/2,amids[i]-b1/2,amids[i]-b1/2)
y1 <- c(xrek1[i],xrek1[i],xrek1[i]-acount[i],xrek1[i]-acount[i],xrek1[i])
polygon(x1,y1)
}
abline(0,0,lty=3)
dev.off()
|