jilm/spm_math.html

## spm_math.html
<!DOCTYPE html>
<html>

    <head>
        <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
        <meta charset="UTF-8">
        <link rel="stylesheet" href="main.css">
        <script type="text/x-mathjax-config">
            MathJax.Hub.Config({
                TeX: { equationNumbers: { autoNumber: "all" } }
            });
        </script>
        <script type="text/javascript"
              src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_CHTML" async>
        </script>
        <title>Poznámky.</title>
    </head>

    <body>
        <div class="book">

            <p class="title">Poznámky k teorii spolehlivosti</p>

            <h1>Definice</h1>

            <p>
                Jestliže \(P(T \le t)\) je pravděpodobnost bezporuchového
                provozu nějakého zařízení po dobu \(t\), potom funkci:

                $$F(t) = P(T \le t),$$

                nazýváme <em>distribuční funkcí</em> a funkci:

                $$R(t) = P(T &gt; t) = 1 - F(t),$$

                nazýváme <em>survival funkce.</em>
            </p>

        <h1>Poissonovo rozdělení pravděpodobnosti</h1>

        <p>
            V případě, že střední doba mezi poruchami \(m\) je konstantní,
            pak počet poruch \(f\) za nějaké období \(T\) je náhodná veličina
            s poissonovým rozdělením pravděpodobnosti:

            $$P(X=f,\lambda) = p(f) = e^{-\lambda}{\lambda^f\over f!},\label{poisson}$$

            kde \(\lambda\) je parametr ve kterém se skrývá předpokládaná
            kumulovaná doba provozu \(T^*\) a střední doba mezi poruchami \(m\):

            $$\lambda = \frac{T^*}{m} = \frac{dT}{m}.$$
        </p>

        <img src="plot1.svg" >
        <div class="caption">Obrázek 1: pravděpodobnost z poissonova rozdělení</div>

        <p>
            Ve vzorci (\ref{poisson}) je součin dvou exponenciál
            a jednoho faktoriálu. Tyto funkce jsou známy
            pro poměrně rychlý růst hodnoty už pro malou hodnotu argumentu.
            Výpočet podle tohoto vzorce tak nemusí být numericky stabilní.
            Lépe je počítat logaritmus:

            $$p(k) = \exp[-\lambda + k \ln\lambda - \ln\Gamma(k+1)]$$

            Zejména proto, že funkce \(\ln\Gamma(x)\) bývá v řadě matematických
            knihoven implementována přímo.
        </p>

        <p>
            Například graf na obtrázku 1 byl nakreslen v gnuplotu za použití:
        </p>
        <code>
            poisson(k, lmbd) = exp(-lmbd+k*log(lmbd)-lgamma(k+1))
        </code>
        <p>
            Pro pravděpodobnost, že za období \(T\) spotřebujeme skladovou
            zásobu dílů \(s\) dostaneme:

            $$P(X \ge s) = F(s) = \sum_{k=0}^s e^{-\lambda}{\lambda^k\over k!},$$
        </p>

        <img src="plot2.svg">
        <!--gnuplot
            set term svg size 600, 400 font "Times New Roman,14"
            set output "plot2.svg"
            set xrange [0:60]
            set yrange [0:1]
            set xlabel "lambda [-]"
            set ylabel "Pravděpodobnost [-]"
            poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))
            cdfpoisson(s, lmbd) = 1.0 - sum [k=0:s] poisson(k, lmbd)
            plot for [s=1:4] cdfpoisson(s*10, x)
        -->
        </img>
        <div class="caption">Distribuční funkce pro různá s</div>

        <img src="plot3.svg">
        <!--gnuplot
set term svg size 650, 450 font "Times New Roman,14"
set output "plot3.svg"
set xrange [0:60]
set sample 61
set yrange [0:1]
set xlabel "lambda [-]"
set ylabel "Pravděpodobnost [-]"
poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))
cdfpoisson(s, lmbd) = 1.0 - sum [k=0:floor(s)] poisson(k, lmbd)
plot for [i=1:4] cdfpoisson(x, 10*i) with points
        -->
        </img>
        <div class="caption">Distribuční funkce pro různé lambda</div>

        <h1>Teorie her</h1>

        <p>
            Ke stanovení odhadu potřebného počtu náhradních dílů potřebujeme
            mimo jiné jakousi kritičnost, což je docela subjektivní údaj.
            Pokud ale známe nákupní cenu náhradního dílu \(c_1\) a navíc známe
            náklady \(c_\dagger\) tzv. alternativy, kterou je potřeba zaplatit
            v případě, že díly dojdou dříve, než je stanovené období \(T\),
            můžeme využít tzv. teorii her, viz [ORR].
        </p>

        <p>
            V tom případě se jedná o tzv. rozhodování při riziku. To znamená,
            že jde o nalezení optimální strategie ve hře dvou hráčů, kdy ten
            první je inteligentní a volí svou strategii, kterou označme
            \(u \in U\), ten druhý volí své strategie náhodně se známým
            rozdělením pravděpodobnosti \(P\) a je lhostejný ke své výhře.
        </p>

        <p>
            Při tomto typu hry je nutno pro každou strategii \(u\) prvního
            hráče vypočítat střední hodnotu výhry:

            $$E[J(u)] = \sum_{v \in V} J(u, v) p(v),$$

            kde \(J(u, v)\) je výhra prvního hráče při zvolených strategiích
            \(u\) a \(v\) a \(p(v)\) je pravděpodobnost strategie \(v\).
            Strategie prvního hráče je optimální, pokud je střední hodnota
            jeho výhry maximální:

            $$u^* = \arg \max_{u \in U} E[J(u)].$$
        </p>

        <p>
            V našem případě představují strategie \(u\) počty náhradních
            dílů \(s\), které se první hráč rozhodne nakoupit do skladu.
            Strategie \(v\) pak představují počty spotřebovaných dílů za
            období \(T\). Velikost výhry pak můžeme zvolit například takto:

            $$
                J(u, v) = \begin{cases}
                    -(u - v) c_1 & \textrm{pro}\quad 0 \leq v \leq u,\\
                    - c_\dagger & \textrm{pro}\quad u &lt; v.
                \end{cases}
            $$

            To znamená, že v případě, kdy náhradních dílů je dostatek,
            představuje výhra souhrnnou cenu náhradních dílů, které zůstaly
            na skladě po uplynutí doby \(T\), v opačném případě představuje
            cenu \(c_\dagger\)
        </p>

        <p>
            V případě, že  počet poruch pro dané období je dán poissonovým
            rozdělením pravděpodobnosti s parametrem \(\lambda\). Potom střední
            hodnota výhry pro každou ze strategií \(u\) je:

            $$\begin{align}
                E[u] &= \sum_{v \in V} J(u, v) p(v)=\\
                &= -c_1 \sum_{v = 0}^u (u-v) p(v) - c_\dagger R(v)=\\
                -c_\dagger E[J(u)] &= \frac{c_1}{c_\dagger}\sum_{v = 0}^u (u-v) p(v) + R(u)= \frac{c_1}{c_\dagger} E[J_1(u)] + R(u)\\
            \end{align}$$

            přitom je zřejmé, že:

            $$\arg\max E[\cdot] = \arg\min -c E[\cdot]\qquad\textrm{pro}\qquad c &gt; 0.$$
        </p>

        <p>
            Pro poissonovo rozdělení můžeme vzorec xxx dále rozepsat:

            $$\begin{align}
                \frac{c_1}{c_\dagger} E[J_1(u)] + R(u)
                = \frac{c_1}{c_\dagger} \sum_{v = 0}^u (u-v) \textrm{e}^{-\lambda} \frac{\lambda^v}{v!}
            + \left(1 - \sum_{v = 0}^u\textrm{e}^{-\lambda} \frac{\lambda^v}{v!}\right)
            \end{align}$$
        </p>

        <p>
            vzorec xxx je součet dvou členů. Funkce \(R(\cdot)\) je funkce
            neklesající a druhý člen je naopak funkce rostoucí a pro větší
            \(u\) je to funkce téměř lineární se směrnicí \(c_1/c_\dagger\).
            Pomyslný průsečík pak tvoří hledané optimum.
        </p>

        <img src="plot4.svg">


        <img src="plot5.svg">
        <!--gnuplot
set term svg size 650, 450 font "Times New Roman,14"
set output "plot4.svg"
set xrange [0:20]
set sample 21
set yrange [0:5]
set xlabel "Počet nd [ks]"
set ylabel "Relatuvní náklady [-]"
poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))
cdfpoisson(s, lmbd) = 1.0 - sum [k=0:floor(s)] poisson(k, lmbd)
strategy(cr, u, lmbd) = cdfpoisson(u, lmbd) + cr * sum [v=0:floor(u)] (u-v) * poisson(v, lmbd)
plot strategy(0.1, x, 5) with points, strategy(0.01, x, 5) with points, strategy(0.5, x, 5) with points
        -->
        </img>

        <img src="plot6.svg">
<!--
        set term svg size 650, 450 font "Times New Roman,14pt"
        set output "plot6.svg"
        #set term aqua
        set xrange [0:50]
        set sample 51
        set yrange [0:1]
        set xlabel "Počet nd [ks]"
        set ylabel "Relatuvní náklady [-]"

        # pravděpodobnost
        p_poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))

        f_poisson(k, lmbd) = sum [i=0:floor(k)] p_poisson(i, lmbd)

        r_poisson(k, lmbd) = 1.0 - f_poisson(k, lmbd)

        mean_cost1(u, cost_ratio, lmbd) = cost_ratio * sum [v=0:floor(u)] (u-v) * p_poisson(v, lmbd)

        plot r_poisson(x, 5) with points, mean_cost1(x, 0.01, 5) with points, mean_cost1(x, 0.05, 5) with points
-->
        <p>
            Pro vlastní nalezení minima té střední výhry by se mi i hodilo
            znát přírůstek \(E[u] - E[u+1] \). V tom případě bych mohl
            použít tzv. algoritmus binárního půlení s logaritmickou složitostí
            \(O(\log n)\).
        </p>

        <p>
            Pro distribuční funkci obecně platí:

            $$F(u+1) - F(u) = p(u+1)$$

            Druhý ze sčítanců můžeme po substituci \(u+1 \rightarrow u\) rozepsat
            jako:

            $$\frac{c_1}{c_\dagger}[(u+1-u-1)p(u+1)]
            + \frac{c_1}{c_\dagger}\sum_{v=0}^u(u-v)p(v)
            + \frac{c_1}{c_\dagger}\sum_{v=0}^u p(v) =  E_1[u] + F(u)$$
        </p>


<h2> Exponenciální rozdělení pravděpodobnosti </h2>

Hustota pravděpodobnosti exponenciálního rozdělení (probability density
function) se dá vyjádřit jako:

$$ f(x|\lambda)=\lambda e^{-\lambda x}; {\rm pro\; } x \geq 0 $$

### Konfidenční interval

$${2n \over {\lambda\chi^2_{{\alpha \over 2}, 2n}}} < {1 \over \lambda} < {2n \over {\lambda\chi^2_{1-{n\over 2}, 2n}}}$$

### Odhad parametrů exponenciálního rozdělení

Bodový maximum likelihood odhad parametru \lambda se vypočítá jako aritmetický průměr, viz [2].

$$ \hat{\lambda} = {f \over Td}, {\rm nebo} \hat{m} = {Td \over f} $$

<h2> Examples </h2>

<h3>s = f</h3>

člověk se občas setká s následujícím názorem: mám v provozu \(d\) dílů, které provozuji po dobu \(T\) a za tu dobu
jsem zaznamenal \(f\) poruch, chci-li provozovat stejné zařízení po stejnou dobu \(T\), musí mi stačit \(f\) dílů
na skladě. Proč jich kupovat více?

### Citlivost potřebného počtu nd na změnu parametru lambda

Už víme, že počet náhradních dílů, které je potřeba nakoupit na sklad je možno odhadnout z poissonova rozdělení
pravděpodobnosti. Jeho parametr lambda pak v sobě skrývá jak počty nasazených dílů, potřebnou dobu provozu i
spolehlivost dílu.

$$ \lambda = \frac{dT}{MTBF} $$

Ze vzorce pro poissonovo rozdělení není závislost mezi \lambda a potřebným počtem dílů ihned patrná. Někdy se
ale také používá zjednodušený vztah kde ta závislost je mnohem jasnější:

$$ s = \lambda + N^{-1}(P) * \sqrt(\lambda) $$

Závislost počtu dílů ve skladu s je lineárně závislá na \lambda, zejména pro větší \lamba.


<h1> Reference </h1>

[1] [Poisson distribution. (15 August 2018 05:20 UTC). In Wikipedia, The Free Encyclopedia. Retrieved  15 August 2018 20:09 UTC]( https://en.wikipedia.org/w/index.php?title=Poisson_distribution&oldid=854992626)

[ČSN IEC 60605-4] [Zkoušení bezporuchovosti zařízení -- Část 4: Statistické postupy pro exponenciální rozdělení -- Bodové odhady, konfidenční intervaly, předpovědní intervaly a toleranční intervaly. Vydal Český Normalizační Institut, Praha, 2002]

[ORR]

  <form action="/action_page.php" oninput="x.value=parseInt(d.value)*parseInt(f.value) / parseFloat(T.value)">
     Deployed: <input name="d" type="text" /> </br>
     Poruch: <input name="f" type="text" /> </br>
     Délka zkoušky:  <input name="T" type="text" /> </br>
      <output name="x" for="d f T"></output>
   </form>

<!--

# copyright 2019 Jiří Lidinský

# ---------------------------------------------
# Definice funkcí poissonova rozdělení

# pravděpodobnost, že nastane právě k událostí
# k >= 0, int
# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
# vrací real z intervalu <0, 1.0>
p_poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))

# distribuční funkce poissonova rozdělení
# k >= 0, int
# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
f_poisson(k, lmbd) = sum [i=0:floor(k)] p_poisson(i, lmbd)

# survival funkce R() = 1-F()
r_poisson(k, lmbd) = 1.0 - f_poisson(k, lmbd)

mean_cost1(u, cost_ratio, lmbd) = cost_ratio * sum [v=0:floor(u)] (u-v) * p_poisson(v, lmbd)

diff_mean_cost(u, cost_ratio, lmbd) = cost_ratio * f_poisson(u, lmbd) - p_poisson(u+1, lmbd)

arg_min_step(a, b, cost_ratio, lmbd) \
    = (b - a) <= 1 \
    ? diff_mean_cost(a, cost_ratio, lmbd) < 0 ? b : a \
    : diff_mean_cost(floor((a+b)/2), cost_ratio, lmbd) < 0 \
        ? arg_min_step(floor((a+b)/2), b, cost_ratio, lmbd) \
        : arg_min_step(a, floor((a+b)/2), cost_ratio, lmbd)

arg_min(cost_ratio, lmbd) = arg_min_step(0, floor(1.0/cost_ratio), cost_ratio, lmbd)

# ---------------------------------------------
# Společná nastavení

set term svg size 650, 450 font "Palatino,14pt"
set grid
set style line 5 linecolor rgbcolor "black"
main_color = "linecolor rgbcolor 'red'"
aux_color = "linecolor rgbcolor 'gray'"
main_points = "points linestyle 6 "
main_line = "lines"
aux_line = "lines dashtype 2"

# ---------------------------------------------
# Obrazek c. 1 - Poissonova Pravděpodobnost

xmax = 30       # max value on axis x
ymax = 0.2      # max value on y axis
aux1 = 7        # points for auxiliary lines
aux2 = 23
aux_offset = 0.01

set output "plot1.svg"
set xrange [0:xmax]
set yrange [0:ymax]
set sample xmax + 1
set xlabel "Počet událostí [-]"
set ylabel "Pravděpodobnost [-]"
set arrow from 8, p_poisson(8, 8) to 7, p_poisson(7, 7)
set label "menší" at aux1, p_poisson(aux1, aux1) + 2*aux_offset left \
    rotate by atan((p_poisson(aux1+1, aux1+1) - p_poisson(aux1, aux1))*xmax/1/ymax) *180/pi
set label "větší" at aux2, p_poisson(aux2, aux2) + 2*aux_offset right \
    rotate by atan((p_poisson(aux2, aux2) - p_poisson(aux2-1, aux2-1))*xmax/1/ymax) *180/pi
set arrow \
    from aux2-1, p_poisson(aux2-1, aux2-1) \
    to aux2, p_poisson(aux2, aux2)
unset grid    # remove unnecesary components of the plot
unset xtics
unset ytics
plot p_poisson(x, 15) with @main_points @main_color notitle, \
    [17:23] p_poisson(x, 20) with @aux_line @aux_color notitle, \
    [7:12] p_poisson(x, 10) with @aux_line @aux_color notitle, \
    [aux1:aux2] p_poisson(x, x) + aux_offset with @aux_line @aux_color notitle
# some cleenup
unset label
unset arrow

# Obrazek c. 2
set output "plot2.svg"
set xrange [0:60]
unset sample
set yrange [0:1]
set ytics
set xtics
set xlabel "lambda [-]"
set ylabel "Pravděpodobnost [-]"
plot f_poisson(35, x) with lines @main_color notitle,\
    [25:30] f_poisson(25, x) with @aux_line @aux_color notitle,\
    f_poisson(15, x) with @aux_line @aux_color notitle,\
    f_poisson(5, x) with @aux_line @aux_color notitle,\
    f_poisson(45, x) with @aux_line @aux_color notitle, \
    f_poisson(x, x)


# Obrazek c. 3
set output "plot3.svg"
set xrange [0:60]
set sample 61
set yrange [0:1]
set ytics
set xlabel "lambda [-]"
set ylabel "Pravděpodobnost [-]"
plot f_poisson(x, 35) with @main_points @main_color notitle, \
    f_poisson(x, 25) with @main_points @aux_color notitle, \
    f_poisson(x, 45) with @main_points @aux_color notitle

# Obrazek c. 4
set output "plot4.svg"
set xrange [0:20]
set sample 21
set yrange [0:1.2]
set xlabel "Počet nd [ks]"
set ylabel "Relatuvní náklady [-]"
plot r_poisson(x, 5) ls 7, \
    mean_cost1(x, 0.1, 5) ls 7, \
    [5:8] mean_cost1(x, 0.19, 5) ls 7, \
    [6:10] mean_cost1(x, 0.05, 5) ls 7, \
    [7:12] mean_cost1(x, 0.01, 5) ls 7, \
    mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points ls 6

# Obrazek c. 5
set output "plot5.svg"
set xrange [0:20]
set sample 21
set yrange [0:5]
set xlabel "Počet nd [ks]"
set ylabel "Relatuvní náklady [-]"
plot mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points, \
    mean_cost1(x, 0.01, 5) + r_poisson(x, 5) with points, \
    mean_cost1(x, 0.5, 5) + r_poisson(x, 5) with points
# copyright 2019 Jiří Lidinský

# ---------------------------------------------
# Definice funkcí poissonova rozdělení

# pravděpodobnost, že nastane právě k událostí
# k >= 0, int
# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
# vrací real z intervalu <0, 1.0>
p_poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))

# distribuční funkce poissonova rozdělení
# k >= 0, int
# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
f_poisson(k, lmbd) = sum [i=0:floor(k)] p_poisson(i, lmbd)

# survival funkce R() = 1-F()
r_poisson(k, lmbd) = 1.0 - f_poisson(k, lmbd)

mean_cost1(u, cost_ratio, lmbd) = cost_ratio * sum [v=0:floor(u)] (u-v) * p_poisson(v, lmbd)

diff_mean_cost(u, cost_ratio, lmbd) = cost_ratio * f_poisson(u, lmbd) - p_poisson(u+1, lmbd)

arg_min_step(a, b, cost_ratio, lmbd) \
    = (b - a) <= 1 \
    ? diff_mean_cost(a, cost_ratio, lmbd) < 0 ? b : a \
    : diff_mean_cost(floor((a+b)/2), cost_ratio, lmbd) < 0 \
        ? arg_min_step(floor((a+b)/2), b, cost_ratio, lmbd) \
        : arg_min_step(a, floor((a+b)/2), cost_ratio, lmbd)

arg_min(cost_ratio, lmbd) = arg_min_step(0, floor(1.0/cost_ratio), cost_ratio, lmbd)

# ---------------------------------------------
# Společná nastavení

set term svg size 650, 450 font "Palatino,14pt"
set grid
set style line 5 linecolor rgbcolor "black"
main_color = "linecolor rgbcolor 'red'"
aux_color = "linecolor rgbcolor 'gray'"
main_points = "points linestyle 6 "
main_line = "lines"
aux_line = "lines dashtype 2"

# ---------------------------------------------
# Obrazek c. 1 - Poissonova Pravděpodobnost

xmax = 30       # max value on axis x
ymax = 0.2      # max value on y axis
aux1 = 7        # points for auxiliary lines
aux2 = 23
aux_offset = 0.01

set output "plot1.svg"
set xrange [0:xmax]
set yrange [0:ymax]
set sample xmax + 1
set xlabel "Počet událostí [-]"
set ylabel "Pravděpodobnost [-]"
set arrow from 8, p_poisson(8, 8) to 7, p_poisson(7, 7)
set label "menší" at aux1, p_poisson(aux1, aux1) + 2*aux_offset left \
    rotate by atan((p_poisson(aux1+1, aux1+1) - p_poisson(aux1, aux1))*xmax/1/ymax) *180/pi
set label "větší" at aux2, p_poisson(aux2, aux2) + 2*aux_offset right \
    rotate by atan((p_poisson(aux2, aux2) - p_poisson(aux2-1, aux2-1))*xmax/1/ymax) *180/pi
set arrow \
    from aux2-1, p_poisson(aux2-1, aux2-1) \
    to aux2, p_poisson(aux2, aux2)
unset grid    # remove unnecesary components of the plot
unset xtics
unset ytics
plot p_poisson(x, 15) with @main_points @main_color notitle, \
    [17:23] p_poisson(x, 20) with @aux_line @aux_color notitle, \
    [7:12] p_poisson(x, 10) with @aux_line @aux_color notitle, \
    [aux1:aux2] p_poisson(x, x) + aux_offset with @aux_line @aux_color notitle
# some cleenup
unset label
unset arrow

# Obrazek c. 2
set output "plot2.svg"
set xrange [0:60]
unset sample
set yrange [0:1]
set ytics
set xtics
set xlabel "lambda [-]"
set ylabel "Pravděpodobnost [-]"
plot f_poisson(35, x) with lines @main_color notitle,\
    [25:30] f_poisson(25, x) with @aux_line @aux_color notitle,\
    f_poisson(15, x) with @aux_line @aux_color notitle,\
    f_poisson(5, x) with @aux_line @aux_color notitle,\
    f_poisson(45, x) with @aux_line @aux_color notitle, \
    f_poisson(x, x)


# Obrazek c. 3
set output "plot3.svg"
set xrange [0:60]
set sample 61
set yrange [0:1]
set ytics
set xlabel "lambda [-]"
set ylabel "Pravděpodobnost [-]"
plot f_poisson(x, 35) with @main_points @main_color notitle, \
    f_poisson(x, 25) with @main_points @aux_color notitle, \
    f_poisson(x, 45) with @main_points @aux_color notitle

# Obrazek c. 4
set output "plot4.svg"
set xrange [0:20]
set sample 21
set yrange [0:1.2]
set xlabel "Počet nd [ks]"
set ylabel "Relatuvní náklady [-]"
plot r_poisson(x, 5) ls 7, \
    mean_cost1(x, 0.1, 5) ls 7, \
    [5:8] mean_cost1(x, 0.19, 5) ls 7, \
    [6:10] mean_cost1(x, 0.05, 5) ls 7, \
    [7:12] mean_cost1(x, 0.01, 5) ls 7, \
    mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points ls 6

# Obrazek c. 5
set output "plot5.svg"
set xrange [0:20]
set sample 21
set yrange [0:5]
set xlabel "Počet nd [ks]"
set ylabel "Relatuvní náklady [-]"
plot mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points, \
    mean_cost1(x, 0.01, 5) + r_poisson(x, 5) with points, \
    mean_cost1(x, 0.5, 5) + r_poisson(x, 5) with points

-->

  </body>
</html>
	<!DOCTYPE html>
	<html>

	<head>
	<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
	<meta charset="UTF-8">
	<link rel="stylesheet" href="main.css">
	<script type="text/x-mathjax-config">
	MathJax.Hub.Config({
	TeX: { equationNumbers: { autoNumber: "all" } }
	});
	</script>
	<script type="text/javascript"
	src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_CHTML" async>
	</script>
	<title>Poznámky.</title>
	</head>

	<body>
	<div class="book">

	<p class="title">Poznámky k teorii spolehlivosti</p>

	<h1>Definice</h1>

	<p>
	Jestliže \(P(T \le t)\) je pravděpodobnost bezporuchového
	provozu nějakého zařízení po dobu \(t\), potom funkci:

	$$F(t) = P(T \le t),$$

	nazýváme <em>distribuční funkcí</em> a funkci:

	$$R(t) = P(T > t) = 1 - F(t),$$

	nazýváme <em>survival funkce.</em>
	</p>

	<h1>Poissonovo rozdělení pravděpodobnosti</h1>

	<p>
	V případě, že střední doba mezi poruchami \(m\) je konstantní,
	pak počet poruch \(f\) za nějaké období \(T\) je náhodná veličina
	s poissonovým rozdělením pravděpodobnosti:

	$$P(X=f,\lambda) = p(f) = e^{-\lambda}{\lambda^f\over f!},\label{poisson}$$

	kde \(\lambda\) je parametr ve kterém se skrývá předpokládaná
	kumulovaná doba provozu \(T^*\) a střední doba mezi poruchami \(m\):

	$$\lambda = \frac{T^*}{m} = \frac{dT}{m}.$$
	</p>

	<img src="plot1.svg" >
	<div class="caption">Obrázek 1: pravděpodobnost z poissonova rozdělení</div>

	<p>
	Ve vzorci (\ref{poisson}) je součin dvou exponenciál
	a jednoho faktoriálu. Tyto funkce jsou známy
	pro poměrně rychlý růst hodnoty už pro malou hodnotu argumentu.
	Výpočet podle tohoto vzorce tak nemusí být numericky stabilní.
	Lépe je počítat logaritmus:

	$$p(k) = \exp[-\lambda + k \ln\lambda - \ln\Gamma(k+1)]$$

	Zejména proto, že funkce \(\ln\Gamma(x)\) bývá v řadě matematických
	knihoven implementována přímo.
	</p>

	<p>
	Například graf na obtrázku 1 byl nakreslen v gnuplotu za použití:
	</p>
	<code>
	poisson(k, lmbd) = exp(-lmbd+k*log(lmbd)-lgamma(k+1))
	</code>
	<p>
	Pro pravděpodobnost, že za období \(T\) spotřebujeme skladovou
	zásobu dílů \(s\) dostaneme:

	$$P(X \ge s) = F(s) = \sum_{k=0}^s e^{-\lambda}{\lambda^k\over k!},$$
	</p>

	<img src="plot2.svg">
	<!--gnuplot
	set term svg size 600, 400 font "Times New Roman,14"
	set output "plot2.svg"
	set xrange [0:60]
	set yrange [0:1]
	set xlabel "lambda [-]"
	set ylabel "Pravděpodobnost [-]"
	poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))
	cdfpoisson(s, lmbd) = 1.0 - sum [k=0:s] poisson(k, lmbd)
	plot for [s=1:4] cdfpoisson(s*10, x)
	-->
	</img>
	<div class="caption">Distribuční funkce pro různá s</div>

	<img src="plot3.svg">
	<!--gnuplot
	set term svg size 650, 450 font "Times New Roman,14"
	set output "plot3.svg"
	set xrange [0:60]
	set sample 61
	set yrange [0:1]
	set xlabel "lambda [-]"
	set ylabel "Pravděpodobnost [-]"
	poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))
	cdfpoisson(s, lmbd) = 1.0 - sum [k=0:floor(s)] poisson(k, lmbd)
	plot for [i=1:4] cdfpoisson(x, 10*i) with points
	-->
	</img>
	<div class="caption">Distribuční funkce pro různé lambda</div>

	<h1>Teorie her</h1>

	<p>
	Ke stanovení odhadu potřebného počtu náhradních dílů potřebujeme
	mimo jiné jakousi kritičnost, což je docela subjektivní údaj.
	Pokud ale známe nákupní cenu náhradního dílu \(c_1\) a navíc známe
	náklady \(c_\dagger\) tzv. alternativy, kterou je potřeba zaplatit
	v případě, že díly dojdou dříve, než je stanovené období \(T\),
	můžeme využít tzv. teorii her, viz [ORR].
	</p>

	<p>
	V tom případě se jedná o tzv. rozhodování při riziku. To znamená,
	že jde o nalezení optimální strategie ve hře dvou hráčů, kdy ten
	první je inteligentní a volí svou strategii, kterou označme
	\(u \in U\), ten druhý volí své strategie náhodně se známým
	rozdělením pravděpodobnosti \(P\) a je lhostejný ke své výhře.
	</p>

	<p>
	Při tomto typu hry je nutno pro každou strategii \(u\) prvního
	hráče vypočítat střední hodnotu výhry:

	$$E[J(u)] = \sum_{v \in V} J(u, v) p(v),$$

	kde \(J(u, v)\) je výhra prvního hráče při zvolených strategiích
	\(u\) a \(v\) a \(p(v)\) je pravděpodobnost strategie \(v\).
	Strategie prvního hráče je optimální, pokud je střední hodnota
	jeho výhry maximální:

	$$u^* = \arg \max_{u \in U} E[J(u)].$$
	</p>

	<p>
	V našem případě představují strategie \(u\) počty náhradních
	dílů \(s\), které se první hráč rozhodne nakoupit do skladu.
	Strategie \(v\) pak představují počty spotřebovaných dílů za
	období \(T\). Velikost výhry pak můžeme zvolit například takto:

	$$
	J(u, v) = \begin{cases}
	-(u - v) c_1 & \textrm{pro}\quad 0 \leq v \leq u,\\
	- c_\dagger & \textrm{pro}\quad u < v.
	\end{cases}
	$$

	To znamená, že v případě, kdy náhradních dílů je dostatek,
	představuje výhra souhrnnou cenu náhradních dílů, které zůstaly
	na skladě po uplynutí doby \(T\), v opačném případě představuje
	cenu \(c_\dagger\)
	</p>

	<p>
	V případě, že počet poruch pro dané období je dán poissonovým
	rozdělením pravděpodobnosti s parametrem \(\lambda\). Potom střední
	hodnota výhry pro každou ze strategií \(u\) je:

	$$\begin{align}
	E[u] &= \sum_{v \in V} J(u, v) p(v)=\\
	&= -c_1 \sum_{v = 0}^u (u-v) p(v) - c_\dagger R(v)=\\
	-c_\dagger E[J(u)] &= \frac{c_1}{c_\dagger}\sum_{v = 0}^u (u-v) p(v) + R(u)= \frac{c_1}{c_\dagger} E[J_1(u)] + R(u)\\
	\end{align}$$

	přitom je zřejmé, že:

	$$\arg\max E[\cdot] = \arg\min -c E[\cdot]\qquad\textrm{pro}\qquad c > 0.$$
	</p>

	<p>
	Pro poissonovo rozdělení můžeme vzorec xxx dále rozepsat:

	$$\begin{align}
	\frac{c_1}{c_\dagger} E[J_1(u)] + R(u)
	= \frac{c_1}{c_\dagger} \sum_{v = 0}^u (u-v) \textrm{e}^{-\lambda} \frac{\lambda^v}{v!}
	+ \left(1 - \sum_{v = 0}^u\textrm{e}^{-\lambda} \frac{\lambda^v}{v!}\right)
	\end{align}$$
	</p>

	<p>
	vzorec xxx je součet dvou členů. Funkce \(R(\cdot)\) je funkce
	neklesající a druhý člen je naopak funkce rostoucí a pro větší
	\(u\) je to funkce téměř lineární se směrnicí \(c_1/c_\dagger\).
	Pomyslný průsečík pak tvoří hledané optimum.
	</p>

	<img src="plot4.svg">


	<img src="plot5.svg">
	<!--gnuplot
	set term svg size 650, 450 font "Times New Roman,14"
	set output "plot4.svg"
	set xrange [0:20]
	set sample 21
	set yrange [0:5]
	set xlabel "Počet nd [ks]"
	set ylabel "Relatuvní náklady [-]"
	poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))
	cdfpoisson(s, lmbd) = 1.0 - sum [k=0:floor(s)] poisson(k, lmbd)
	strategy(cr, u, lmbd) = cdfpoisson(u, lmbd) + cr * sum [v=0:floor(u)] (u-v) * poisson(v, lmbd)
	plot strategy(0.1, x, 5) with points, strategy(0.01, x, 5) with points, strategy(0.5, x, 5) with points
	-->
	</img>

	<img src="plot6.svg">
	<!--
	set term svg size 650, 450 font "Times New Roman,14pt"
	set output "plot6.svg"
	#set term aqua
	set xrange [0:50]
	set sample 51
	set yrange [0:1]
	set xlabel "Počet nd [ks]"
	set ylabel "Relatuvní náklady [-]"

	# pravděpodobnost
	p_poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))

	f_poisson(k, lmbd) = sum [i=0:floor(k)] p_poisson(i, lmbd)

	r_poisson(k, lmbd) = 1.0 - f_poisson(k, lmbd)

	mean_cost1(u, cost_ratio, lmbd) = cost_ratio * sum [v=0:floor(u)] (u-v) * p_poisson(v, lmbd)

	plot r_poisson(x, 5) with points, mean_cost1(x, 0.01, 5) with points, mean_cost1(x, 0.05, 5) with points
	-->
	<p>
	Pro vlastní nalezení minima té střední výhry by se mi i hodilo
	znát přírůstek \(E[u] - E[u+1] \). V tom případě bych mohl
	použít tzv. algoritmus binárního půlení s logaritmickou složitostí
	\(O(\log n)\).
	</p>

	<p>
	Pro distribuční funkci obecně platí:

	$$F(u+1) - F(u) = p(u+1)$$

	Druhý ze sčítanců můžeme po substituci \(u+1 \rightarrow u\) rozepsat
	jako:

	$$\frac{c_1}{c_\dagger}[(u+1-u-1)p(u+1)]
	+ \frac{c_1}{c_\dagger}\sum_{v=0}^u(u-v)p(v)
	+ \frac{c_1}{c_\dagger}\sum_{v=0}^u p(v) = E_1[u] + F(u)$$
	</p>


	<h2> Exponenciální rozdělení pravděpodobnosti </h2>

	Hustota pravděpodobnosti exponenciálního rozdělení (probability density
	function) se dá vyjádřit jako:

	$$ f(x\|\lambda)=\lambda e^{-\lambda x}; {\rm pro\; } x \geq 0 $$

	### Konfidenční interval

	$${2n \over {\lambda\chi^2_{{\alpha \over 2}, 2n}}} < {1 \over \lambda} < {2n \over {\lambda\chi^2_{1-{n\over 2}, 2n}}}$$

	### Odhad parametrů exponenciálního rozdělení

	Bodový maximum likelihood odhad parametru \lambda se vypočítá jako aritmetický průměr, viz [2].

	$$ \hat{\lambda} = {f \over Td}, {\rm nebo} \hat{m} = {Td \over f} $$

	<h2> Examples </h2>

	<h3>s = f</h3>

	člověk se občas setká s následujícím názorem: mám v provozu \(d\) dílů, které provozuji po dobu \(T\) a za tu dobu
	jsem zaznamenal \(f\) poruch, chci-li provozovat stejné zařízení po stejnou dobu \(T\), musí mi stačit \(f\) dílů
	na skladě. Proč jich kupovat více?

	### Citlivost potřebného počtu nd na změnu parametru lambda

	Už víme, že počet náhradních dílů, které je potřeba nakoupit na sklad je možno odhadnout z poissonova rozdělení
	pravděpodobnosti. Jeho parametr lambda pak v sobě skrývá jak počty nasazených dílů, potřebnou dobu provozu i
	spolehlivost dílu.

	$$ \lambda = \frac{dT}{MTBF} $$

	Ze vzorce pro poissonovo rozdělení není závislost mezi \lambda a potřebným počtem dílů ihned patrná. Někdy se
	ale také používá zjednodušený vztah kde ta závislost je mnohem jasnější:

	$$ s = \lambda + N^{-1}(P) * \sqrt(\lambda) $$

	Závislost počtu dílů ve skladu s je lineárně závislá na \lambda, zejména pro větší \lamba.


	<h1> Reference </h1>

	[1] [Poisson distribution. (15 August 2018 05:20 UTC). In Wikipedia, The Free Encyclopedia. Retrieved 15 August 2018 20:09 UTC]( https://en.wikipedia.org/w/index.php?title=Poisson_distribution&oldid=854992626)

	[ČSN IEC 60605-4] [Zkoušení bezporuchovosti zařízení -- Část 4: Statistické postupy pro exponenciální rozdělení -- Bodové odhady, konfidenční intervaly, předpovědní intervaly a toleranční intervaly. Vydal Český Normalizační Institut, Praha, 2002]

	[ORR]

	<form action="/action_page.php" oninput="x.value=parseInt(d.value)*parseInt(f.value) / parseFloat(T.value)">
	Deployed: <input name="d" type="text" /> </br>
	Poruch: <input name="f" type="text" /> </br>
	Délka zkoušky: <input name="T" type="text" /> </br>
	<output name="x" for="d f T"></output>
	</form>

	<!--

	# copyright 2019 Jiří Lidinský

	# ---------------------------------------------
	# Definice funkcí poissonova rozdělení

	# pravděpodobnost, že nastane právě k událostí
	# k >= 0, int
	# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
	# vrací real z intervalu <0, 1.0>
	p_poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))

	# distribuční funkce poissonova rozdělení
	# k >= 0, int
	# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
	f_poisson(k, lmbd) = sum [i=0:floor(k)] p_poisson(i, lmbd)

	# survival funkce R() = 1-F()
	r_poisson(k, lmbd) = 1.0 - f_poisson(k, lmbd)

	mean_cost1(u, cost_ratio, lmbd) = cost_ratio * sum [v=0:floor(u)] (u-v) * p_poisson(v, lmbd)

	diff_mean_cost(u, cost_ratio, lmbd) = cost_ratio * f_poisson(u, lmbd) - p_poisson(u+1, lmbd)

	arg_min_step(a, b, cost_ratio, lmbd) \
	= (b - a) <= 1 \
	? diff_mean_cost(a, cost_ratio, lmbd) < 0 ? b : a \
	: diff_mean_cost(floor((a+b)/2), cost_ratio, lmbd) < 0 \
	? arg_min_step(floor((a+b)/2), b, cost_ratio, lmbd) \
	: arg_min_step(a, floor((a+b)/2), cost_ratio, lmbd)

	arg_min(cost_ratio, lmbd) = arg_min_step(0, floor(1.0/cost_ratio), cost_ratio, lmbd)

	# ---------------------------------------------
	# Společná nastavení

	set term svg size 650, 450 font "Palatino,14pt"
	set grid
	set style line 5 linecolor rgbcolor "black"
	main_color = "linecolor rgbcolor 'red'"
	aux_color = "linecolor rgbcolor 'gray'"
	main_points = "points linestyle 6 "
	main_line = "lines"
	aux_line = "lines dashtype 2"

	# ---------------------------------------------
	# Obrazek c. 1 - Poissonova Pravděpodobnost

	xmax = 30 # max value on axis x
	ymax = 0.2 # max value on y axis
	aux1 = 7 # points for auxiliary lines
	aux2 = 23
	aux_offset = 0.01

	set output "plot1.svg"
	set xrange [0:xmax]
	set yrange [0:ymax]
	set sample xmax + 1
	set xlabel "Počet událostí [-]"
	set ylabel "Pravděpodobnost [-]"
	set arrow from 8, p_poisson(8, 8) to 7, p_poisson(7, 7)
	set label "menší" at aux1, p_poisson(aux1, aux1) + 2*aux_offset left \
	rotate by atan((p_poisson(aux1+1, aux1+1) - p_poisson(aux1, aux1))xmax/1/ymax) 180/pi
	set label "větší" at aux2, p_poisson(aux2, aux2) + 2*aux_offset right \
	rotate by atan((p_poisson(aux2, aux2) - p_poisson(aux2-1, aux2-1))xmax/1/ymax) 180/pi
	set arrow \
	from aux2-1, p_poisson(aux2-1, aux2-1) \
	to aux2, p_poisson(aux2, aux2)
	unset grid # remove unnecesary components of the plot
	unset xtics
	unset ytics
	plot p_poisson(x, 15) with @main_points @main_color notitle, \
	[17:23] p_poisson(x, 20) with @aux_line @aux_color notitle, \
	[7:12] p_poisson(x, 10) with @aux_line @aux_color notitle, \
	[aux1:aux2] p_poisson(x, x) + aux_offset with @aux_line @aux_color notitle
	# some cleenup
	unset label
	unset arrow

	# Obrazek c. 2
	set output "plot2.svg"
	set xrange [0:60]
	unset sample
	set yrange [0:1]
	set ytics
	set xtics
	set xlabel "lambda [-]"
	set ylabel "Pravděpodobnost [-]"
	plot f_poisson(35, x) with lines @main_color notitle,\
	[25:30] f_poisson(25, x) with @aux_line @aux_color notitle,\
	f_poisson(15, x) with @aux_line @aux_color notitle,\
	f_poisson(5, x) with @aux_line @aux_color notitle,\
	f_poisson(45, x) with @aux_line @aux_color notitle, \
	f_poisson(x, x)


	# Obrazek c. 3
	set output "plot3.svg"
	set xrange [0:60]
	set sample 61
	set yrange [0:1]
	set ytics
	set xlabel "lambda [-]"
	set ylabel "Pravděpodobnost [-]"
	plot f_poisson(x, 35) with @main_points @main_color notitle, \
	f_poisson(x, 25) with @main_points @aux_color notitle, \
	f_poisson(x, 45) with @main_points @aux_color notitle

	# Obrazek c. 4
	set output "plot4.svg"
	set xrange [0:20]
	set sample 21
	set yrange [0:1.2]
	set xlabel "Počet nd [ks]"
	set ylabel "Relatuvní náklady [-]"
	plot r_poisson(x, 5) ls 7, \
	mean_cost1(x, 0.1, 5) ls 7, \
	[5:8] mean_cost1(x, 0.19, 5) ls 7, \
	[6:10] mean_cost1(x, 0.05, 5) ls 7, \
	[7:12] mean_cost1(x, 0.01, 5) ls 7, \
	mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points ls 6

	# Obrazek c. 5
	set output "plot5.svg"
	set xrange [0:20]
	set sample 21
	set yrange [0:5]
	set xlabel "Počet nd [ks]"
	set ylabel "Relatuvní náklady [-]"
	plot mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points, \
	mean_cost1(x, 0.01, 5) + r_poisson(x, 5) with points, \
	mean_cost1(x, 0.5, 5) + r_poisson(x, 5) with points
	# copyright 2019 Jiří Lidinský

	# ---------------------------------------------
	# Definice funkcí poissonova rozdělení

	# pravděpodobnost, že nastane právě k událostí
	# k >= 0, int
	# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
	# vrací real z intervalu <0, 1.0>
	p_poisson(k, lmbd) = exp(-lmbd + k * log(lmbd) - lgamma(k+1))

	# distribuční funkce poissonova rozdělení
	# k >= 0, int
	# lmbd > 0, real, parametr poissonova rozdělení, bezrozměrný
	f_poisson(k, lmbd) = sum [i=0:floor(k)] p_poisson(i, lmbd)

	# survival funkce R() = 1-F()
	r_poisson(k, lmbd) = 1.0 - f_poisson(k, lmbd)

	mean_cost1(u, cost_ratio, lmbd) = cost_ratio * sum [v=0:floor(u)] (u-v) * p_poisson(v, lmbd)

	diff_mean_cost(u, cost_ratio, lmbd) = cost_ratio * f_poisson(u, lmbd) - p_poisson(u+1, lmbd)

	arg_min_step(a, b, cost_ratio, lmbd) \
	= (b - a) <= 1 \
	? diff_mean_cost(a, cost_ratio, lmbd) < 0 ? b : a \
	: diff_mean_cost(floor((a+b)/2), cost_ratio, lmbd) < 0 \
	? arg_min_step(floor((a+b)/2), b, cost_ratio, lmbd) \
	: arg_min_step(a, floor((a+b)/2), cost_ratio, lmbd)

	arg_min(cost_ratio, lmbd) = arg_min_step(0, floor(1.0/cost_ratio), cost_ratio, lmbd)

	# ---------------------------------------------
	# Společná nastavení

	set term svg size 650, 450 font "Palatino,14pt"
	set grid
	set style line 5 linecolor rgbcolor "black"
	main_color = "linecolor rgbcolor 'red'"
	aux_color = "linecolor rgbcolor 'gray'"
	main_points = "points linestyle 6 "
	main_line = "lines"
	aux_line = "lines dashtype 2"

	# ---------------------------------------------
	# Obrazek c. 1 - Poissonova Pravděpodobnost

	xmax = 30 # max value on axis x
	ymax = 0.2 # max value on y axis
	aux1 = 7 # points for auxiliary lines
	aux2 = 23
	aux_offset = 0.01

	set output "plot1.svg"
	set xrange [0:xmax]
	set yrange [0:ymax]
	set sample xmax + 1
	set xlabel "Počet událostí [-]"
	set ylabel "Pravděpodobnost [-]"
	set arrow from 8, p_poisson(8, 8) to 7, p_poisson(7, 7)
	set label "menší" at aux1, p_poisson(aux1, aux1) + 2*aux_offset left \
	rotate by atan((p_poisson(aux1+1, aux1+1) - p_poisson(aux1, aux1))xmax/1/ymax) 180/pi
	set label "větší" at aux2, p_poisson(aux2, aux2) + 2*aux_offset right \
	rotate by atan((p_poisson(aux2, aux2) - p_poisson(aux2-1, aux2-1))xmax/1/ymax) 180/pi
	set arrow \
	from aux2-1, p_poisson(aux2-1, aux2-1) \
	to aux2, p_poisson(aux2, aux2)
	unset grid # remove unnecesary components of the plot
	unset xtics
	unset ytics
	plot p_poisson(x, 15) with @main_points @main_color notitle, \
	[17:23] p_poisson(x, 20) with @aux_line @aux_color notitle, \
	[7:12] p_poisson(x, 10) with @aux_line @aux_color notitle, \
	[aux1:aux2] p_poisson(x, x) + aux_offset with @aux_line @aux_color notitle
	# some cleenup
	unset label
	unset arrow

	# Obrazek c. 2
	set output "plot2.svg"
	set xrange [0:60]
	unset sample
	set yrange [0:1]
	set ytics
	set xtics
	set xlabel "lambda [-]"
	set ylabel "Pravděpodobnost [-]"
	plot f_poisson(35, x) with lines @main_color notitle,\
	[25:30] f_poisson(25, x) with @aux_line @aux_color notitle,\
	f_poisson(15, x) with @aux_line @aux_color notitle,\
	f_poisson(5, x) with @aux_line @aux_color notitle,\
	f_poisson(45, x) with @aux_line @aux_color notitle, \
	f_poisson(x, x)


	# Obrazek c. 3
	set output "plot3.svg"
	set xrange [0:60]
	set sample 61
	set yrange [0:1]
	set ytics
	set xlabel "lambda [-]"
	set ylabel "Pravděpodobnost [-]"
	plot f_poisson(x, 35) with @main_points @main_color notitle, \
	f_poisson(x, 25) with @main_points @aux_color notitle, \
	f_poisson(x, 45) with @main_points @aux_color notitle

	# Obrazek c. 4
	set output "plot4.svg"
	set xrange [0:20]
	set sample 21
	set yrange [0:1.2]
	set xlabel "Počet nd [ks]"
	set ylabel "Relatuvní náklady [-]"
	plot r_poisson(x, 5) ls 7, \
	mean_cost1(x, 0.1, 5) ls 7, \
	[5:8] mean_cost1(x, 0.19, 5) ls 7, \
	[6:10] mean_cost1(x, 0.05, 5) ls 7, \
	[7:12] mean_cost1(x, 0.01, 5) ls 7, \
	mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points ls 6

	# Obrazek c. 5
	set output "plot5.svg"
	set xrange [0:20]
	set sample 21
	set yrange [0:5]
	set xlabel "Počet nd [ks]"
	set ylabel "Relatuvní náklady [-]"
	plot mean_cost1(x, 0.1, 5) + r_poisson(x, 5) with points, \
	mean_cost1(x, 0.01, 5) + r_poisson(x, 5) with points, \
	mean_cost1(x, 0.5, 5) + r_poisson(x, 5) with points

	-->

	</body>
	</html>