Str. 13

Rozkłady prawdopodobieństwa w programie R – wybrane przykłady

?Distributions

Rozkład dwumianowy

?Binomial

Prawdopodobieństwo uzyskania określonej liczby sukcesów w n próbach Bernoulliego

Obliczenia realizuje funkcja dbinom(x,size,prob), gdzie x=liczba sukcesów, size=liczba prób, prob=prawdopodobieństwo sukcesu w pojedynczej próbie:

Przykład. Niech X~Bin(10,1/2). Obliczyć P(X=2), czyli prawdopodobieństwo, że w 10 próbach będą 2 sukcesy, jeśli prawdopodobieństwo sukcesu w jednej próbie wynosi 0,5.

dbinom(x=2, size=10,prob=0.5)

„Hurtowe” obliczenia prawdopodobieństw dla każdej możliwej liczby sukcesów, jaką można uzyskać w 10 próbach:

c(0,1,2,3,4,5,6,7,8,9,10) # kolejne liczby całkowite od 0 do 10

#to samo można uzyskać za pomocą funkcji seq()

seq(from=0,to=10,by=1)

#to samo można uzyskać za pomocą operatora dwukropek

0:10 # kolejne liczby całkowite od 0 do 10

dbinom(x=0:10, size=10,prob=0.5)

#wynik obliczeń można zapisać jako obiekt (wektor) np. o nazwie y

y<-dbinom(x=0:10, size=10,prob=0.5)

y

Funkcja plot() służy do wykonywania różnych wykresów na płaszczyźnie

plot(x=0:10,y=y)

Funkcja plot() ma wiele opcji, zob. ?plot oraz ?plot.default

plot(x=0:10,y=y, type="h",xlab="Liczba sukcesów", ylab="Prawdopodobieństwo")

Kształt rozkładu dwumianowego zależy od wielkości prawdopodobieństwa „sukcesu” w pojedynczej próbie:

plot(x=0:10,y=y, type="h",xlab="Liczba sukcesów", ylab="Prawdopodobieństwo",main="p=0,5")

y<-dbinom(x=0:10, size=10,prob=0.25)

plot(x=0:10,y=y,type="h",xlab="Liczba sukcesów", ylab="Prawdopodobieństwo",main="p=0,25")

y<-dbinom(x=0:10, size=10,prob=0.75)

plot(x=0:10,y=y,type="h",xlab="Liczba sukcesów", ylab="Prawdopodobieństwo",main="p=0,75")

Można spowodować, żeby wykresy pojawiły się w jednym oknie graficznym podzielonym na części:

#wykresy będą umieszczone kolejno w wierszach macierzy 2x2

par(mfrow=c(2,2))

y<-dbinom(x=0:10, size=10,prob=0.5)

plot(x=0:10,y=y,type="h",xlab="Liczba sukcesów", ylab="Prawdopodobieństwo",main="p=0,5")

y<-dbinom(x=0:10, size=10,prob=0.25)

plot(x=0:10,y=y,type="h",xlab="Liczba sukcesów", ylab="Prawdopodobieństwo",main="p=0,25")

y<-dbinom(x=0:10, size=10,prob=0.75)

plot(x=0:10,y=y,type="h",xlab="Liczba sukcesów", ylab="Prawdopodobieństwo",main="p=0,75")

par(mfrow=c(1,1))#powrót do ustawienia domyślnego okna graficznego

Dystrybuanta rozkładu dwumianowego

Obliczenia realizuje funkcja pbinom(q, size, prob), gdzie q=argument, dla którego obliczana jest dystrybuanta, size=liczba prób, prob=prawdopodobieństwo sukcesu w pojedynczej próbie

Przykład.

Niech X~Bin(10,1/2). Obliczyć wartość dystrybuanty tego rozkładu w punkcie x=2,7.

Niech F oznacza dystrybuantę rozkładu Bin(10,1/2). Trzeba zatem obliczyć wartość funkcji F w punkcie 2,7, czyli F(2,7).

pbinom(q=2.7, size=10, prob=0.5)

Bezpośrednio z definicji dystrybuanty wynika, że F(2,7) = P(X ≤ 2,7) = P(X=0)+P(X=1)+P(X=2).

dbinom(x=0:2, size=10, prob=0.5)

sum(dbinom(x=0:2, size=10, prob=0.5))

Generowanie liczb losowych z rozkładu dwumianowego

rbinom(n, size, prob), gdzie n=liczba wylosowanych wartości, size=liczba prób (liczba powtórzeń próby Bernoulliego), prob=prawdopodobieństwo sukcesu w pojedynczej próbie

rbinom(n=5, size=10, prob=0.5)

set.seed(112)#zapewnia powtarzalność wygenerowanych liczb losowych

rbinom(n=5, size=10, prob=0.5)

Rozkład Poissona

?Poisson

Przykład.

Niech X~Poiss(3). Obliczyć P(X=4).

dpois(x=4, lambda=3)

Niech X~Poiss(3). Obliczyć P(X=0), P(X=1), …, P(X=20)

dpois(x=0:20, lambda=3)

y<-dpois(x=0:20, lambda=3)

plot(x=0:20,y=y,type="h",xlab="Wartości zmiennej losowej", ylab="Prawdopodobieństwa")

Niech X~Poiss(3). Obliczyć F(15) i F(-10) (wartość dystrybuanty w punkcie 15 i w punkcie -10):

ppois(q=15, lambda=3)

ppois(q=-10, lambda=3)

Generowanie liczb losowych z rozkładu Poissona

rpois(n, lambda), gdzie n=liczba wylosowanych wartości, lambda=parametr rozkładu

rpois(n=4,lambda=0.5)

Rozkład wykładniczy

?Exponential

Wyznaczanie wartości funkcji gęstości rozkładu Exp(lambda) - funkcja dexp(x, rate)

gdzie x=liczba, dla której należy wyznaczyć wartość funkcji gęstości , rate=lambda

#gęstość rozkładu Exp(2) w punkcie 0

dexp(x=0, rate = 2)

dexp(x=0.5, rate = 2)

Wykresy gęstości dla różnych wartości lambda

curve(expr=dexp(x, rate = 2),from=0,to=10,ylab="Gęstość", main="Gęstość rozkładu wykładniczego")

curve(expr=dexp(x, rate = 1),from=0,to=10,add=T,col="blue")

curve(expr=dexp(x, rate = 0.5),from=0,to=10,add=T,col="red")

legend("topright",legend=c("lambda=2","lambda=1","lambda=0,5"),

lty=1,col=c("black","blue","red"))

Sprawdzenie ścieżki do katalogu roboczego programu R:

getwd()

Zapis rysunku w formacie jpeg (można zapisywać także w formacie bmp, png, tiff, pdf, eps)

jpeg(file="wykladniczy.jpeg")# zapis pliku do katalogu roboczego

curve(expr=dexp(x, rate = 2),from=0,to=10,ylab="Gęstość", main="Gęstość rozkładu wykładniczego")

curve(expr=dexp(x, rate = 1),from=0,to=10,add=T,col="blue")

curve(expr=dexp(x, rate = 0.5),from=0,to=10,add=T,col="red")

legend("topright",lty=1, col=c("black","blue","red"),

legend=c(expression(lambda==2),expression(lambda==1),

expression(lambda==0.5)))

dev.off()

Generowanie liczb losowych z rozkładu wykładniczego

rexp(n, rate), gdzie n=liczba wylosowanych wartości, rate=parametr rozkładu

rexp(n=6,rate=0.5)

Rozkład normalny

?Normal

Wyznaczanie wartości funkcji gęstości rozkładu N(mi, sigma^2) – funkcja
dnorm(x, mean, sd), gdzie x=liczba, dla której należy wyznaczyć wartość funkcji gęstości , mean=mi, sd=sigma

#gęstość rozkładu N(0,1) w punkcie 0

dnorm(x=0,mean=0,sd=1)

dnorm(x=c(-1,0,1),mean=0,sd=1)

curve(expr=dnorm(x,mean=0,sd=1),from=-4,to=4, ylab="Gęstość", main="Gęstość rozkładu N(0,1)")

curve(expr=dnorm(x,mean=0,sd=1),from=-4,to=4,ylab="Gęstość", main="Gęstość rozkładu normalnego \n w zależności od parametrów")

curve(expr=dnorm(x,mean=0,sd=2),from=-4,to=4, ylab="Gęstość", add=T,lty=2)

legend("topright",lty=c(1,2),legend=c(expression(paste(mu==0, ", ",sigma==1)),expression(paste(mu==0, ", ",sigma==2))))

Wyznaczanie dystrybuanty: funkcja pnorm(q, mean, sd), gdzie q=liczba, dla której należy wyznaczyć wartość dystrybuanty , mean=mi, sd=sigma

Przykład. Niech X~N(3,2^2). Wyznaczyć:

P(X < 2)

pnorm(q=2, mean=3,sd=2)

P(X > 1)

1-pnorm(q=1, mean=3,sd=2)

#inaczej:

pnorm(q=1, mean=3,sd=2,lower.tail=FALSE)

Wyznaczanie kwantyli rozkładu: funkcja qnorm(p, mean, sd), gdzie p=rząd kwantyla , mean=mi, sd=sigma

#wyznaczenie kwantyla rzędu 0,5 rozkładu N(0,1)

qnorm(p=0.5,mean=0,sd=1)

#wyznaczenie kwantyla rzędu 0,975 rozkładu N(0,1)

qnorm(p=0.975,mean=0,sd=1)

jednoczesne wyznaczenie kwantyli rzędów 0,9; 0,95; 0,975; 0,99; 0,995 rozkładu N(0,1)

qnorm(p=c(0.9, 0.95, 0.975, 0.99, 0.995), mean=0, sd=1)

Nadanie nazw poszczególnym współrzędnym wektora

kwantyle<-qnorm(p=c(0.9, 0.95, 0.975, 0.99, 0.995), mean=0, sd=1)

names(kwantyle)<-paste("kwantyl rzędu",c(0.9, 0.95, 0.975, 0.99, 0.995))

kwantyle

Generowanie liczb losowych z rozkładu normalnego

rnorm(n, mean, sd), gdzie n=liczba wylosowanych wartości, mean=mi, sd=sigma

rnorm(n=4, mean=10, sd=2)

set.seed(120)

losowe<-rnorm(n=1000,mean=10,sd=2)

#wyświetlenie 10 pierwszych wartości

losowe[1:10]

summary(losowe)

sd(losowe)

hist(losowe)

hist(losowe,prob=TRUE)

lines(density(losowe))#gęstość empiryczna

x<-seq(3,20,by=0.1)

hist(losowe,prob=TRUE)

lines(density(losowe))#gęstość empiryczna

lines(x=x,y=dnorm(x=x,mean=10,sd=2),col="red",lwd=2)#gęstość teoretyczna

Rozkład t-Studenta

?TDist

df=liczba stopni swobody

Funkcja gęstości dt(x, df, ncp, log = FALSE)

Dystrybuanta pt(q, df, ncp, lower.tail = TRUE, log.p = FALSE)

Kwantyl qt(p, df, ncp, lower.tail = TRUE, log.p = FALSE)

Liczby losowe rt(n, df, ncp)

#kwantyl rzędu 0,975 dla rozkładu t-Studenta z 10 stopniami swobody

qt(p=0.975,df=10)

Rozkład chi-kwadrat

? Chisquare

df=liczba stopni swobody

Funkcja gęstości dchisq(x, df, ncp=0, log = FALSE)

Dystrybuanta pchisq(q, df, ncp=0, lower.tail = TRUE, log.p = FALSE)

Kwantyl qchisq(p, df, ncp=0, lower.tail = TRUE, log.p = FALSE)

Liczby losowe rchisq(n, df, ncp=0)

#kwantyl rzędu 0,95 dla rozkładu chi-kwadrat z 8 stopniami swobody

qchisq(p=0.95,df=8)

Rozkład F-Snedecora

? FDist

df1=liczba stopni swobody „licznika”

df2=liczba stopni swobody „mianownika”

Funkcja gęstości df(x, df1, df2, ncp, log = FALSE)

Dystrybuanta pf(q, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)

Kwantyl qf(p, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)

Liczby losowe rf(n, df1, df2, ncp)

#kwantyl rzędu 0,9 dla rozkładu F-Snedecora z df1=3,df2=8

qf(p=0.9,df1=3,df2=8)

Rozkład logarytmiczno-normalny

?Lognormal

Gęstość rozkładu

The mean is E(X) = exp(μ + 1/2 σ^2), the median is med(X) = exp(μ), and the variance Var(X) = exp(2*μ + σ^2)*(exp(σ^2) - 1) and hence the coefficient of variation is sqrt(exp(σ^2) - 1) which is approximately σ when that is small (e.g., σ < 1/2).

Gęstość rozkładu logarytmiczno-normalnego

dlnorm(x, meanlog, sdlog)

curve(expr=dlnorm(x,meanlog=0,sdlog=1),from=0,to=8,ylab="Gęstość")

curve(expr=dlnorm(x,meanlog=2,sdlog=0.5),from=0,to=40,ylab="Gęstość")

curve(expr=dlnorm(x,meanlog=4,sdlog=2),from=0,to=40,ylab="Gęstość",lty=2,

add=T)