Bakgrund

I boken A History of Mathematical Statistics From 1750 to 1930 finns det på sidan 54 en intressant bild på en sannolikhetsfördelning, se nedan. Fördelningen gäller spridningen av medelvärdet av tre mätvärden, likformigt fördelade på intervallet [0, a].
Det är ganska uppenbart att det är något fel på kurvan. Den består av tre andragradskurvor - en från 0 till a/3, en från a/3 till 2a/3 och till sist en från 2a/3 till a på X-axeln. På sidan 53 finns det matematiska uttryck som ger tre korrekta kurvor, R-koder finns nedan, men det är lite otydligt varför diagrammet har ritats som på bilden.
























Nedan anges tre olika illustrationer:

   1. Laplaces uttryck och Bates fördelning

   2. En simulering av Bates fördelning

   3. Normalfördelningen som tre andragradsfunktioner



 1. Laplaces uttryck och Bates fördelning

Laplaces andragradsuttryck finns på sidan 53 i boken och är programmerade i kodavsnittet nedan. ‘Bates distribution’ (se t.ex. Wikipedia) handlar om medelvärdet av n värden från en likformig fördelning på intervallet [0, a].
För beräkning av Bates fördelning har jag skapat en funktion (Batesfxn) som tar emot ett antal parametrar, se kodningen. 


# Se Bates Distribution i Wikipedia

library(ggplot2)

Batesfxn <- function(n, x, a, b) {
           zsum <- c()
            for(k in c(0:n)) {
             x1 <- (x - a)/(b - a)
             z1 <- 1/(b-a) * n/(2*factorial(n-1)) * (-1)^k * choose(n, k) * (n*x1 - k)^(n-1) * sign(n*x1 - k)
            zsum <- c(zsum, z1)
  }
            return(sum(zsum))
}

n <- 3                           # Antal termer. OBS i denna övning är n bara tre (p.gr.a. Laplaces uttryck)
a <- 0                           # Undre och övre interval på den likformiga fördelningen.
b <- 1

xvarde <- seq(a, b, 0.01)        # Skapar en lista med x-värden.
yvarde <- c()                    # En kolumn som tar emot y-värdena (dvs f(x))

  for (i in c(1:length(xvarde))) {
       yvarde <- c(yvarde, Batesfxn(n, xvarde[i], a, b))
  }

# 'xvarde' och 'yvarde' är den korrekta sannolikhetsfördelningen, ritas nedan.
# Skapats enligt 'Bates distribution'
# --------------------------------------------------------------------

# Nedan skapas tre kurvor enligt uttrycken på sidan 53 i boken.

a <- 1
n <- 3

x1 <- seq(0, a/3, 0.01)
r <- 1
z <- x1 - (r-1)*a/n
y1 <- 27/(2*a^3)*(z)^2                           # Vänstra (blå) kurvsegmentet.

x2 <- seq(a/3, 2*a/3, 0.01)
r <- 2
z <- x2 - (r-1)*a/n
y2 <- 3/(2*a^3)*(a^2 + 6*a*(z) - 18*(z)^2)       # Röda kurvsegmentet.

x3 <- seq(2*a/3, a, 0.01)
r <- 3
z <- x3 - (r-1)*a/n
y3 <- 3/(2*a^3)*(a^2 - 6*a*(z) + 9*z^2)          # Högra (blå) kurvsegmentet.

fordelning <- ggplot() +
  geom_line(aes(x=xvarde, y=yvarde), linewidth=3, linetype=2, color=("green")) +   # Enligt Bates Distribution.
  geom_line(aes(x=x1, y=y1), linewidth=1.5, color=("blue")) +                      # Vänstra delen enl boken.
  geom_line(aes(x=x2, y=y2), linewidth=1.5, color=("red"))  +                      # Mittendel enl boken.
  geom_line(aes(x=x3, y=y3), linewidth=1.5, color=("blue")) +                      # Högra delen enl boken.
  xlab("x") + ylab("f(x)")
fordelning

Den gröna streckade linjen är korrekt sannolikhetsfördelning beräknad enligt Bates fördelning. De två blå och det röda avsnitten är andragradsfunktionerna enligt Laplaces uttryck enigt sidan 53 i boken. Alla avsnitten sammanfaller bra med den korrekta funktionen.




 2. En simulering av Bates fördelning

Nedan simuleras tre kolumner med värden från en likformig fördelning [0, 1] där medelvärdet beräknas radvis. Grafen innehåller ett histogram samt en kurva (som beräknats ovan): 


antal <- 5000

a <- 0
b <- 1

data1 <- runif(antal, a, b) 
data2 <- runif(antal, a, b) 
data3 <- runif(antal, a, b) 

dataSim<- data.frame(cbind(data1, data2, data3))
dataSim$medelv <- rowMeans(dataSim) 

fordelnkurva <- data.frame(xvarde, yvarde)

histogram <- ggplot(dataSim, aes(medelv)) + geom_histogram(bins=15, color='yellow', alpha = 0.5, 
             aes(y = after_stat(density)), position = 'identity') + labs(x = "Simulerade värden", y = "") +
             geom_line(data = fordelnkurva, aes(x = xvarde, y = yvarde), linewidth = 1.5, col = "green")
histogram

Kommentar till grafen. Histogrammet är som förväntat centrerat runt medelvärdet (0.5). Ju fler termer som ingår i medelvärdet desto mer normalfördelningslik blir fördelningen, dock med undantaget att svansarna inte går mot oändligheten.




 3. Normalfördelningen som tre andragradsfunktioner

Det verkar en smula onödigt att försöka anppassa tre andragradskurvor till normal­för­delningens sannolikhetsfunktion. Denna är ju välkänd och med många datorkommandon till hands.
X-axeln delades upp i tre sektioner och för varje sektion beräknas f(x). Därefter används glm-kommandot för att anpassa en andragradskurva till värdena. De erhållna parametervärdena utnyttjas sedan för att rita kurvan: 


library(ggplot2)

xaxel1 <- seq(-3, -1, 0.2)                  # Left part of X-axis.
yaxel1 <- dnorm(xaxel1)                     # Calculating f(x).

ana1 <- glm(yaxel1 ~ xaxel1 + I(xaxel1^2))  # Performing the analysis.
b01 <- ana1$coef[1]                         # Extracts the three coefficients.
b11 <- ana1$coef[2]
b21 <- ana1$coef[3]

xaxel2 <- seq(-1, 1, 0.2)                   # Middle part of the X-axis.
yaxel2 <- dnorm(xaxel2)                     # Calculating f(x).

ana2 <- glm(yaxel2 ~ xaxel2 + I(xaxel2^2))  # Performing the analysis.
b02 <- ana2$coef[1]                         # Extracts the three coefficients.
b12 <- ana2$coef[2]
b22 <- ana2$coef[3]

xaxel3 <- seq(1, 3, 0.2)                    # Right part of the X-axis.
yaxel3 <- dnorm(xaxel3)                     # Calculating f(x).

ana3 <- glm(yaxel3 ~ xaxel3 + I(xaxel3^2))  # Performing the analysis.
b03 <- ana3$coef[1]                         # Extracts the three coefficients.
b13 <- ana3$coef[2]
b23 <- ana3$coef[3]

del1 <- function(x) {b01 + b11*x + b21*x*x} # Creating functions of the three parts.
del2 <- function(x) {b02 + b12*x + b22*x*x}
del3 <- function(x) {b03 + b13*x + b23*x*x}
del4 <- function(x) {dnorm(x)}              # f(x) according to 'R'-code.

approxNormal <- ggplot() +
            geom_function(fun=del1, colour='blue',  linewidth=1.5, xlim=c(-3, -1)) +            # Left 2-degree
            geom_function(fun=del2, colour='red',   linewidth=1.5, xlim=c(-1,  1)) +            # Middle 2-degree
            geom_function(fun=del3, colour='blue',  linewidth=1.5, xlim=c( 1,  3)) +            # Right 2-degree
            geom_function(fun=del4, colour='green', linewidth=1.5, linetype=2, xlim=c(-3, 3)) + # normal curve
            xlab("x") + ylab("f(x)")
approxNormal

Kommentar till grafen. Grön linje är den korrekta kurvan men med ögonmått passar andragradskurvorna bra. Alla andragradskurvor är ju parabler och vid -3 and +3 på X-axeln syns kurvorna vända uppåt.



Avslutningsvis. Ovanstående text är föranledd av en bild i boken som omnämnts och innehåller i sak inget nytt.

(Analyserna har gjorts med datorprogrammet ‘R’ och med det grafiska gränssnittet ‘R-studio’ och bägge är gratis tillgängliga på nätet. Se https://www.indstat.se och knappen [Statistikprogram - R] för installation.)

(Se https://www.indstat.se för många andra simuleringsövningar.)