Det är kanske inte så ofta man behöver integrera en funktion, dvs beräkna ytan under en kurva. I statistisk analys är det antagligen vanligast då man söker en sannolikhet i en kontinuerlig fördelning. I dylika fall fanns det ‘förr i tiden’ tabeller men nu används enkla kommandon i något statistikprogram. Koderna nedan visar dock med några exempel hur detta kan göras med R-kommandon.

För flera av de vanliga fördelningarna - men ej för t.ex. normalfördelningen - finns det relativt enkla formler för att beräkna väntevärde och varians, men för andra är man hänvisad till någon numerisk funktion via dator.

Nedan sker beräkningar både med integrering men också med de vanliga formlerna och naturligtvis blir resultatet detsamma (för åtminstone ett antal decimaler). Följande uttryck är generella för alla kontinuerliga fördelningar:

Väntevärdet (μ), varians (σ-kvadrat) samt ytan (F(z)) från x = 0 till x = z (Generellt gäller integreringen från negativ oändlighet och till z men här diskuteras fördelningar enbart på den positiva x-axeln):


\[ \begin{align} \large \mu =\int_{}^{} x\cdot f(x)dx & \phantom{litet tomt ut} \large \sigma^{2} = \int_{}^{} (x-\mu)^{2}\cdot f(x) \; dx & \phantom{litet tomt ut} \large F(z) = \int_{0}^{z} f(x) \; dx \end{align} \]


Nedan visas tre fall - en integrering av exponentialfördelningens funktion, en integrering av Weibullfördelningens funktion samt konvertering av en funktion till en sannolikhets- och fördelningsfunktion.


Detta dokument delas in i tre delar:


   1. Integrering av exponentialfördelningen

   2. Integrering av Weibullfördelningen

   3. Skapa och simulera en sannolikhets- och fördelningsfunktion





 1. Integrering av exponentialfördelningen

Exemplet nedan visar en vanlig exponentialfördelning vars egenskaper är välkända och uttrycken för väntevärde och varians är enkla:

\[ \begin{align} \large f(x) = \lambda e^{-\lambda x} & \phantom{litet tomt ut} \large \mu = \frac{1}{\lambda} & \phantom{litet tomt ut} \large \sigma^{2} = \frac{1}{\lambda^{2}} \end{align} \] 

library(ggplot2)

expfunktion <- function(x) {lam*exp(-lam*x)}            # Originalfunktionen.

lam    = 1.7
grans1 = 0.5
grans2 = 3.1

yta1 <- integrate(expfunktion, grans1, grans2)

yta  <- sprintf("%.3f", round(yta1$value, 3))
yta2 <- sprintf("%.3f", round(pexp(grans2, lam) - pexp(grans1, lam), 3))

text1 <- paste0(" Exponentialfördelningen")
text2 <- paste0(" Yta (integration) = ", yta)
text3 <- paste0(" Yta (R-kommandon) = ", yta2)

mingraf <- ggplot() + 
  xlim(0, 4) + 
  ylim(0, lam) +
  geom_vline(xintercept = c(grans1, grans2)) + 
  geom_function(fun = expfunktion, colour='blue', linewidth=1.5) +
  annotate("text", x=-Inf, y=Inf, label=text1, vjust = 2, hjust=-1, size=4.5, colour="blue") +
  annotate("text", x=-Inf, y=Inf, label=text2, vjust = 5, hjust=-1.2, size=4., colour="red") +
  annotate("text", x=-Inf, y=Inf, label=text3, vjust = 8, hjust=-1.02, size=4., colour="red") +
  xlab('x') + ylab("f(x)")
mingraf

Grafen ovan visar en exponentialfördelning med två gränser och med motsvarande yta. Det är som förväntat samma resultat med bägge beräkningsmetoderna, integrering och R-kommandon.


Beräkning av väntevärde och varians. Beräkningen genomförs både med de generella uttrycket för kontinuerliga fördelningar och formlerna för exponentialfördelningen. Nederst visas resultaten i en tabell: 


library(kableExtra) 

my        <- function(x) {x * expfunktion(x)}                    # Väntevärde 
myinter   <- integrate(my, 0, Inf)$value
vari      <- function(x){(x-myinter)^2 * expfunktion(x)}         # Varians
variinter <- integrate(vari, 0, Inf)$value

myinter   <- round(myinter, 3)
variinter <- round(variinter, 3)

medelvarde <- round(1/lam, 3)                        # Väntevärde o varians i en exp-fördelning.
varians    <- round((1/lam)*(1/lam), 3)

a <- c("", "Yta", "Väntevärde", "Varians")
b <- c("Integration", yta, myinter, variinter)
c <- c("R-kommandon", yta2, medelvarde, varians)

df <- data.frame(a, b, c)
names(df) <- NULL

df %>%
  kbl(booktabs = T, align = "l") %>%
  kable_styling( font_size = 17) %>%
  kable_paper(full_width = T) %>%
  column_spec(1, bold=TRUE) %>%
  column_spec(2, bold=TRUE) %>%
  column_spec(3, bold=TRUE)
Integration R-kommandon
Yta 0.422 0.422
Väntevärde 0.588 0.588
Varians 0.346 0.346




 2. Integrering av Weibullfördelningen

Exemplet nedan visar en vanlig Weibullfördelning vars egenskaper är välkända men uttrycken för väntevärde och varians är en smula komplicerade:


\[ \begin{align} \large \mu = \lambda \cdot \Gamma \left( 1+\frac{1}{k} \right) & \phantom{litet tomt ut} \large \sigma^{2} = \lambda^{2} \cdot \left( \Gamma \left( 1+\frac{2}{k} \right) -\Gamma \left( 1+\frac{1}{k}\right)^{2} \right) \end{align} \] 

library(ggplot2)

Weibfunktion <- function(x) {k/lam * (x/lam)^(k-1) * exp(-(x/lam)^k)}    # Originalfunktionen.

k      <- 1.5
lam    <- 3.7
grans1 <-1.1
grans2 <- 5.3

yta1 <- integrate(Weibfunktion, grans1, grans2)
yta  <- sprintf("%.3f", round(yta1$value, 3))
yta2 <- sprintf("%.3f", round(pweibull(grans2, k, lam) - pweibull(grans1, k, lam), 3))

text1 <- paste0(" Weibullfördelningen")
text2 <- paste0(" Yta (integration) = ", yta)
text3 <- paste0(" Yta (R-kommandon) = ", yta2)

mingraf <- ggplot() + 
  xlim(0, 14) + 
  geom_vline(xintercept = c(grans1, grans2)) + 
  geom_function(fun = Weibfunktion, colour='blue', linewidth=1.5) +
  annotate("text", x=-Inf, y=Inf, label=text1, vjust = 2, hjust=-1, size=4.5, colour="blue") +
  annotate("text", x=-Inf, y=Inf, label=text2, vjust = 5, hjust=-1.2, size=4., colour="red") +
  annotate("text", x=-Inf, y=Inf, label=text3, vjust = 8, hjust=-1.02, size=4., colour="red") +
  xlab('x') + ylab("f(x)")
mingraf

Grafen ovan visar en Weibullfördelning med två gränser och med motsvarande yta. Det är som förväntat samma resultat med bägge beräkningsmetoderna, integrering och R-kommandon.


Beräkning av väntevärde och varians. Beräkningen genomförs både med de generella uttrycket för kontinuerliga fördelningar och formlerna för Weibullfördelningen. Nederst visas resultaten i en tabell: 


library(kableExtra) 

my        <- function(x) {x * Weibfunktion(x)}              # Väntevärde 
myinter   <- integrate(my, 0, Inf)$value
vari      <- function(x){(x-myinter)^2 * Weibfunktion(x)}   # Varians
variinter <- integrate(vari, 0, Inf)$value

myinter   <- sprintf("%.3f", round(myinter, 3))
variinter <- sprintf("%.3f", round(variinter, 3))

myWeibull  <- lam*(gamma(1 + 1/k))
varWeibull <- lam*lam*(gamma(1 + 2/k) - (gamma(1 + 1/k))^2)

myWeibull  <- sprintf("%.3f", round(myWeibull, 3))          # Väntevärde o varians i en Weibull-fördelning.
varWeibull <- sprintf("%.3f", round(varWeibull, 3))

a <- c("", "Yta", "Väntevärde", "Varians")
b <- c("Integration", yta, myinter, variinter)
c <- c("R-kommandon", yta2, myWeibull, varWeibull)

df <- data.frame(a, b, c)
names(df) <- NULL

df %>%
  kbl(booktabs = T, align = "l") %>%
  kable_styling( font_size = 17) %>%
  kable_paper(full_width = T) %>%
  column_spec(1, bold=TRUE) %>%
  column_spec(2, bold=TRUE) %>%
  column_spec(3, bold=TRUE)
Integration R-kommandon
Yta 0.670 0.670
Väntevärde 3.340 3.340
Varians 5.143 5.143



 3. Skapa och simulera en sannolikhets- och fördelningsfunktion

Exemplet nedan visar i täljaren en trigonomisk funktion. Ytan under funktionen är 1.455 och såleden divideras funktionen med detta värde för att ytan skall bli 1 och kunna vara en sannolikhetsfördelning (ursprunget till funktionen redovisas inte här):


\[ \begin{align} \large f(x) = \frac{(sinx)^2 \cdot cosx+0.4} {1.455} \end{align} \] 

library(ggplot2)

# Antag att följande funktion skall utgöra en sannolikhetsfördelning:
#
#      sin(x)^2 * cos(x) + 0.4)
#
# Men ytan under kurvan är inte 1 utan först måste ytan beräknas och
# sedan divideras funktionen med ytan:

f <- function(x) {  (sin(x)^2 * cos(x) + 0.4) }             # Originalfunktionen.

grans1 = 0
grans2 = 4
area <- integrate(f, grans1, grans2)$value

fordelning <- function(x) {((sin(x)^2 * cos(x) + 0.4))/area}  # Sann.fördelning.

mingraf <- ggplot()
mingraf <- mingraf + xlim(grans1, grans2) +xlab('x') + ylab('f(x)')
mingraf <- mingraf + geom_function(fun = fordelning, colour='blue', linewidth=1.5) +
           annotate("text", x = 2, y = Inf, label=" En fördelning", vjust = 2, size = 5.5, colour = "blue")
mingraf

Grafen ovan visar fördelningen med två tydliga toppar.


Simulering av data. Data från de vanligaste sannolikhetsgfördelningarna kan simuleras med hjälp av relativt enkla kommandon. Vid ovanligare fördelningar anger litteraturen några olika sätt. Ett av dem kallas ‘The Rejection Method’ (se t.ex. boken SIMULATION av Sheldon M. Ross). Denna metod används nedan genom ett antal rader kod: 


# Hitta fördelningens högsta punkt. Skall användas för simulering av
# med den s.k. 'The Rejection Method'.
# Metoden behöver ytterligare en fördelning att användas vid simuleringen.
# Här används en likformig fördelning (g(x)) på intervallet [0, 4] och som
# sålunda har höjden 1/4 (0.25).

#    f(x)     (sin(x)^2 * cos(x) + 0.4))/area}
#   ------ = ----------------------------------  -> c-värdet (sid 66, SIMULATION)
#    g(x)               1/4

fxdivgx    <- function(x) {0.25*((sin(x)^2 * cos(x) + 0.4))/area}  # Sann.fördelning.
maxfxdivgx <- optimize(fxdivgx, interval = c(0, 4), maximum = TRUE)
c          <- as.numeric(maxfxdivgx[1])


# f(x)/(c*g(x)) används för att avgöra om ett slumptal skall anses tillhöra
# ovanstående 'fordelning'.

antslump <- 10000

U1 <- runif(antslump, 0, 4)
U2 <- runif(antslump, 0, 4)

z <- fordelning(U1)/(c*0.25)

resultat <- ifelse(U2 < z, U1, NA)          # Sid 66-67 i SIMULATION

data <- data.frame(U1, U2, z, resultat)

resultat <- resultat[!is.na(resultat)]

hist <- ggplot() + geom_histogram(aes(resultat, after_stat(density)), bins=20,color='black', alpha = 0.3) +
                   geom_function(fun = fordelning, colour='blue', linewidth=1.5)
hist

Ovanstående graf visar det simulerade utfallet tillsammans med fördelningen.


Beräkning av väntevärde och varians. Eftersom det inte finns enkla formler för beräkningarna måste den ursprungliga definition används, dvs en integrering av funktionen. Tabellen visar god överensstämmelse mellan teoretiska och simulerade värden: 


library(ggplot2)

# beräkna my och varians

totyta  <- integrate(fordelning, 0, 4)

myfun   <- function(x) {x*fordelning(x)}                         # my
#myvarde <- round(integrate(myfun, 0, 4)$value, 3)

myvarde <- integrate(myfun, 0, 4)$value

varfun  <- function(x){(x-myvarde)^2 * fordelning(x) }           # Variance
varians <- round(integrate(varfun, 0, 4)$value, 3)

myvarde <- sprintf("%.3f", round(myvarde, 3))
varians <- sprintf("%.3f", round(varians, 3))

simmedel   <- mean(resultat)
simvarians <- var(resultat)

simmedel   <- sprintf("%.3f", round(simmedel, 3))
simvarians <- sprintf("%.3f", round(simvarians, 3))

a <- c("", "Väntevärde", "Varians")
b <- c("Integration", myvarde, varians)
c <- c("Simulering", simmedel, simvarians)

df <- data.frame(a, b, c)
names(df) <- NULL

df %>%
  kbl(booktabs = T, align = "l") %>%
  kable_styling( font_size = 17) %>%
  kable_paper(full_width = T) %>%
  column_spec(1, bold=TRUE) %>%
  column_spec(2, bold=TRUE) %>%
  column_spec(3, bold=TRUE)
Integration Simulering
Väntevärde 1.520 1.542
Varians 1.189 1.187


Skapa fördelningsfunktionen - F(x). Fördelningsfunktionens Y-axel visar hur ‘sannolikhetsytan’ ökar då värdena på X-axeln ökas och funktionen är alltid växande. Ofta anpassar man mätdata till en tänkt fördelningsfunktion för att undersöka överensstämmelse mellan data och model (se t.ex. ‘normalfördelningstest’): 


library(ggplot2)

Fx     <- c()
Xskala <- c()

for (i in seq(0, 4, by = 0.04)) {
  svar   <- integrate(fordelning, 0, i)$value
  Fx     <- append(Fx, svar)
  Xskala <- append(Xskala, i) 
}

Fxdata <- data.frame(Fx, Xskala)

mingraf <- ggplot() +
           xlim(0, 4) +
           xlab('x') +
           ylab('F(x)') + 
           geom_point(data=Fxdata, aes(Xskala, Fx), size=1, col="blue") +
           geom_path(data=Fxdata,  aes(Xskala, Fx), color="blue") +
           annotate("text", x = 2, y = Inf, label=" F(x)", vjust = 2, size = 5.5, colour = "blue")
mingraf



Avslutningsvis
Ovanstående koder och texter avser inte att ge en uttömmande diskussion om integrering av en funktion. Men här har framställningen poängterat enkelhet.

(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.)