Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > pl.comp.programming > #29076 > unrolled thread

kontynuacja generatory: mersen vs ranlux

Started by"M.M." <mmarszik@gmail.com>
First post2016-10-09 06:26 -0700
Last post2016-10-11 21:47 +0200
Articles 10 — 3 participants

Back to article view | Back to pl.comp.programming

Contents

  kontynuacja generatory: mersen vs ranlux "M.M." <mmarszik@gmail.com> - 2016-10-09 06:26 -0700
    Re: kontynuacja generatory: mersen vs ranlux bartekltg <bartekltg@gmail.com> - 2016-10-09 15:55 +0200
      Re: kontynuacja generatory: mersen vs ranlux "M.M." <mmarszik@gmail.com> - 2016-10-09 14:12 -0700
        Re: kontynuacja generatory: mersen vs ranlux "M.M." <mmarszik@gmail.com> - 2016-10-10 00:31 -0700
        Re: kontynuacja generatory: mersen vs ranlux bartekltg <bartekltg@gmail.com> - 2016-10-11 22:00 +0200
          Re: kontynuacja generatory: mersen vs ranlux "M.M." <mmarszik@gmail.com> - 2016-10-11 19:14 -0700
    Re: kontynuacja generatory: mersen vs ranlux Borneq <borneq@antyspam.hidden.pl> - 2016-10-11 17:44 +0200
      Re: kontynuacja generatory: mersen vs ranlux "M.M." <mmarszik@gmail.com> - 2016-10-11 11:20 -0700
        Re: kontynuacja generatory: mersen vs ranlux "M.M." <mmarszik@gmail.com> - 2016-10-11 11:53 -0700
      Re: kontynuacja generatory: mersen vs ranlux bartekltg <bartekltg@gmail.com> - 2016-10-11 21:47 +0200

#29076 — kontynuacja generatory: mersen vs ranlux

From"M.M." <mmarszik@gmail.com>
Date2016-10-09 06:26 -0700
Subjectkontynuacja generatory: mersen vs ranlux
Message-ID<f646bd6b-6d53-4bdc-9ca2-3c1579d20e7e@googlegroups.com>
Wywaliłem diehardera. Napisałem na kolanie kilka swoich
testów, pewnie moje są gorsze, ale przynajmniej mam jasność
co do zasady działania i wiem kiedy ryzykuję przekroczenie 
wartości parametrów.

Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT 
oblał stabilny test. Odpaliłem kilka testów dla ranluxa,
na oko bylo widać że ranlux zalicza testy lepiej. Odpaliłem
dłuższy test, jak ktoś ma ochotę, niech odpali u siebie.

Wyniki dla std::mt19937_64
test1 24.592034%
test1  1.487763%
test1  4.233741%
test1 26.478715%
test1  4.174345%
test1 27.751500%
test1 18.710311%
test1 27.819010%
test1 39.071225%
test2 23.595284%
test2 48.625664%
test3 45.054474%
test3  0.498713%
test4 37.237054%
test4  5.056990%
test4 11.755519%
test4 45.547699%
test5 15.588600%
test5 20.128564%
test5  9.747300%
test5 41.895282%
test5 25.154676%
test5 17.381929%
test5 20.064510%
test6  2.690305%
test6 14.040318%
test6 15.422749%
test6 17.341370%
test6  7.701529%
test6  6.545549%
test6 23.770691%
test6 18.814217%
test6  8.335626%
test6  0.035693%

Ostatni wynik niby jest niski, ale paradoks dnia urodzin dla 100 dni
wymaga bardzo dużo liczb losowych, a po drugie ten kod może już
tracić stabilność numeryczną dla 99 punktów swobody. No i jeszcze
mogą być jakieś błędy w kodzie, zobaczymy jak długi test
przejdzie ranlux, gdy przejdzie źle, to bedzie wskazówka do
szukania błędu.






#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <random>

typedef double ftyp;
typedef const ftyp cftyp;

typedef __float128 dtyp;
typedef const dtyp cdtyp;


typedef int ityp;
typedef const ityp cityp;

typedef unsigned int utyp;
typedef const utyp cutyp;

typedef long long ltyp;
typedef const long long cltyp;

typedef unsigned long long ultyp;
typedef const unsigned long long cultyp;

#define CHI_STEPS (512)

//Przykłowy generator
class RLin {
private:
    ultyp x;
public:
    RLin() { sRand(12347512); }
    RLin( cultyp seed ) { sRand(seed); }
    void sRand( cultyp seed ) { x = seed; }
    ultyp getUL() { return x = x * 25214903917ULL + 11ULL; }
    utyp getU()  { return (utyp)getUL(); }
    utyp operator()() { return getU(); }
    ftyp getF()  { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};


// Prawie liniowy generator ;-)
class RLin2 {
private:
    ultyp a,b,c,d;
    utyp n;
    static cultyp m[20];
public:
    RLin2( cultyp seed=12347512 ) { sRand(seed); }
    void sRand( ultyp seed ) {
        a = (seed = seed * 25214903917ULL + 11ULL);
        b = (seed = seed * 25214903917ULL + 11ULL);
        c = (seed = seed * 25214903917ULL + 11ULL);
        d = (seed = seed * 25214903917ULL + 11ULL);
        n = 0xFFFFFFFFu;
    }
    ultyp getUL() {
        if( ++n >= 20 ) n = 0;
        cultyp r = (a * 0xEC5ECAB872DF70B2ull + b * 0x36285E818B75145Full + c * 0xEA8DB7BF70F68277ull + d * 0x32592BEF5EC46BA8ull + 0x1D1BF584DBC7A169ull) % m[n];
        d = c;
        c = b;
        b = a;
        return a = r;
    }
    utyp getU()  { return (utyp)getUL(); }
    utyp operator()() { return getU(); }
    ftyp getF()  { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};

cultyp RLin2::m[] = {
0x3EDB94BD5CE7041Eull,
0xCF64D8021AFBE487ull,
0xD7BE02752CCAD058ull,
0xC258C6DFB3B92DD8ull,
0x6DA82BC9D51F561Eull,
0x0CF26705AFC7640Eull,
0xAC25BA196CCC5A2Aull,
0x7E52DA42E37C8323ull,
0xBD72B519BBF6F08Aull,
0x68B84AF1BBA868A6ull,
0x8BC55D980AC0C053ull,
0x639366BAE2A9815Bull,
0x4E0BECBDECFFC92Dull,
0x6FCE8E25E5C82797ull,
0xDBD42E002C7B6D6Dull,
0xBDBFFA31AA65CDE1ull,
0xFE69A7DA2830FA95ull,
0xEA119C87B81B3C45ull,
0x818F544C89F8F908ull,
0x3C968AEF7ED95B40ull,
};


class RMT {
private:
    std::mt19937_64 re;
public:
    RMT( cultyp seed=12347512 ) { sRand(seed); }
    void sRand( cultyp seed ) { re.seed(seed); }
    ultyp getUL() { return re(); }
    utyp getU()  { return (utyp)getUL(); }
    utyp operator()() { return getU(); }
    ftyp getF()  { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};


class RMTSu {
private:
    std::shuffle_order_engine< std::mt19937_64 , 20000 > re;
public:
    RMTSu( cultyp seed=12347512 ) { sRand(seed); }
    void sRand( cultyp seed ) { re.seed(seed); }
    ultyp getUL() { return re(); }
    utyp getU()  { return (utyp)getUL(); }
    utyp operator()() { return getU(); }
    ftyp getF()  { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};



class RRnx {
private:
    std::ranlux48 re;
public:
    RRnx( cultyp seed=12347512 ) { sRand(seed); }
    void sRand( cultyp seed ) { re.seed(seed); }
    ultyp getUL() { return re(); }
    utyp getU()  { return (utyp)getUL(); }
    utyp operator()() { return getU(); }
    ftyp getF()  { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};


//testowany generator
typedef RRnx RndTest;

static ftyp f2( cftyp x ) { return x*x; }

// rozkład x^2
static dtyp dCh2( cdtyp x , cdtyp k ) {
    return pow( 0.5 , 0.5*k ) * pow( x  , 0.5*k-1 ) * exp( -0.5*x ) / tgamma( 0.5*k );
}

// dystrybuanta x^2 numeryczna
static ftyp ddCh2( cdtyp x , cdtyp k ) {
    dtyp step = x / CHI_STEPS;
    if( step == 0 )
        return 0;
    dtyp sum=0;
    dtyp old, v = 0;
    for( dtyp i=0 ; i<=x ; i+=step ) {
        old = v;
        v = dCh2( i , k );
        if( isinf(v) ) v = 0;
        sum += (v + old) * 0.5 * step; // trapezami
    }
    return (ftyp)sum;
//    return (ftyp)(sum < 1-sum ? sum : 1-sum) ; // przedział lewy lub prawy
}


// N-kubełków wypełniamy równomiernie
static ftyp test1( RndTest &rnd , cultyp SIZE, cutyp BUCKET, cutyp SKIP=0 ) {
    ultyp *const freq = new ultyp[BUCKET];
    memset( freq , 0 , sizeof(ultyp) * BUCKET );

    for( ultyp i=0 ; i<SIZE ; i++ ) {
        freq[ rnd() % BUCKET ] ++;
        for( utyp j=0 ; j<SKIP ; j++ ) rnd();
    }
    ftyp sum = 0;
    for( utyp i=0 ; i<BUCKET ; i++ ) {
//        printf("%llu,",freq[i]);
        sum += f2(freq[i] - (ftyp)SIZE/BUCKET) / ((ftyp)SIZE/BUCKET);
    }
//    printf("\n");
    delete[] freq;
    return ddCh2( sum , BUCKET-1 );
}


// całka wielomianu piątego stopnia
static ftyp test2( RndTest &rnd , cultyp SIZE, cutyp SKIP1=0 , cutyp SKIP2=0 ) {
    ultyp ok = 0;
    for( ultyp i=0 ; i<SIZE ; i++ ) {
        cftyp x = rnd.getF() * 2.0 - 1.0;
        for( utyp j=0 ; j<SKIP1 ; j++ ) rnd();
        cftyp y = ((((4.03 * x + 0.35) * x - 5.83) * x - 0.63) * x + 1.81) * x + 0.55; // horner
        ok += y >= rnd.getF();
        for( utyp j=0 ; j<SKIP2 ; j++ ) rnd();
    }
//    printf("%.14lf\n" , ok / (double)SIZE * 2 - 0.82 );
    ftyp sum = 0;
    sum += f2(ok - SIZE*0.41) / (SIZE*0.41);
    ok = SIZE - ok;
    sum += f2(ok - SIZE*0.59) / (SIZE*0.59);
    return ddCh2( sum , 1 );
}


// pięć 7-wymiarowych sfer
static ftyp test3( RndTest &rnd , cultyp SIZE ) {
    cftyp sf[5][8] = {
        {0.6909275387413800,0.2643721119500700,0.3555511147249490,0.5734918751753870,0.2892100607510660,0.3206429267302160,0.6386669614352290,0.1846687629586090*0.1846687629586090},
        {0.7423938774969430,0.2995040355715900,0.4184010286349800,0.6758922111708670,0.3378314392175530,0.6362593002617360,0.3229378267191350,0.1131459328345950*0.1131459328345950},
        {0.3506312669720500,0.6557646378409120,0.3557951149996370,0.7002436001319440,0.3976298568770290,0.6021499079186470,0.2164524640422310,0.1645823939819820*0.1645823939819820},
        {0.4523639138322320,0.3345844920724630,0.3528572631068530,0.7479594776872550,0.6326152258086950,0.2310057792346920,0.3267629486043010,0.1733424318023030*0.1733424318023030},
        {0.3796008176170290,0.4209442320279780,0.6877630952280020,0.3223769583739340,0.4397911571431910,0.7923348305746910,0.2916006030049170,0.1399564120918510*0.1399564120918510}
    };
    ultyp freq[6] = {0,0,0,0,0,0};
    for( ultyp i=0 ; i<SIZE ; i++ ) {
        utyp k=0;
        while( k<5 && sf[k][7] < f2(rnd.getF()-sf[k][0]) + f2(rnd.getF()-sf[k][1]) + f2(rnd.getF()-sf[k][2]) + f2(rnd.getF()-sf[k][3]) + f2(rnd.getF()-sf[k][4]) + f2(rnd.getF()-sf[k][5]) + f2(rnd.getF()-sf[k][6]))
            k++;
        freq[k] ++ ;
    }
    ftyp sum = 0;
    ftyp sump = 0;
    ftyp p;
    for( utyp i=0 ; i<5 ; i++ ) {
        cftyp r = sqrt(sf[i][7]);
        p = 16 * pow(M_PI,3) / 105 * pow(r,7);
        sum += f2( freq[i] - p * SIZE ) / (p * SIZE);
        sump += p;
    }
    p = (1.0 - sump);
    sum += f2( freq[5] - p * SIZE ) / (p * SIZE);
    return ddCh2( sum , 5 );
}

//przedziały
static ftyp test4( RndTest &rnd , cultyp SIZE , cutyp SKIP=0) {
    cftyp rngs[]  = {0,0.1249265665,0.1729261032,0.2464769642,0.3476090933,0.3792258092,0.4183281037,0.432585527,0.5611377299,0.6213443817,0.6453121416,0.680396048,0.7939877372,0.8844306844,0.9670833906,1.00000000001};
    cutyp BUCKETS = sizeof(rngs) / sizeof(rngs[0]) - 1;
    ultyp freq[BUCKETS];
    memset( freq , 0 , sizeof(freq) );
    for( ultyp i=0 ; i<SIZE ; i++ ) {
        cftyp v = rnd.getF();
        for( utyp j=0 ; j<SKIP ; j++ ) rnd();
        utyp j=1;
        while( v > rngs[j] )
            j++;
        freq[j-1] ++ ;
    }
    ftyp sum = 0;
    for( utyp i=0 ; i<BUCKETS ; i++ )
        sum += f2( freq[i] - (rngs[i+1] - rngs[i]) * SIZE ) / ((rngs[i+1] - rngs[i]) * SIZE);
    return ddCh2( sum , BUCKETS-1 );
}


//cztery bity
static ftyp test5( RndTest &rnd , cultyp SIZE , utyp B1=0, utyp B2=4, utyp B3=16, utyp B4=24, cutyp SKIP=0 ) {
    cutyp BUCKETS = 16;
    ultyp freq[BUCKETS];
    memset( freq , 0 , sizeof(freq) );
    cutyp bits1[4] = {B1,B2,B3,B4};
    cutyp bits2[4] = {1u<<B1,1u<<B2,1u<<B3,1u<<B4};
    for( ultyp i=0 ; i<SIZE ; i++ ) {
        utyp v=0;
        for( utyp j=0 ; j<4 ; j++ ) {
            cutyp r = rnd();
            v = (v << 1) | (( r & bits2[j] ) >> bits1[j]);
            for( utyp k=0 ; k<SKIP ; k++ )
                rnd();
        }
        freq[v]++;
    }
    ftyp sum = 0;
    for( utyp i=0 ; i<BUCKETS ; i++ ) {
//        printf("%llu,",freq[i]);
        sum += f2( freq[i] - (ftyp)SIZE/BUCKETS ) / ((ftyp)SIZE/BUCKETS);
    }
    return ddCh2( sum , BUCKETS-1 );

}


//BIRTHDAY PROBLEM
static ftyp test6( RndTest &rnd , cultyp SIZE , cutyp BUCKETS, cutyp SKIP=0 ) {
    ultyp *const freq = new ultyp[BUCKETS];
    memset( freq , 0 , sizeof(ultyp) * BUCKETS );
    ftyp *const p = new ftyp[BUCKETS]; // rozkład dla paradoksu dnia urodziń
    char *const days = new char[BUCKETS];
    for( utyp i=0 ; i<BUCKETS ; i++ ) {
        p[i] = (ftyp)(i+1) / BUCKETS;
        for( utyp j=0 ; j<i ; j++ ) {
            p[i] *= 1.0 - ( (ftyp)(j+1) / BUCKETS );
        }
//        printf("%.12lf\n",(double)p[i]);
    }

    ftyp sum = 0;
    for( utyp i=0 ; i<=BUCKETS ; i++ )
        sum += p[i];
    if( fabs( sum - 1.0 ) > 0.000000001 )
        abort();


    for( ultyp i=0 ; i<SIZE ; i++ ) {
        memset( days , 0 , sizeof(char)*BUCKETS );
        utyp v=0;
        while(true) {
            cutyp r = rnd() % BUCKETS;
            if( days[r] )
                break;
            days[r] = 1;
            v++;
        }
        if( v-1 >= BUCKETS ) abort();
        freq[v-1] ++;
    }
//    for( utyp i=0 ; i<BUCKETS ; i++ )
//        printf("%llu," , freq[i] );

    sum = 0;
    for( utyp i=0 ; i<BUCKETS ; i++ )
        sum += f2( freq[i] - SIZE * p[i] ) / (SIZE*p[i]);
//    printf("\nsum %lf\n" , (double)sum );

    delete[] p;
    delete[] days;
    delete[] freq;
    return ddCh2( sum , BUCKETS-1 );
}



ftyp rl( cftyp chi ) {
    return (chi < 1-chi ? chi : 1-chi) * 100;
}

int main( int argc , char *argv[] ) {
    RndTest rnd( time(NULL) );

    cultyp m=1000;

    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*1000000 , 16 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*1000000 , 10 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*1000000 , 16 , 1 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*1000000 , 16 , 2 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*1000000 , 16 , 3 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*1000000 , 20 , 4 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*1000000 , 32 , 5 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*300000  , 16 , 20 ) )); fflush(stdout);
    printf("test1 %9.6lf%%\n" , rl(test1( rnd  , m*100000  , 16 , 80 ) )); fflush(stdout);

    printf("test2 %9.6lf%%\n" , rl(test2( rnd  , m*100000 ) ));fflush(stdout);
    printf("test2 %9.6lf%%\n" , rl(test2( rnd  , m*100000 ) ));fflush(stdout);

    printf("test3 %9.6lf%%\n" , rl(test3( rnd  , m*100000 ) ));fflush(stdout);
    printf("test3 %9.6lf%%\n" , rl(test3( rnd  , m*100000 ) ));fflush(stdout);

    printf("test4 %9.6lf%%\n" , rl(test4( rnd  , m*100000 ) ));fflush(stdout);
    printf("test4 %9.6lf%%\n" , rl(test4( rnd  , m*100000 , 1 ) ));fflush(stdout);
    printf("test4 %9.6lf%%\n" , rl(test4( rnd  , m*100000 , 2 ) ));fflush(stdout);
    printf("test4 %9.6lf%%\n" , rl(test4( rnd  , m*100000 , 3 ) ));fflush(stdout);

    printf("test5 %9.6lf%%\n" , rl(test5( rnd  , m*100000 ) ));fflush(stdout);
    printf("test5 %9.6lf%%\n" , rl(test5( rnd  , m*100000 , 0 , 1 , 2 , 3 ) ));fflush(stdout);
    printf("test5 %9.6lf%%\n" , rl(test5( rnd  , m*100000 , 0 , 1 , 2 , 3 , 20 ) ));fflush(stdout);
    printf("test5 %9.6lf%%\n" , rl(test5( rnd  , m*100000 , 0 , 1 , 2 , 3 , 80 ) ));fflush(stdout);
    printf("test5 %9.6lf%%\n" , rl(test5( rnd  , m*100000 , 2 , 4 , 6 , 8 ) ));fflush(stdout);
    printf("test5 %9.6lf%%\n" , rl(test5( rnd  , m*100000 , 2 , 4 , 6 , 8 , 20 ) ));fflush(stdout);
    printf("test5 %9.6lf%%\n" , rl(test5( rnd  , m*100000 , 2 , 4 , 6 , 8 , 80 ) ));fflush(stdout);

    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 ,  6 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 ,  7 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 10 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 11 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 14 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 15 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 16 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 17 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 31 ) )); fflush(stdout);
    printf("test6 %9.6lf%%\n" , rl(test6( rnd  , m*10000000 , 100 ) )); fflush(stdout);

    return 0;
}

[toc] | [next] | [standalone]

#29077

Frombartekltg <bartekltg@gmail.com>
Date2016-10-09 15:55 +0200
Message-ID<ntdi7m$boj$1@node1.news.atman.pl>
In reply to#29076
On 09.10.2016 15:26, M.M. wrote:
> Wywaliłem diehardera. Napisałem na kolanie kilka swoich
> testów, pewnie moje są gorsze,

Jak wspominałęm, cppCon w "telewizji" leci.
Oglądam do kotlera powolutku co IMHO ciekawsze odcinki.
M.in ten:
https://www.youtube.com/watch?v=6DPkyvkMkk8
W sumie niewiele się dowiedziałem, to wstęp do użycia
<random> (choć jak ktoś nie zna za dobrze, pewnie warto
obejrzeć). Wytępuje tam jednak jedna powracająca dygresja.

Nie jesteś spacjalistą, nie rób własnego PRNG. Popsujesz.
Jesteś zajebistym programistą? A jesteś przy okazji
matematykiem/statystykiem? Albo chociaż przeczytałęś
z pełnym zrozumieniem tę górę prac na ten temat Nie?
Prawdopodobnie coś zwalisz.
Historia usiana jest dobrymi na pozór generatorami,
które był← powszechnie używane, a potem okazywało się,
że są OKDP. Z RANDU na czele ;-)

Moja intuicja podpowiada mi, że z testowaniem jest
identycznie ;-)
Zauważ*) że i w dieharder znajdują się tsty określane
jako niepoprawne! A ktoś je tam kiedyś wsadził myśląc,
że są dobre. Spac z dziedziny!


Tak sobie marudzę.
Przez cały post!


 > ale przynajmniej mam jasność
 > co do zasady działania

<złośliwość> Jakbyś dokłądnie p[rzeczytał dokumentację,
wiedziałbyś też jak działa dieharder <\złośliwosć>
;-)

*) mozna to zauważyć czytając pdfy od dieharder;>



> Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
> oblał stabilny test.

A MT idealny nie jest... ;>

> Odpaliłem kilka testów dla ranluxa,
> na oko bylo widać że ranlux zalicza testy lepiej.

Co to znaczy lepiej?
Nie opisałeś kryterium 'dobroci' wyniku.
To, że daje wynik bliższy całce nie musi oznaczać,
że jest lepszy ;>


>  Odpaliłem
> dłuższy test, jak ktoś ma ochotę, niech odpali u siebie.
>
> Wyniki dla std::mt19937_64
> test1 24.592034%
...
> test6  0.035693%

Co oznaczają te procenty?

Co dokładnie robią te testy?

Jasna, mogę poczytać kod. Ale mi się nie chce ;>

Chyba nie wkleiłąś wyników ranluxa.

> Ostatni wynik niby jest niski, ale paradoks dnia urodzin dla 100 dni
> wymaga bardzo dużo liczb losowych,

Gdyby tylko istniał jakiś sposób na zweryfikowanie, czy taki
wynik mocno odstaje od normy... gdyby mo można było
przypisać p-value, przedział ufnośći czy coś takiego...
  ;-)


> a po drugie ten kod może już
> tracić stabilność numeryczną dla 99 punktów swobody. No i jeszcze
> mogą być jakieś błędy w kodzie, zobaczymy jak długi test
> przejdzie ranlux, gdy przejdzie źle, to bedzie wskazówka do
> szukania błędu.

Zaraz, bo się gubię. Jeszcze raz: to jakie są zalety używania
własnego kodu w stosunku do dieharder? ;-)
Myślałem, ze zaletami, o których pisałeś na początku są
właśnie stabilność i dokłądna zjanomość granic, w jakich
testy zachowują się porządnie.

pzdr
bartekltg

[toc] | [prev] | [next] | [standalone]

#29078

From"M.M." <mmarszik@gmail.com>
Date2016-10-09 14:12 -0700
Message-ID<36cac237-fa1a-4622-93ff-1cf4bfb42696@googlegroups.com>
In reply to#29077
On Sunday, October 9, 2016 at 3:55:03 PM UTC+2, bartekltg wrote:
> On 09.10.2016 15:26, M.M. wrote:
> > Wywaliłem diehardera. Napisałem na kolanie kilka swoich
> > testów, pewnie moje są gorsze,
> 
> Jak wspominałęm, cppCon w "telewizji" leci.
> Oglądam do kotlera powolutku co IMHO ciekawsze odcinki.
> M.in ten:
> https://www.youtube.com/watch?v=6DPkyvkMkk8
> W sumie niewiele się dowiedziałem, to wstęp do użycia
> <random> (choć jak ktoś nie zna za dobrze, pewnie warto
> obejrzeć). Wytępuje tam jednak jedna powracająca dygresja.
> 
> Nie jesteś spacjalistą, nie rób własnego PRNG. Popsujesz.
> Jesteś zajebistym programistą? A jesteś przy okazji
> matematykiem/statystykiem? Albo chociaż przeczytałęś
> z pełnym zrozumieniem tę górę prac na ten temat Nie?
> Prawdopodobnie coś zwalisz.

Nie można użyć gotowca, dieharder ma zwalone testy a
dla tych poprawnych wyświetla że generator ich nie
zalicza. Im dłuższy test, tym większe prawdopodobieństwo
że dieharder wyświetli failed. Nie twierdzę że przez
dwie godziny napisałem coś super rewelacyjnego, ale
chyba nie mam wyboru.


> Historia usiana jest dobrymi na pozór generatorami,
> które był← powszechnie używane, a potem okazywało się,
> że są OKDP. Z RANDU na czele ;-)
Nie mówię że nie.


> Moja intuicja podpowiada mi, że z testowaniem jest
> identycznie ;-)
Gdy wady wyjdą, to pomyśli się nad poprawkami. Gdy
wady wyszły w dieharderze to nie ma sensu abym nawet
myślał nad poprawkami.


> Zauważ*) że i w dieharder znajdują się tsty określane
> jako niepoprawne! A ktoś je tam kiedyś wsadził myśląc,
> że są dobre. Spac z dziedziny!
No wiem, tak to już bywa w tworzeniu oprogramowaniu. 
Pewne jest, że ja tego kodu nie poprawię. A jak są
nie działa, to używać też nie mogę.


> Tak sobie marudzę.
> Przez cały post!
Lubię konstruktywną krytykę. 


> 
> 
>  > ale przynajmniej mam jasność
>  > co do zasady działania
> 
> <złośliwość> Jakbyś dokłądnie p[rzeczytał dokumentację,
> wiedziałbyś też jak działa dieharder <\złośliwosć>
> ;-)
> 
> *) mozna to zauważyć czytając pdfy od dieharder;>
Nie sądzę abym zrozumiał w kilka chwil na tyle, aby
go poprawić. Swój kod mogę w kilka godzin napisać.


> > Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
> > oblał stabilny test.
> 
> A MT idealny nie jest... ;>
Idealny nie jest, ale żeby FAILED?



> > Odpaliłem kilka testów dla ranluxa,
> > na oko bylo widać że ranlux zalicza testy lepiej.
> 
> Co to znaczy lepiej?
Tylko to co napisałem - na oko. Częściej było prawdopodobieństwo
z przedziału 10-40% niż w MT.


> Nie opisałeś kryterium 'dobroci' wyniku.
Porównanie z oczekiwanym rozkładem i całka po rozkładzie 
chi-kwadrat.


> To, że daje wynik bliższy całce nie musi oznaczać,
> że jest lepszy ;>
I tak i nie. Generalnie porównanie rozkładu jest
lepsze. Ale przy pewnych założeniach porównanie do
wartości oczekiwanej całki też jest niczego sobie.


> >  Odpaliłem
> > dłuższy test, jak ktoś ma ochotę, niech odpali u siebie.
> >
> > Wyniki dla std::mt19937_64
> > test1 24.592034%
> ...
> > test6  0.035693%
> 
> Co oznaczają te procenty?
Oznaczają prawdopodobieństwo że rozkład był uzyskany generatorem 
równomiernym.


> Co dokładnie robią te testy?
Na razie sześć prościutkich testów. Generalnie wszystkie budują
jakiś rozkład i porównują uzyskany z teoretycznym. Jest po
jednym zdaniu komentarza do każdego testu, nazwy parametrów są
intuicyjne.


> Jasna, mogę poczytać kod. Ale mi się nie chce ;>
> 
> Chyba nie wkleiłąś wyników ranluxa.
Robią się bardzo długo. Bronek zaczął dyskusję, mógłby 
też odpalić u siebie dla ranluxa - wystarczy skompilować i
uruchomić. Potrzebuję w nocy kompa, więc chyba na jutro też
nie dam wyników, może za dwa dni.



> > Ostatni wynik niby jest niski, ale paradoks dnia urodzin dla 100 dni
> > wymaga bardzo dużo liczb losowych,
> 
> Gdyby tylko istniał jakiś sposób na zweryfikowanie, czy taki
> wynik mocno odstaje od normy... gdyby mo można było
> przypisać p-value, przedział ufnośći czy coś takiego...
>   ;-)
Nie rozumiem. Można policzyć jakie jest prawdopodobieństwo uzyskania
danej odchyłki od rozkładu teoretycznego.



> > a po drugie ten kod może już
> > tracić stabilność numeryczną dla 99 punktów swobody. No i jeszcze
> > mogą być jakieś błędy w kodzie, zobaczymy jak długi test
> > przejdzie ranlux, gdy przejdzie źle, to bedzie wskazówka do
> > szukania błędu.
> 
> Zaraz, bo się gubię. Jeszcze raz: to jakie są zalety używania
> własnego kodu w stosunku do dieharder? ;-)
Będę wiedział co poprawić gdy się okaże że mam literówkę w kodzie.
Po drugie, po 10 godzinach testowania nie ma FAILED dla
dobrych generatorów. Po trzecie nie muszę czytać pdfów i 
blokować testów złych z definicji.


> Myślałem, ze zaletami, o których pisałeś na początku są
> właśnie stabilność i dokłądna zjanomość granic, w jakich
> testy zachowują się porządnie.
Bo póki co widzę, że zachowują się porządnie, a 100 dni
dałem specjalnie aby sprawdzić co się dzieje na krawędzi
możliwości.


Pozdrawiam

[toc] | [prev] | [next] | [standalone]

#29079

From"M.M." <mmarszik@gmail.com>
Date2016-10-10 00:31 -0700
Message-ID<0782abd7-c3fa-40b4-bf55-e07c7756923b@googlegroups.com>
In reply to#29078
ranlux też nie przechodzi niektórych testów z diehardera?



time dieharder -a -g 43 -k 2 -m 1              
#=============================================================================#
#            dieharder version 3.31.1 Copyright 2003 Robert G. Brown          #
#=============================================================================#
   rng_name    |rands/second|   Seed   |
         ranlux|  1.70e+07  |2462611344|
#=============================================================================#
        test_name   |ntup| tsamples |psamples|  p-value |Assessment
#=============================================================================#
   diehard_birthdays|   0|       100|     100|0.88022236|  PASSED  
      diehard_operm5|   0|   1000000|     100|0.90030699|  PASSED  
  diehard_rank_32x32|   0|     40000|     100|0.34202277|  PASSED  
    diehard_rank_6x8|   0|    100000|     100|0.90004615|  PASSED  
   diehard_bitstream|   0|   2097152|     100|0.58330877|  PASSED  
        diehard_opso|   0|   2097152|     100|0.25110439|  PASSED  
        diehard_oqso|   0|   2097152|     100|0.00000000|  FAILED  
         diehard_dna|   0|   2097152|     100|0.00000000|  FAILED  
diehard_count_1s_str|   0|    256000|     100|0.49913662|  PASSED  
diehard_count_1s_byt|   0|    256000|     100|0.10809117|  PASSED  
 diehard_parking_lot|   0|     12000|     100|0.77995133|  PASSED  
    diehard_2dsphere|   2|      8000|     100|0.77360600|  PASSED  
    diehard_3dsphere|   3|      4000|     100|0.90895850|  PASSED  
     diehard_squeeze|   0|    100000|     100|0.88240148|  PASSED  
        diehard_sums|   0|       100|     100|0.10298532|  PASSED  
        diehard_runs|   0|    100000|     100|0.95445937|  PASSED  
        diehard_runs|   0|    100000|     100|0.91343512|  PASSED  
       diehard_craps|   0|    200000|     100|0.35637141|  PASSED  
       diehard_craps|   0|    200000|     100|0.97035063|  PASSED  
 marsaglia_tsang_gcd|   0|  10000000|     100|0.85802568|  PASSED  
 marsaglia_tsang_gcd|   0|  10000000|     100|0.61368471|  PASSED  
         sts_monobit|   1|    100000|     100|0.42974295|  PASSED  
            sts_runs|   2|    100000|     100|0.68531642|  PASSED  
          sts_serial|   1|    100000|     100|0.95842043|  PASSED  
          sts_serial|   2|    100000|     100|0.33993473|  PASSED  
          sts_serial|   3|    100000|     100|0.12476066|  PASSED  
          sts_serial|   3|    100000|     100|0.63641557|  PASSED  
          sts_serial|   4|    100000|     100|0.80343420|  PASSED  
          sts_serial|   4|    100000|     100|0.21776213|  PASSED  
          sts_serial|   5|    100000|     100|0.30776040|  PASSED  
          sts_serial|   5|    100000|     100|0.54439159|  PASSED  
          sts_serial|   6|    100000|     100|0.45257163|  PASSED  
          sts_serial|   6|    100000|     100|0.36794052|  PASSED  
          sts_serial|   7|    100000|     100|0.91790726|  PASSED  
          sts_serial|   7|    100000|     100|0.99814453|   WEAK   
          sts_serial|   8|    100000|     100|0.69159010|  PASSED  
          sts_serial|   8|    100000|     100|0.83825976|  PASSED  
          sts_serial|   9|    100000|     100|0.97426109|  PASSED  
          sts_serial|   9|    100000|     100|0.41829717|  PASSED  
          sts_serial|  10|    100000|     100|0.70102580|  PASSED  
          sts_serial|  10|    100000|     100|0.25888646|  PASSED  
          sts_serial|  11|    100000|     100|0.80886103|  PASSED  
          sts_serial|  11|    100000|     100|0.31149446|  PASSED  
          sts_serial|  12|    100000|     100|0.19057194|  PASSED  
          sts_serial|  12|    100000|     100|0.93981744|  PASSED  
          sts_serial|  13|    100000|     100|0.71842057|  PASSED  
          sts_serial|  13|    100000|     100|0.82889716|  PASSED  
          sts_serial|  14|    100000|     100|0.10300738|  PASSED  
          sts_serial|  14|    100000|     100|0.35830770|  PASSED  
          sts_serial|  15|    100000|     100|0.84164959|  PASSED  
          sts_serial|  15|    100000|     100|0.31474944|  PASSED  
          sts_serial|  16|    100000|     100|0.54017219|  PASSED  
          sts_serial|  16|    100000|     100|0.28910527|  PASSED  
         rgb_bitdist|   1|    100000|     100|0.77142701|  PASSED  
         rgb_bitdist|   2|    100000|     100|0.56358554|  PASSED  
         rgb_bitdist|   3|    100000|     100|0.98327710|  PASSED  
         rgb_bitdist|   4|    100000|     100|0.37516059|  PASSED  
         rgb_bitdist|   5|    100000|     100|0.86456045|  PASSED  
         rgb_bitdist|   6|    100000|     100|0.62414793|  PASSED  
         rgb_bitdist|   7|    100000|     100|0.81959304|  PASSED  
         rgb_bitdist|   8|    100000|     100|0.10188415|  PASSED  
         rgb_bitdist|   9|    100000|     100|0.30269712|  PASSED  
         rgb_bitdist|  10|    100000|     100|0.93517990|  PASSED  
         rgb_bitdist|  11|    100000|     100|0.73689472|  PASSED  
         rgb_bitdist|  12|    100000|     100|0.05511034|  PASSED  
rgb_minimum_distance|   2|     10000|    1000|0.42810939|  PASSED  
rgb_minimum_distance|   3|     10000|    1000|0.04064362|  PASSED  
rgb_minimum_distance|   4|     10000|    1000|0.90358592|  PASSED  
rgb_minimum_distance|   5|     10000|    1000|0.18022751|  PASSED  
    rgb_permutations|   2|    100000|     100|0.94235533|  PASSED  
    rgb_permutations|   3|    100000|     100|0.86690842|  PASSED  
    rgb_permutations|   4|    100000|     100|0.20977449|  PASSED  
    rgb_permutations|   5|    100000|     100|0.80389819|  PASSED  
      rgb_lagged_sum|   0|   1000000|     100|0.16755942|  PASSED  
      rgb_lagged_sum|   1|   1000000|     100|0.55384844|  PASSED  
      rgb_lagged_sum|   2|   1000000|     100|0.37733022|  PASSED  
      rgb_lagged_sum|   3|   1000000|     100|0.36760463|  PASSED  
      rgb_lagged_sum|   4|   1000000|     100|0.47523472|  PASSED  
      rgb_lagged_sum|   5|   1000000|     100|0.00842941|  PASSED  
      rgb_lagged_sum|   6|   1000000|     100|0.34494944|  PASSED  
      rgb_lagged_sum|   7|   1000000|     100|0.63678282|  PASSED  
      rgb_lagged_sum|   8|   1000000|     100|0.32773103|  PASSED  
      rgb_lagged_sum|   9|   1000000|     100|0.99367979|  PASSED  
      rgb_lagged_sum|  10|   1000000|     100|0.11164313|  PASSED  
      rgb_lagged_sum|  11|   1000000|     100|0.06842134|  PASSED  
      rgb_lagged_sum|  12|   1000000|     100|0.61437768|  PASSED  
      rgb_lagged_sum|  13|   1000000|     100|0.66287614|  PASSED  
      rgb_lagged_sum|  14|   1000000|     100|0.94178390|  PASSED  
      rgb_lagged_sum|  15|   1000000|     100|0.55486272|  PASSED  
      rgb_lagged_sum|  16|   1000000|     100|0.44590029|  PASSED  
      rgb_lagged_sum|  17|   1000000|     100|0.93088384|  PASSED  
      rgb_lagged_sum|  18|   1000000|     100|0.82484685|  PASSED  
      rgb_lagged_sum|  19|   1000000|     100|0.52468877|  PASSED  
      rgb_lagged_sum|  20|   1000000|     100|0.23220840|  PASSED  
      rgb_lagged_sum|  21|   1000000|     100|0.40290270|  PASSED  
      rgb_lagged_sum|  22|   1000000|     100|0.19995720|  PASSED  
      rgb_lagged_sum|  23|   1000000|     100|0.31981343|  PASSED  
      rgb_lagged_sum|  24|   1000000|     100|0.68635252|  PASSED  
      rgb_lagged_sum|  25|   1000000|     100|0.01243818|  PASSED  
      rgb_lagged_sum|  26|   1000000|     100|0.40787048|  PASSED  
      rgb_lagged_sum|  27|   1000000|     100|0.72885975|  PASSED  
      rgb_lagged_sum|  28|   1000000|     100|0.16945152|  PASSED  
      rgb_lagged_sum|  29|   1000000|     100|0.57625574|  PASSED  
      rgb_lagged_sum|  30|   1000000|     100|0.48784614|  PASSED  
      rgb_lagged_sum|  31|   1000000|     100|0.59446903|  PASSED  
      rgb_lagged_sum|  32|   1000000|     100|0.63657413|  PASSED  
     rgb_kstest_test|   0|     10000|    1000|0.57321002|  PASSED  
     dab_bytedistrib|   0|  51200000|       1|0.21742195|  PASSED  
             dab_dct| 256|     50000|       1|0.74014075|  PASSED  
Preparing to run test 207.  ntuple = 0
        dab_filltree|  32|  15000000|       1|0.23455416|  PASSED  
        dab_filltree|  32|  15000000|       1|0.09528911|  PASSED  
Preparing to run test 208.  ntuple = 0
       dab_filltree2|   0|   5000000|       1|0.56167378|  PASSED  
       dab_filltree2|   1|   5000000|       1|0.74798469|  PASSED  
Preparing to run test 209.  ntuple = 0
        dab_monobit2|  12|  65000000|       1|0.65575155|  PASSED  

real    82m22.710s
user    82m15.539s
sys     0m5.135s

[toc] | [prev] | [next] | [standalone]

#29086

Frombartekltg <bartekltg@gmail.com>
Date2016-10-11 22:00 +0200
Message-ID<ntjgcu$sgu$1@node2.news.atman.pl>
In reply to#29078
On 09.10.2016 23:12, M.M. wrote:
> On Sunday, October 9, 2016 at 3:55:03 PM UTC+2, bartekltg wrote:
>> On 09.10.2016 15:26, M.M. wrote:
>>> Wywaliłem diehardera. Napisałem na kolanie kilka swoich
>>> testów, pewnie moje są gorsze,
>>
>> Jak wspominałęm, cppCon w "telewizji" leci.
>> Oglądam do kotlera powolutku co IMHO ciekawsze odcinki.
>> M.in ten:
>> https://www.youtube.com/watch?v=6DPkyvkMkk8
>> W sumie niewiele się dowiedziałem, to wstęp do użycia
>> <random> (choć jak ktoś nie zna za dobrze, pewnie warto
>> obejrzeć). Wytępuje tam jednak jedna powracająca dygresja.
>>
>> Nie jesteś spacjalistą, nie rób własnego PRNG. Popsujesz.
>> Jesteś zajebistym programistą? A jesteś przy okazji
>> matematykiem/statystykiem? Albo chociaż przeczytałęś
>> z pełnym zrozumieniem tę górę prac na ten temat Nie?
>> Prawdopodobnie coś zwalisz.
>
> Nie można użyć gotowca, dieharder ma zwalone testy a
> dla tych poprawnych wyświetla że generator ich nie
> zalicza. Im dłuższy test, tym większe prawdopodobieństwo
> że dieharder wyświetli failed. Nie twierdzę że przez
> dwie godziny napisałem coś super rewelacyjnego, ale
> chyba nie mam wyboru.


Co to znaczy "nie mam wyboru"?

Uważasz, że zmajstrujesz test który nie będzie się wykrzaczał
dla rzeczywistych liczb losowych?
Bardzo możliwe.
Ja też umiem napsiać taki test. Licze zera i jedynki, mają
być z grubsza po połowie ;-)

Nie o to chodzi, by test się nie wykrzaczał.
Ważniejsze jest, jak _silny_ jest ten test w obszarze
jego stosowalności (zanim sie wykrzaczy i dla losowych).

Twoim zadaniem byłoby więc nie napisanie testu, który
się nie wykrzacza, ale ma małą 'moc sprawdzającą',
ale testu, który się nie wykrzacza i równie dobrze
jak dieharder wykrywa złe generatory.

To nie jest proste i to nie jest nawet robota na dwa tygodnie;>

Zresztą, procedura postępowania jest w dieharder wyjaśniona.

Jest tryb testowania 'do zepsucia'. Jeśli testowany generator
psuje się dla tego samego rzędu #p_samples co najlepsze
generatory, to uznajemy, że test wyczerpał swoje możliwości
i test należy zaliczyć pozytywnie.




> Gdy wady wyjdą, to pomyśli się nad poprawkami. Gdy
> wady wyszły w dieharderze to nie ma sensu abym nawet
> myślał nad poprawkami.

- To wady w sensie 'test nie ejst dowolnie mocny',
a nie 'test jest popsuty'.




>
>
>> Zauważ*) że i w dieharder znajdują się tsty określane
>> jako niepoprawne! A ktoś je tam kiedyś wsadził myśląc,
>> że są dobre. Spac z dziedziny!
> No wiem, tak to już bywa w tworzeniu oprogramowaniu.
> Pewne jest, że ja tego kodu nie poprawię. A jak są
> nie działa, to używać też nie mogę.

Przeciez działa.
Dokłądnei tak jak opisano w dokumenbtacji;>

>> <złośliwość> Jakbyś dokłądnie p[rzeczytał dokumentację,
>> wiedziałbyś też jak działa dieharder <\złośliwosć>
>> ;-)
>>
>> *) mozna to zauważyć czytając pdfy od dieharder;>
> Nie sądzę abym zrozumiał w kilka chwil na tyle, aby
> go poprawić. Swój kod mogę w kilka godzin napisać.

I co z tego?
Prawdopodobnie będzie to słabsz yest niż dieharder
na domyślnych ustawianiach czy od biedy -m 100.



>>> Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
>>> oblał stabilny test.
>>
>> A MT idealny nie jest... ;>
> Idealny nie jest, ale żeby FAILED?

Pamiętaj, zę ejst tam kilka generatorów MT.
Kilka, bop je modyfikowano. Z jakiś powodów;>

>
>>> Odpaliłem kilka testów dla ranluxa,
>>> na oko bylo widać że ranlux zalicza testy lepiej.
>>
>> Co to znaczy lepiej?
> Tylko to co napisałem - na oko. Częściej było prawdopodobieństwo
> z przedziału 10-40% niż w MT.

To równie dobrze możę znaczyć "gorzej";>

>> Nie opisałeś kryterium 'dobroci' wyniku.
> Porównanie z oczekiwanym rozkładem i całka po rozkładzie
> chi-kwadrat.

Całka?



>
>> To, że daje wynik bliższy całce nie musi oznaczać,
>> że jest lepszy ;>
> I tak i nie. Generalnie porównanie rozkładu jest
> lepsze. Ale przy pewnych założeniach porównanie do
> wartości oczekiwanej całki też jest niczego sobie.

Tak, przy testowaniu genertora do obliczenia tej konkretnej całki ;-)


pzdr
bartekltg

[toc] | [prev] | [next] | [standalone]

#29089

From"M.M." <mmarszik@gmail.com>
Date2016-10-11 19:14 -0700
Message-ID<81728cb8-62e9-481e-b42c-2707c5b0cb1c@googlegroups.com>
In reply to#29086
On Tuesday, October 11, 2016 at 10:00:31 PM UTC+2, bartekltg wrote:
> On 09.10.2016 23:12, M.M. wrote:
> > On Sunday, October 9, 2016 at 3:55:03 PM UTC+2, bartekltg wrote:
> >> On 09.10.2016 15:26, M.M. wrote:
> >>> Wywaliłem diehardera. Napisałem na kolanie kilka swoich
> >>> testów, pewnie moje są gorsze,
> >>
> >> Jak wspominałęm, cppCon w "telewizji" leci.
> >> Oglądam do kotlera powolutku co IMHO ciekawsze odcinki.
> >> M.in ten:
> >> https://www.youtube.com/watch?v=6DPkyvkMkk8
> >> W sumie niewiele się dowiedziałem, to wstęp do użycia
> >> <random> (choć jak ktoś nie zna za dobrze, pewnie warto
> >> obejrzeć). Wytępuje tam jednak jedna powracająca dygresja.
> >>
> >> Nie jesteś spacjalistą, nie rób własnego PRNG. Popsujesz.
> >> Jesteś zajebistym programistą? A jesteś przy okazji
> >> matematykiem/statystykiem? Albo chociaż przeczytałęś
> >> z pełnym zrozumieniem tę górę prac na ten temat Nie?
> >> Prawdopodobnie coś zwalisz.
> >
> > Nie można użyć gotowca, dieharder ma zwalone testy a
> > dla tych poprawnych wyświetla że generator ich nie
> > zalicza. Im dłuższy test, tym większe prawdopodobieństwo
> > że dieharder wyświetli failed. Nie twierdzę że przez
> > dwie godziny napisałem coś super rewelacyjnego, ale
> > chyba nie mam wyboru.
> 
> 
> Co to znaczy "nie mam wyboru"?
Dieharder ma problemy. Ciężko jest określić jakie. Kilka
testów nie działa. Jakie nie działają, inaczej wynika z
dokumentacji, inaczej z pomocy podręcznej. Gdy zwiększam
n-tuples, nie wiem co dokładnie robię, ale podejrzanie
dużo testów nie wychodzi na dobrych generatorach. Nie wiem
czy problemem jest stabilność numeryczna, czy błąd testów,
czy faktycznie generator ma problem. Wolę mieć tester o
mniejszym obszarze stosowalności w którym wiadomo co się
dzieje.
 

> Uważasz, że zmajstrujesz test który nie będzie się wykrzaczał
> dla rzeczywistych liczb losowych?
> Bardzo możliwe.
> Ja też umiem napsiać taki test. Licze zera i jedynki, mają
> być z grubsza po połowie ;-)
No ale aż tak prostego to ja nie mam ;-) Na bazie tamtego
kodu można zrobić kilkadziesiąt-kilkaset testów z sześciu
grup. Wada jest taka, że obsługuje małą ilość punktów 
swobody - ale o tym potem.


> Nie o to chodzi, by test się nie wykrzaczał.
> Ważniejsze jest, jak _silny_ jest ten test w obszarze
> jego stosowalności (zanim sie wykrzaczy i dla losowych).
Nie wiem czy rozumiem. Dla losowych składniki sumy:
(x[i] - E[i])^2 / E[i]
Będą najczęściej dążyły do małych wartości, a małych 
wartości, o ile rozumiem, rozkład i jego całki łatwo
i stabilnie się liczy. Może chodzi o to, że dla prawie 
losowych generatorów typu generator liniowy się wywalą.
Moim zdaniem też nie powinien. Nie podoba mi się to że taki 
program uruchomiony z domyślnymi wartościami ma
prawo się wywalić. Powinien wyświetlić 
prawdopodobieństwo równe zero i tyle. Dopiero dla
zaawansowanych, którzy statystykę mają w paluszku, 
powinien mieć np. niebezpieczne opcje dzięki którym
program działa 10 razy szybciej.


> Twoim zadaniem byłoby więc nie napisanie testu, który
> się nie wykrzacza, ale ma małą 'moc sprawdzającą',
> ale testu, który się nie wykrzacza i równie dobrze
> jak dieharder wykrywa złe generatory.
Największą wadą mojego testera jest brak sprytnej implementacji
całki z rozkładu chi. To pociąga za sobą brak możliwości
policzenia testów z dużą ilością punktów swobody. W sieci
nie znalazłem gotowca, samemu nie chce mi się ulepszać,
bo nie mam motywacji. Bronek ostatnio-często pisał  o 
testach generatorów, a że kiedyś odrobinę się tematem
ciekawiłem, nawet mam generator swojego autorstwa, to dałem
się sprowokować do napisania prostego testu - sorki, ale innej
motywacji do implementacji rozkładu nie mam. Jak podacie
dobrego liba, to podmienię i udostępnię kod. Wracając, 
poza tą wadą mój tester powinien też bardzo dobrze 
wykrywać złe generatory, liniowe od razu wykrywa.


> To nie jest proste i to nie jest nawet robota na dwa tygodnie;>
Hmmmm. Jak ktoś dobrze programuje i może uzyskać odpowiedź
od kolegi dobrego z matematyki i statystyki, to powinien zrobić w
20 dni. Wiadomo, im więcej typów testów i im testy bardziej
sparametryzowane, tym więcej pracy. Dużo zależy od tego, czy ktoś
się tym wcześniej interesował. 

Najpierw sprytna implementacja testu-chi. Jak nie znajdzie w sieci, to
powinien zrobić w 5-7 dni. Potem trzeba opracować testy, jak ktoś się
tym nie interesował wcześniej, to i miesiąc będzie opracowywał.
Do testów trzeba znać rozkłady - więc musi być konsultacja z kimś
od statystyki. Potem implementacja testów, jakiegoś interfejsu...
pesymistycznie 2 miesiące, optymistycznie trzy tygodnie roboty.


> Zresztą, procedura postępowania jest w dieharder wyjaśniona.
Niby jest, ale ja jakoś nadal nie jestem pewny co jest
problemem gdy w podaję duże n-tuples. Nie wiem czy generator
jest zepsuty, czy jest błąd w dieharderze, czy test jest
zepsuty, czy liczy gamma (a więc prawie silia) z 0.5*n-tuples i
numerycznie pada.



> Jest tryb testowania 'do zepsucia'. Jeśli testowany generator
> psuje się dla tego samego rzędu #p_samples co najlepsze
> generatory, to uznajemy, że test wyczerpał swoje możliwości
> i test należy zaliczyć pozytywnie.
Nie rozumiem tego. Najlepszy generator powinien zawsze
przejść test. Z ciekawości włączyłem swój test generatora
fionacciego dla około 40 bilionów liczb. Test to paradoks
dnia urodzin. Metakod testu:

kubełki[30] = {zera};
for( pięć bilionow razy ) {
   dni[30] = {zera};
   cnt = 0;
   while(true) {
     r = rand() % 30;
     if( dni[r] ) break;
     v ++ ;
     dni[r] = 1;
   } 
   kubełki[ v-1 ] ++ ;
}
return sprawdz_rozklad( kubelki );

Spodziewam się, że generator przejdzie test i nie
mam obaw co do stabilności numerycznej - no chyba
że wypełni się tylko jeden kubełek - to tylko
bignum nie straci stabilności.
 



> > Gdy wady wyjdą, to pomyśli się nad poprawkami. Gdy
> > wady wyszły w dieharderze to nie ma sensu abym nawet
> > myślał nad poprawkami.
> 
> - To wady w sensie 'test nie ejst dowolnie mocny',
> a nie 'test jest popsuty'.
Nie rozumiem czemu test nie jest dowolnie mocny. Nie rozumiem z
tych samych powodów co powyżej: składniki 
(x[i] - E[i])^2 / E[i]
w przeciętnym generatorze lepszym od liniowego nie rozjeżdżają
się do +-inf.

> >> Zauważ*) że i w dieharder znajdują się tsty określane
> >> jako niepoprawne! A ktoś je tam kiedyś wsadził myśląc,
> >> że są dobre. Spac z dziedziny!
> > No wiem, tak to już bywa w tworzeniu oprogramowaniu.
> > Pewne jest, że ja tego kodu nie poprawię. A jak są
> > nie działa, to używać też nie mogę.
> 
> Przeciez działa.
> Dokłądnei tak jak opisano w dokumenbtacji;>
Gdy daję m=1 i odpalam 1000 razy test dla mersena, to
wszystkie 1000 testów przechodzi. Gdy daję m=1000 i
jeden test, to nie przechodzi. Dlaczego?
 

> 
> >> <złośliwość> Jakbyś dokłądnie p[rzeczytał dokumentację,
> >> wiedziałbyś też jak działa dieharder <\złośliwosć>
> >> ;-)
> >>
> >> *) mozna to zauważyć czytając pdfy od dieharder;>
> > Nie sądzę abym zrozumiał w kilka chwil na tyle, aby
> > go poprawić. Swój kod mogę w kilka godzin napisać.
> 
> I co z tego?
> Prawdopodobnie będzie to słabsz yest niż dieharder
> na domyślnych ustawianiach czy od biedy -m 100.
Czemu -m 100 to jest od biedy? A jak chcę przetestować
na bilionach liczb? Czy będzie słabszy... no bez 
lepszej implementacji chi-kwadrat nie przetestuje się 
tym programem na dużej ilości kubełków. Poza tym 
myslę że jest dobry. 


> >>> Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
> >>> oblał stabilny test.
> >>
> >> A MT idealny nie jest... ;>
> > Idealny nie jest, ale żeby FAILED?
> 
> Pamiętaj, zę ejst tam kilka generatorów MT.
> Kilka, bop je modyfikowano. Z jakiś powodów;>
Hmmmm może faktycznie trafiłem na jakiś MT do bani...
Ale czemu tysiąc testów na m=1 przeszedł?



> >>> Odpaliłem kilka testów dla ranluxa,
> >>> na oko bylo widać że ranlux zalicza testy lepiej.
> >>
> >> Co to znaczy lepiej?
> > Tylko to co napisałem - na oko. Częściej było prawdopodobieństwo
> > z przedziału 10-40% niż w MT.
> 
> To równie dobrze możę znaczyć "gorzej";>
Na oko to na oko.



> >> Nie opisałeś kryterium 'dobroci' wyniku.
> > Porównanie z oczekiwanym rozkładem i całka po rozkładzie
> > chi-kwadrat.
> 
> Całka?
Dystrybuanta.


> >> To, że daje wynik bliższy całce nie musi oznaczać,
> >> że jest lepszy ;>
> > I tak i nie. Generalnie porównanie rozkładu jest
> > lepsze. Ale przy pewnych założeniach porównanie do
> > wartości oczekiwanej całki też jest niczego sobie.
> 
> Tak, przy testowaniu genertora do obliczenia tej konkretnej całki ;-)
No nie tylko :) Jeśli policzy się wiele różnorodnych całek już
test jest silniejszy. Jeśli całki policzy się przy pomocy każdej
liczby, potem co drugiej, co trzeciej, itd, to znowu wzmacniamy test.
Jeśli całki policzy się xorem z N kolejnych liczb to znowu wzmacniamy
test. Z każdą taką sztuczką wzmacniamy test, a pozbywamy się
numerycznej niestabilności. W końcu każdy generator programy dla
jakiejś maski źle policzy którąś z całek, ciekawe jakby mersen 
policzył calkę gdyby brać co około 623 liczbę.


Pozdrawiam

[toc] | [prev] | [next] | [standalone]

#29082

FromBorneq <borneq@antyspam.hidden.pl>
Date2016-10-11 17:44 +0200
Message-ID<ntj1d1$uah$1@node1.news.atman.pl>
In reply to#29076
W dniu 09.10.2016 o 15:26, M.M. pisze:
> Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
> oblał stabilny test. Odpaliłem kilka testów dla ranluxa,

Test urodzinowy to przybliżenie:
po pierwsze jest rozkład Poissona a ang. wiki mówi że ma być wykładniczy
po drugie - wartość oczekiwana to m^3/4n a gdy m=n to wychodzi n^2/4 
podczas gdy maksymalnie może być n. Z tego wynika że ten wzór to 
przybliżenie.

[toc] | [prev] | [next] | [standalone]

#29083

From"M.M." <mmarszik@gmail.com>
Date2016-10-11 11:20 -0700
Message-ID<8049b394-58c4-472b-a485-23e07eb556af@googlegroups.com>
In reply to#29082
On Tuesday, October 11, 2016 at 5:44:35 PM UTC+2, Borneq wrote:
> W dniu 09.10.2016 o 15:26, M.M. pisze:
> > Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
> > oblał stabilny test. Odpaliłem kilka testów dla ranluxa,
> 
> Test urodzinowy to przybliżenie:
> po pierwsze jest rozkład Poissona a ang. wiki mówi że ma być wykładniczy
> po drugie - wartość oczekiwana to m^3/4n a gdy m=n to wychodzi n^2/4 
> podczas gdy maksymalnie może być n. Z tego wynika że ten wzór to 
> przybliżenie.

Zgłupiałem całkowicie, bo inne testy są opisanie jako podejrzanie na
stronie diehardera i inne wyskoczyły gdy wyświetliłem z linii komend -
to zamieszanie to kolejny argument za tym, aby napisać coś swojego do
testowania generatorów.

Odpowiadając na to co napisałeś o teście urodzinowym, to można napisać
wiele różnych testów, inspirując się faktem, że dzieci rodzą się w 
jakimś dniu tygodnia, miesiąca, czy roku. Bardzo słynny jest paradoks
dnia urodzin (który w istocie nie jest żadnym paradoksem), więc
ja mam test ściśle dotyczący paradoksu dnia urodzin. W rozkładzie
dla paradoksu dnia urodzin prawdopodobieństwo bardzo szybko maleje,
dla 30 dni (w roku) przegranie zakładu gdy w sali jest 29 osób, zdarza
się gdzieś raz na 100 miliardów razy. Dlatego testy powinny być
wykonywane dla bilionów losowań. 

U mnie dla miliarda rozkłady wyglądają tak:

http://s5.ifotos.pl/img/screen1pn_arenenn.png

Jak widać MT przechodzi test z shufflowaniem i bez, fibbonaci też
przechodzi, liniowe absolutnie nie, a na wyniki z ranluxa się
nie doczekałem - wolno działa.

Pozdrawiam

[toc] | [prev] | [next] | [standalone]

#29084

From"M.M." <mmarszik@gmail.com>
Date2016-10-11 11:53 -0700
Message-ID<3550a4a0-65d7-45c3-9ac1-8c728e483738@googlegroups.com>
In reply to#29083
On Tuesday, October 11, 2016 at 8:20:05 PM UTC+2, M.M. wrote:
> Jak widać MT przechodzi test z shufflowaniem i bez, fibbonaci też
> przechodzi, liniowe absolutnie nie, a na wyniki z ranluxa się
> nie doczekałem - wolno działa.

Doczekałem się i dodałem do tabelki ranluxa. Ponadto wrzuciłem
lekko zoptymalizowaną wersję fibonacciego. Widać że generator
fibonnacciego działa prawie tak wydajnie jak generator liniowy, a
wyniki w tym teście są porównywalne do MT i RANLUX.

RANLUX działa koszmarnie wolno.

http://s10.ifotos.pl/img/screen2pn_arensrh.png

Pozdrawiam

[toc] | [prev] | [next] | [standalone]

#29085

Frombartekltg <bartekltg@gmail.com>
Date2016-10-11 21:47 +0200
Message-ID<ntjfkl$rro$1@node2.news.atman.pl>
In reply to#29082
On 11.10.2016 17:44, Borneq wrote:
> W dniu 09.10.2016 o 15:26, M.M. pisze:
>> Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
>> oblał stabilny test. Odpaliłem kilka testów dla ranluxa,
>
> Test urodzinowy to przybliżenie:
> po pierwsze jest rozkład Poissona a ang. wiki mówi że ma być wykładniczy
> po drugie - wartość oczekiwana to m^3/4n a gdy m=n to wychodzi n^2/4
> podczas gdy maksymalnie może być n. Z tego wynika że ten wzór to
> przybliżenie.

Z tego wynika, zę mieszasz dwa testy o podobnej tej samej nazwie!

Jeden mówi*) o rozkładzie odstępów (i to w przybliżeniu ciagłym,
a w rzeczywistośći losujemy z dyskretnego zbioru!)
Drugi o liczbie _powtórzeń_ na liście odstępów.
m^3 /(4n) dotyczy tego drugiego (też jest tylko asymptotyczne).


*) do tego nieźle oszukuje;> różnice wspołrzednych będą
liczbami z wykładniczego o gęstości λ jesli wygenerujemy
na tym przedziale proces Poissona. Np losując liczbę
N z rozkładu Poissona (λ * dlugość przedzialu), po czym
losując jednostajnie N liczb z tego przedziału.
Wtedy różnice pomiędzy kolejnymi (posortowanymi) liczbami
są niezależnymi zmiennymi z rozkłądu wykładniczego.

Jeśli mamy _zadane_ N i losujemy N liczb jednostajnie na odcinku,
ich różnice:
nie są niezależne!
różnice x1-0, x2-x1, x3-x2.. 1-x_n
mają rozkład Dirichleta Dir(1,1,1,1,1,1).

Zapewne dąży do tego, co trzeba, ale też jest
to rozkłąd przybliżony, i to dość nieźle**, ale nadal
tylko asymptotyczny, tak jak w przypadku testu na
powtórzenia.


**)
n=100000000;
r=rand(n,1);  %jednorodny na [0,1]
k=1000;
r=sort(r);
[h,xx]=hist(diff(r),k);
r=[];
semilogy(xx,h/max(xx)/n*k,'.',xx, n*exp(-(n)*xx),'-')

https://www.dropbox.com/s/jlltp0p2chh18bo/b_test.eps?dl=0

Układa się doćś ładnie na krzywej odpowiedniego exp,
ale oczywisćie nie idealnie.




Losowe rzeczy, które się otworzyły z googla

https://www.jstatsoft.org/article/view/v007i03/tuftests.pdf

https://www.cs.indiana.edu/~kapadia/project2/node21.html

https://en.wikipedia.org/wiki/Dirichlet_distribution

pzdr
bartekltg

[toc] | [prev] | [standalone]

Back to top | Article view | pl.comp.programming

csiph-web