Fit Verteilung auf empirische Daten

Question

Ich versuche, eine Beta-Verteilung auf ein Histogramm aus empirischen Daten anzupassen.

Das Problem, auf das ich stoße, ist, dass die angepasste Verteilung viel höher ist als die Balken im ursprünglichen Histogramm.

Die Originaldaten liegen außerhalb des Bereichs von [0,1]. Dies ist der Bereich, in dem die Beta-Verteilung ausgewertet werden kann. Daher skaliere ich die Originaldaten so, dass sie im Intervall [0,1] liegt.

Hier ist mein Code:

load("https://www.dropbox.com/s/c3psxx8jjbc20mo/data.Rdata?dl=0") 

    #create histogram with values normalized between 0 and 1 
    h <- hist((data-min(data))/(max(data)-min(data)),lty="blank",col="grey") 
    #normalize the density so the y-axis goes from 0 to 1 
    h$density <- h$counts/max(h$counts) 
    #plot the results 
    plot(h,freq=FALSE,cex.main=1,cex.axis=1,yaxt='n',ylim=c(0,1.5),col='grey',lty='blank',xaxt='n') 
    axis(2,at=seq(0,1,0.5),labels=seq(0,1,0.5)) 
    axis(1,at=seq(0,1,0.5),labels=seq(0,1,0.5)) 

    #fit beta distribution 
    a <- (data-min(data))/(max(data)-min(data)) 
    a[a==1] <- 0.9999 
    a[a==0] <- 0.0001 
    fit.beta <- suppressWarnings(fitdistr(a, "beta", start = list(shape1=0.1, shape2=0.1))) 

    #overlay curve from beta distribution 
    alpha <- fit.beta$estimate[1] 
    beta <- fit.beta$estimate[2] 
    b <- rbeta(length(data),alpha,beta) 
    lines(density(b)) 

Was bin ich?

Answer 1

Zuerst möchten Sie hist(..., freq=TRUE) für das Histogramm verwenden. Um den y-Achsenbereich korrekt festzulegen, können Sie den maximalen Wert Ihrer Beta-Verteilung berechnen (see e.g. here). Schließlich woul es viel besser, dbeta zu verwenden, als eine zufällige Stichprobe durch Schätzen der Dichte gefolgt erzeugen:

maxibeta <- dbeta((alpha-1)/(alpha+beta-2), alpha, beta) 
hist((data-min(data))/(max(data)-min(data)), 
     prob=TRUE, col="grey", border="white", ylim=c(0, maxibeta), 
     main="Histogram + fitted distribution") 
plot(function(x) dbeta(x,alpha,beta), add=TRUE, col=2, lwd=2) 

---

EDIT: eine allgemeinere Lösung, aber das macht mich ein wenig traurig

fbeta <- function(x) dbeta(x,alpha,beta) 
maxibeta <- optimize(fbeta, interval = c(0,1), maximum = TRUE)$objective 

histo <- hist((data-min(data))/(max(data)-min(data)), plot = FALSE) 

plot(histo, freq=FALSE, col="grey", border="white", 
    ylim=c(0, max(maxibeta, max(histo$density))), 
    main="Histogram + fitted distribution") 
plot(fbeta, add=TRUE, col=2, lwd=2)