# ......................................................................................
# .........................Cvičení 5 - Spojitá náhodná veličina.........................
# ..................Pavel Jahoda, vytvořeno na základě studijních.......................
#..........materiálů Martina Litschmannová, Adéla Vrtková, Michal Béreš.................
# ......................................................................................

# Nezobrazuje-li se vám text korektně, nastavte File \ Reopen with Encoding... na UTF-8
# Pro zobrazení obsahu skriptu použijte CTRL+SHIFT+O
# Pro spouštění příkazů v jednotlivých řádcích použijte CTRL+ENTER

# **Obsah tohoto skriptu je pouze jako doplňující ilustrace k cvičení, není nutno znát
# ke zkoušce. Důležité je to umět spočítat ručně.**
# * Numerická integrace v Rku ####
# Rkovská funkce **integrate**
# integrate(f, a, b) = $\int_{a}^{b}f(x)dx$
# - **f** je Rková funkce (námi definovaná) která má jeden vstupní argument - vektor
# hodnot ve kterém má vrátit své hodnoty
# - **a** dolní integrační mez
# - **b** horní integrační mez


f = function(x){return(x*x)} # x^2
a = -1
b = 2
integrate(f, a, b)



#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 1.Příklad 
# Náhodná veličina $X$ je dána hustotou pravděpodobnosti f:
# $f(x)=\begin{cases}
# 0    &      x < -1 \\
# c(1-x^2) &  -1 \leq x \leq 1 \\
# 0    &  1 < x
# \end{cases}$

#-----------------------------------------------------------------------------------
# \item[a)] Určete konstantu c.
c = 1

f = function(x){
  f = c*(1 - x^2)     # jinde
  f[x< -1] = 0 # 0 pro x< -1, pozor na x<-1 protože '<-' je v rku přiřazení
  f[x>1] = 0  # 0 pro x>1
  return(f)
} #ALT ENTER

integrate(f, -1, 1) # pro c = 1 vychází 4/3, nelze! musí vyjít 1!!!

c = 3/4

integrate(f, -1, 1) # pro c = 3/4 vychází 1, tak jak má!!!


#-----------------------------------------------------------------------------------
# \item[b)] Načrtněte graf hustoty pravděpodobnosti f.

x = seq(from = -5, to = 5, by = 0.01) # tady určíme pro která x spočítáme f(x)
fx = f(x) # tady jsme vytvořili vektor hodnot f(x) pro výše zvolená x
plot(x,fx, type = 'l') # vykreslit jako čáru


#-----------------------------------------------------------------------------------
# \item[c)] Určete pravděpodobnost, $P(-1/2 < X < 1/2)$

integrate(f, -0.5, 0.5) # = 0.6875 přesně


#-----------------------------------------------------------------------------------
# \item[d)] Určete pravděpodobnost, $P(0 < X)$

integrate(f, 0, 1) # = 0.5 přesně

#-----------------------------------------------------------------------------------
# \item[e)] Určete předpis distribuční funkce $F(x))$.

F = function(x){
  F = -(1/4)*(x^3) + (3/4)*(x) + 1/2    # jinde
  F[x< -1] = 0 # 0 pro x< -1, pozor na x<-1 protože '<-' je v rku přiřazení
  F[x>1] = 1  # 0 pro x>1
  return(F)
} #ALT ENTER

#Načrtněte graf distribuční funkce F.

x = seq(from = -5, to = 5, by = 0.01) # tady určíme pro která x spočítáme F(x)
Fx = F(x) # tady jsme vytvořili vektor hodnot F(x) pro výše zvolená x
plot(x,Fx, type = 'l') # vykreslit jako čáru



# srovnání hustoty pravděpodobnosti a distribuční funkce
x = seq(from = -5, to = 5, by = 0.01)  # hodnoty ve kterých budeme počítat f(x) a F(x)
fx = f(x)
plot(x, fx,  ylab='',ylim=c(0,1),main="f(x) a F(x)", col=c("red"))
grid()
par(new=TRUE) # že chceme kreslit do jednoho grafu
Fx = F(x)
plot(x, Fx, pch=4,  ylab='',ylim=c(0,1),main="f(x) a F(x)", col=c("blue"))



#-----------------------------------------------------------------------------------
# \item[f)] Určete pravděpodobnost $P(0,1 < X < 3)$.

# P(0,1 < X < 3) = F(3) - F(0.1) = 

F(3) - F(0.1)

#NEBO podle předpisu F:
1 - (-(1/4)*(0.1^3) + (3/4)*(0.1) + 1/2 )


#-----------------------------------------------------------------------------------
# \item[g)] Určete střední hodnotu náhodné veličiny X (značíme EX)


# E(X) = integrál od -infty do infty z funkce x*f(x)
xf = function(x){
  f = f(x)
  return(x*f)
} 
# stačí integrovat jen tam, kde víme, že je f(x) nenulová
integrate(xf, -1, 1)

#-----------------------------------------------------------------------------------
# \item[h)] Určete rozptyl náhodné veličiny X (značíme DX)


# D(X) = integrál od -infty do infty z funkce (x-EX)^2 *f(x)

EX = 0

# E(X) = integrál od -infty do infty z funkce x*f(x)
d = function(x){
  f = f(x)
  return((x-EX)^2*f)
} 
# stačí integrovat jen tam, kde víme, že je f(x) nenulová
integrate(d, -1, 1) # = 0.2 přesně

#-----------------------------------------------------------------------------------
# \item[i)] Určete směrodatnou odchylku náhodné veličiny X (značíme $\sigma$)

# sigma = odmocnina z rozptylu
sqrt(0.2) # = 0.4472136







#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 2.Příklad 
# Náhodná veličina $X$ je dána distribuční funkcí $F$:
# $F(x)=\begin{cases}
# 0             &         x < 2 \\
# -x^2 + 6x - 8 &  2 \leq x \leq 3 \\
# 1             &     3 < x
# \end{cases}$

#-----------------------------------------------------------------------------------
# \item[a)] Určete pravděpodobnost $P(1 < X < \frac{5}{2})$


F2 = function(x){
  F2 = -(x^2) + 6*(x) - 8 # jinde
  F2[x<2] = 0 # 0 pro x < 2
  F2[x>3] = 1  # 1 pro x > 3
  return(F2)
} #ALT ENTER

# P(1 < X < \frac{5}{2}) = 

Fx(5/2)


#-----------------------------------------------------------------------------------
# \item[b)] Určete předpis hustoty pravděpodobnosti NV $X$

# f je derivací F

f2 = function(x){
  f2 = -2*x + 6 # jinde
  f2[x<2] = 0 # 0 pro x < 2
  f2[x>3] = 0  # 1 pro x > 3
  return(f2)
} #ALT ENTER


# srovnání hustoty pravděpodobnosti a distribuční funkce
x = seq(from = -5, to = 5, by = 0.01)  # hodnoty ve kterých budeme počítat f(x) a F(x)
F2x = F2(x)
plot(x, F2x, type = "l",lwd=5,  ylab='',ylim=c(0,2),main="f(x) a F(x)", col=c("blue"))
grid()
par(new=TRUE) # že chceme kreslit do jednoho grafu
f2x = f2(x)
plot(x, f2x, type = "l",lwd=2,  ylab='',ylim=c(0,2),main="f(x) a F(x)", col=c("red"))



#-----------------------------------------------------------------------------------
# \item[c)] Určete pravděpodobnost $P( X < \frac{12}{5})$


integrate(f2, 2, 12/5) # 16/25 = 0,64

#NEBO:

F2(12/5) # 16/25 = 0,64

#-----------------------------------------------------------------------------------
# d) Určete střední hodnotu náhodné veličiny $X$.

# $EX = \int_{-infty}^{\infty} x.f(x) dx =  \int_{2}^{3} x.f(x) dx $

xf = function(x){
  xf = -2*x^2 + 6*x     #
  xf[x < 2] = 0 # 0 pro x<2
  xf[x > 3] = 0  # 1 pro x>3
  return(xf)
}

integrate(xf,2,3) # = 7/3


#-----------------------------------------------------------------------------------
# e) Určete rozptyl náhodné veličiny $X$.

# $DX = E(X^2) - (EX)^2 $

xxf = function(x){
  xxf = x^2*(-2*x+6 )     #
  xxf[x < 2] = 0 # 0 pro x<2
  xxf[x > 3] = 0  # 1 pro x>3
  return(xxf)
}

integrate(xxf,2,3) # = 11/2 = E(X^2)

11/2 -(7/3)^2 # = 1/18 


# NEBO $DX = \int_{-\infty}^{\infty} (x-EX)^2 f(x) dx$

ef = function(x){
  ef = (x -(7/3))^2*(-2*x+6 )     #
  ef[x < 2] = 0 # 0 pro x<2
  ef[x > 3] = 0  # 1 pro x>3
  return(ef)
}

integrate(ef,2,3) # = 1/18 = DX










###DALE STARE, PREDELAT################################
######################################################

#######################################################
# derovací F(x) získáme hustotu pravd. f(x)
# příslušná hustota pravděpodobnosti na intrvalu <0,1>
f = function(x){return(2*x)} # f(x) = 2x
a = 0
b = 1
integrate(f, a, b)$value

# c = 1, proto distribuční funkce vypadá takto:
F.dist = function(x){
    res = x*x     # x^2
    res[x<=0] = 0 # 0 pro x<=0
    res[x>1] = 1  # 1 pro x>1
    return(res)
}

x = seq(from = -1, to = 2, by = 0.01) # body na ose x
FX = F.dist(x) # hodnoty F(x)
plot(x, FX, type = 'l') # vykreslit jako čáru

# * Příklad 2. ####
# Rozdělení náhodné veličiny X je dáno hustotou
# $f(x)=\begin{cases}
# 2x+2 & x \in <-1;0> \\
# 0    & x \notin <-1;0>
# \end{cases}$
# Určete:
# ** 2. a) ####
# $F(x)$,


f.dens = function(x){
    res = 2*x + 2 
    # pozor na x<-1 protože '<-' je v rku přiřazení
    res[x < -1] = 0 # 0 pro x<=0
    res[x > 0] = 0  # 1 pro x>1
    return(res)
}

x = seq(from = -2, to = 1, by = 0.01) # body na ose x
fx = f.dens(x) # hodnoty f(x)
plot(x, fx, cex=0.2) # vykreslit tečky (cex je velikost)

F.dist = function(x){
    res = x*x+2*x+1     # x^2+2x+1
    res[x < -1] = 0 # 0 pro x<=0
    res[x > 0] = 1  # 1 pro x>1
    return(res)
}

x = seq(from = -2, to = 1, by = 0.01) # body na ose x
FX = F.dist(x) # hodnoty f(x)
plot(x, FX, type='l') # vykreslit tečky (cex je velikost)

# ** 2. b) ####
# P(−2 ≤ X ≤ −0.5), P(−2 ≤ X ≤ −1), P(X > 0.5), P(X = 0.3)


# P(−2 ≤ X ≤ −0.5)
integrate(f.dens, -2, -0.5)$value
integrate(f.dens, -1, -0.5)$value

# P(−2 ≤ X ≤ −1)
integrate(f.dens, -2, -1)$value

# P(X > 0.5)
integrate(f.dens, 0.5, 1e16)$value # tohle nebude vždy fungovat

# P(X = 0.3)
integrate(f.dens, 0.3, 0.3)$value
# je jasné že tato pravděpodobnost je 0
# odpovídá integrálu s a=b tedy s nulovou velikostí na ose x

# ** 2. c) ####
# střední hodnotu, rozptyl a směrodatnou odchylku náhodné veličiny X.


# E(X)
x_fx = function(x){
    fx = f.dens(x)
    return(x*fx)
} 
# integrujeme jen tam kde víme, že je f(x) nenulová
E_X = integrate(x_fx, -1, 0)$value
E_X
-1/3

# E(X^2)
xx_fx = function(x){
    fx = f.dens(x) 
    return(x*x*fx)
} 
# integrujeme jen tam kde víme, že je f(x) nenulová
E_XX = integrate(xx_fx, -1, 0)$value
E_XX
1/6

# D(X)
D_X = E_XX - E_X^2
D_X
1/18

# sigma(x)
std_X = sqrt(D_X)
std_X
sqrt(2)/6

# ** 2. d) ####
# modus $\hat{x}$


# modus = 0 

# ** 2. e) ####
# medián $x_{0,5}$


x = seq(from = -2, to = 1, by = 0.001) # body na ose x
FX = F.dist(x)
plot(x, FX, type='l')
lines(c(-2, 1),c(0.5, 0.5))

x[FX >= 0.5][1] # první prvek z x pro který F(x)>=0.5
(-2+sqrt(2))/2

# * Příklad 3. ####
# Náhodná veličina Y je definována jako: Y = 3X+1, kde X je náhodná veličina z
# předcházejícího
# příkladu. Určete:
# ** 3. a) ####
# $F_Y(y)$


FY.dist = function(y){
    # spočteno ze vztahu FY(y) = P(Y < y) = P(3X + 1 < y) = ...
    x = (y-1)/3 
    FY = F.dist(x)
    return(FY)
}
y = seq(from = -3, to = 2, by = 0.001) # body na ose x
FY = FY.dist(y)
plot(y, FY, type='l')

# ** 3. b) ####
# $f_Y(y)$


# derivace F_Y
fY.dens = function(y){
    res = 2/9*(y + 2) 
    res[y < -2] = 0 # 0 pro x<-2
    res[y > 1] = 0  # 1 pro x>1
    return(res)
}
integrate(fY.dens,-2,1)$value # kontrola celkového integrálu
y = seq(from = -3, to = 2, by = 0.001) # body na ose x
fY = fY.dens(y)
plot(y, fY, cex=0.2)

# ** 3. c) ####
# E(Y), D(Y), σ(Y)


# E(Y)
y_fy = function(y){
    fy = fY.dens(y)
    return(y*fy)
} 
# integrujeme jen tam kde víme, že je f(y) nenulová
E_Y = integrate(y_fy, -2, 1)$value
E_Y
0

# alternativně
E_Y = 3*E_X + 1
E_Y

# E(Y^2)
yy_fy = function(y){
    fy = fY.dens(y)
    return(y*y*fy)
} 
# integrujeme jen tam kde víme, že je f(y) nenulová
E_YY = integrate(yy_fy, -2, 1)$value
E_YY
1/2

# D(Y)
D_Y = E_YY - E_Y^2
D_Y
1/2

# alternativně
D_Y = 3^2*D_X
D_Y

# sigma(Y)
sqrt(D_Y)
sqrt(2)/2

# * Příklad 4. (není ze sbírky) ####
# Spočtěte $\omega$ takové, aby náhodná veličina X s hustotou pravděpodobnosti:
# $f(x)=\begin{cases}
# 0 & x < 0 \\
# 3e^{-3x}    & x \geq 0
# \end{cases}$ 
# byla s pravděpodobností 0.3 větší než $\omega$


F.dist = function(x){
    res = 1 - exp(-3*x)
    res[x < 0] = 0 # 0 pro x<=0
    return(res)
}

x = seq(from = -1, to = 3, by = 0.001) # body na ose x
FX = F.dist(x)
plot(x, FX, type='l')
lines(c(-1, 3),c(0.7, 0.7))

x[FX >= 0.7][1]

-1/3*log(0.3)