Rozkład sumy trzech niezależnych zmiennych losowych

[Ravy]

Hej,

czy ktoś mógłby pomóc w następującym problemie:

\(\displaystyle{ X, Y, Z \sim \mathrm{rand}(a)}\) niezależne, gdzie \(\displaystyle{ a}\) jest stałą, a \(\displaystyle{ \mathrm{rand}(a)}\) oznacza rozkład jednostajny dyskretny na przedziale \(\displaystyle{ \{ 0, 1, \ldots, a-1 \}.}\)

Interesuje mnie rozkład zmiennej \(\displaystyle{ SUMA = X + Y + Z.}\)

W szczególności \(\displaystyle{ P( SUMA > C ).}\)

Co w przypadku gdy \(\displaystyle{ C}\) jest zmienną losową o rozkładzie \(\displaystyle{ b ( \mathrm{const} ) \cdot \mathrm{rand}(x)}\) (\(\displaystyle{ x}\) const).

Tzn wtedy by mnie interesowała \(\displaystyle{ SUMA2 = X + Y + Z - b \cdot U}\)

\(\displaystyle{ P(SUMA2 > 0)}\)

Ktoś ma pomysł jak się do czegoś takiego zabrać?

[Zordon]

Ale w jakiej postaci chcesz uzyskać ten rozkład? Chcesz jakiś wzór, chcesz algorytm na obliczanie tego na komputerze, chcesz mieć wzór przybliżony?

[Ravy]

Potrzebny mi wzór lub niekosztowny algorytm obliczający.

Z uwagi na koszt, symulacje odpadają, musi być to mniej lub bardziej analityczne.
Przybliżony wzór, z niewielkim błędem (błąd musiałby być policzalny) też by dał radę.

Konkretniej, jak mam SUME2, to potrzebuję:

float fun1(a, b, X) \(\displaystyle{ = \mathcal{P} \left( \frac{ \mathrm{rand}(a) + \mathrm{rand}(a) + \mathrm{rand}(a)}{3} > X \cdot \mathrm{rand}(b) \right)}\)

float fun2(a, b, X) \(\displaystyle{ = \mathbb{E} \left( \frac{\mathrm{rand}(a) + \mathrm{rand}(a) + \mathrm{rand}(a)}{3} - X \cdot \mathrm{rand}(b) \right)}\)

Z grubsza coś takiego, z niewielkimi błędami względnymi.

Do tej pory skleciłem coś takiego, ale nie spełnia założeń jeżeli chodzi o wydajność
(dla zmniejszenia kosztu próbowałem trzymać prawdopodobieństwa w tablicy trójwymiarowej)

Kod: Zaznacz cały

float policz_szanse_sumy(int n, int k)
{
    if(k < 0 || k > 3*n || n < 0)
        return 0.0;

    if(k <= n)
        return (itof((k+1)*(k+2)) / 2.0) / pow(itof(n+1), 3.0);
    else if(k < 2*n)
        return itof(6*k*n + 3*n + 2 - 3*n*n - 2*k*k) / 2.0
            / pow(itof(n+1), 3.0);

    return itof((3*n + 1 - k)*(3*n + 2 -k)) / 2.0 / pow(itof(n+1), 3.0);
}

float dystrybuanta(int x, int y)
{
    mapping dystrybuanta;
    float res;

    if(y > x)
        y = x;

    if(y < 0)
        y = -1;

    return itof(y+1) / itof(x+1);
}

float policz_prawdopodobienstwo(int n, int x, int const)
{
    mapping map;
    int a = 3 * n + 1;
    float ret = 0.0;

    if(const <= -200)
        return 0.0;

    if(const > 200)
        return 100.0;


    if(mappingp(map = prawdopodobienstwa[n]))
    {
        if(mappingp(map = map[x]))
        {
            if(ret = map[const])
                return (ret == -1.0 ? 0.0 : ret);
        }
        else
            prawdopodobienstwa[n][x] = ([]);
    }
    else
        prawdopodobienstwa[n] = ([x : ([])]);

    while(a--)
        ret += policz_szanse_sumy(n, a)
            * dystrybuanta(x, a * (200 + const) / (200 - const));

    prawdopodobienstwa[n][x][const] = ret*100.0;

    return ret*100.0;
}

[Dasio11]

Nie próbowałem zrozumieć twojego kodu, ale być może poniższa funkcja oblicza to, czego potrzebujesz:

Kod: Zaznacz cały

#include <vector>

float f( unsigned int a, unsigned int b, float t )
{
    const unsigned int M = 1 + 3*(a-1); // Rozmiar V[3].

    std::vector< int > V[4];
    std::vector< int > S( M, 0 );

    /*
    V[0] będzie zawierać rozkład zmiennej 0    , tzn. V[0][j] =       P( 0     = j )
    V[1] będzie zawierać rozkład zmiennej X    , tzn. V[1][j] = a   * P( X     = j )
    V[2] będzie zawierać rozkład zmiennej X+Y  , tzn. V[2][j] = a^2 * P( X+Y   = j )
    V[3] będzie zawierać rozkład zmiennej X+Y+Z, tzn. V[3][j] = a^3 * P( X+Y+Z = j )
    Liczb tych nie dzielimy przez potęgi a, żeby jak najdłużej pracować na liczbach całkowitych i nie tracić precyzji na działaniach.
    */
    V[0].push_back( 1 );

    for( int i = 0; i < 3; ++i )
    {
        const unsigned int K = 1 +  i   *(a-1); // Rozmiar wektora V[i].
        const unsigned int L = 1 + (i+1)*(a-1); // Rozmiar wektora V[i+1].
        V[i+1] = std::vector< int >( L, 0 );

        /*
        Obliczamy sumy prefiksowe dla wektora V[i], czyli
            S[j] = V[i][0] + ... V[i][j].
        Dzięki temu będzie można w czasie stałym obliczać wartości V[i+1][j], zamiast wzoru
            V[i+1][j] = V[i][j] + V[i][j-1] + ... + V[i][j-(a-1)]
        stosując wzór
            V[i+1][j] = S[j] - S[j-a].
        */
        S[0] = V[i][0];
        for( unsigned int j = 1; j < K; ++j ) S[j] = S[j-1] + V[i][j];
        for( unsigned int j = K; j < L; ++j ) S[j] = S[j-1];

        for( unsigned int j = 0; j < a; ++j ) V[i+1][j] = S[j];
        for( unsigned int j = a; j < L; ++j ) V[i+1][j] = S[j] - S[j-a];
    }

    /*
    Niech S = X + Y + Z oraz niech W będzie zmienną losową o rozkładzie rand( b ). Wtedy
        P( (X+Y+Z)/3 > t*W ) = P( S > 3*t*W ) = P( S > 0 ) * P( W = 0 ) + P( S > 3*t ) * P( W = 1 ) + ... + P( S > 3*t*(b-1) ) * P( W = b-1 )
                             = [ P( S > 0 ) + P( S > 3*t ) + ... + P( S > 3*t*(b-1) ) ]/b.
    Teraz więc obliczymy sumy sufiksowe rozkładu zmiennej S - tak, aby dla float r zachodziło
        S[ (int) r ] = a^3 * P( S > r ) = V[3][ (int) r+1 ] + ... V[3][M-1].
    */

    S[ M-1 ] = 0;
    for( int j = M-2; j >= 0; --j ) S[j] = S[j+1] + V[3][j+1];

    int P = 0;
    for( unsigned int j = 0; j < b; ++j )
    {
        unsigned int R = (unsigned int) ( 3*t*j );
        if( M <= R ) break;
        else P += S[ R ];
    }

    return P * 1.0 / (a*a*a*b);
}

Wynikiem działania funkcji f( a, b, t ) jest liczba

\(\displaystyle{ \mathcal{P} \left( \frac{X+Y+Z}{3} > t \cdot W \right),}\)

gdzie \(\displaystyle{ X, Y, Z}\) są zmiennymi losowymi o rozkładzie \(\displaystyle{ \mathrm{rand}(a),}\) a \(\displaystyle{ W}\) o rozkładzie \(\displaystyle{ \mathrm{rand}(b).}\)

Wielkość

\(\displaystyle{ \mathbb{E} \left( \frac{X+Y+Z}{3} - t \cdot W \right)}\)

można zaś obliczyć matematycznie, bo wartość oczekiwana jest liniowa:

\(\displaystyle{ {\mathbb{E} \left( \frac{X+Y+Z}{3} - t \cdot W \right) = \frac{ \mathbb{E} X + \mathbb{E} Y + \mathbb{E} Z }{3} - t \cdot \mathbb{E} W = \frac{ \frac{a-1}{2} + \frac{a-1}{2} + \frac{a-1}{2} }{3} - t \cdot \frac{b-1}{2}} \\[1ex]
= \frac{ (a-1) - t \cdot (b-1) }{2}.}\)

[Ravy]

Dopiero teraz zauważyłem, że odpisałeś

Będę zerkał i upewnie się czy to co napisałeś liczy poprawnie.

Masz jakąś sugestię jak jeszcze w miarę wygodnie policzyć warunkową wartość oczekiwaną przy zdarzeniu? Tzn. jaka jest wartość oczekiwana tej zmiennej losowej E((X+Y+Z)/3 - t*W | (X+Y+Z)/3 - t*W > 0)?

Pewnie z tego wzoru \(\displaystyle{ E(X|A)=\frac {1} {P(A)} \sum_{x_i \in A}^{} p_i x_i}\)?

Czyli musiałbym znów policzyć wszystkie prawdopodobieństwa, że (X+Y+Z)/3 - t*W = k dla k > 0?