Brownian motion z własnym generatorem PRNG

Mam taki kod, który wykonuje wykres błądzenia losowego:

import matplotlib.pyplot as plt
import numpy as np

a = [0] * 99
lst = []

def myarray():
    k = 399955
    x = 557877
    counter = 1999

    for i in range(99):
        k = np.uint32(k ^ x)
        x = np.uint32(x * (k >> 1))
        counter = np.uint32(counter + 1)
        x = np.uint32(x ^ counter)
        result = np.uint32(k^x)
        a[i] = (result/2**32)-0.5
    numpyarray = np.array(a)
    return numpyarray

def main():
    n = 100
    T = 1.

    times = np.linspace(0., T, n)
    dt = times[1] - times[0]

    a = myarray()

    #dB = np.sqrt(dt) * np.random.normal(size=(n - 1,))
    dB = np.sqrt(dt) * a
    B0 = np.zeros(shape=(1,))
    B = np.concatenate((B0, np.cumsum(dB, axis=0)), axis=0)
    plt.plot(times, B)
    plt.show()


if __name__ == '__main__':
    main()

W pętli funkcji myarray jest pewien generator PRNG, nie za dobry, ale wystarczający. Gdy zainicjujemy go różnymi k, x, counter daje różne wystarczająco chaotyczne wyniki. Po odkomentowaniu dB można użyć też jakiegoś pythonowego generatora.

Ale chcę takich symulacji błądzenia wykonać wiele na jednym wykresie, dlatego pracuję nad takim kodem:

import matplotlib.pyplot as plt
import numpy as np

a = [0] * 100
lst = []

def myarray():
    for j in range(99):
        k = 399955 + j
        x = 557877 + j
        counter = 1999 + j

        for i in range(100):
            k = np.uint32(k ^ x)
            x = np.uint32(x * (k >> 1))
            counter = np.uint32(counter + 1)
            x = np.uint32(x ^ counter)
            result = np.uint32(k^x)
            a[i] = (result/2**32)-0.5
        lst.append(a)
    numpyarray = np.array(lst)
    return numpyarray

def main():
    n = 100
    d = 100
    T = 1.

    times = np.linspace(0., T, n)
    dt = times[1] - times[0]

    a = myarray()

    #dB = np.sqrt(dt) * np.random.normal(size=(n - 1, d))
    dB = np.sqrt(dt) * a
    B0 = np.zeros(shape=(1, d))
    B = np.concatenate((B0, np.cumsum(dB, axis=0)), axis=0)
    plt.plot(times,B)
    plt.show()

if __name__ == '__main__':
    main()

Po odkomentowaniu dB można zobaczyć, co chcę uzyskać. Ale z jakichś powodów ten program z moim własnym generatorem w myarray dołącza mi do listy lst ciągle tę samą listę a, pomimo, że początkowe k, x i counter się zmieniają, więc pętla powinna tworzyć za każdym razem inną listę a i dołączać ją do listy lst.

Wychodzi zbiór prostych kresek. Pomimo, że pierwszy kod generuje prawidłowo błądzenie losowe. Skąd te proste linie, skoro po wartościach w listach a wydaje się, że liczby są przypadkowe? Co prawda wszystkie listy a wychodzą identyczne (nie wiem jakim cudem), ale wartości w listach a są błądzeniem. Nic tu nie ma sensu. Spodziewałbym się na wykresie raczej 100 identycznych wykresów błądzenia losowego, a nie 100 różnych prostych kresek. Może do wykresu są brane kolumny zamiast wierszy? Więc może powinienem wygenerować kolumny liczb losowych, albo transponować macierz, którą mam? Nie wiem.

Tak na oko losujesz cały czas tym samym seedem, ziarnem tym samym za każdym razem, to te same liczby cały czas powstają.
Ziarno nie powinno być ustawiane od nowe przy każdym uruchomieniu funkcji tylko gdzieś na początku programu.

@Szalony Programista2: ale po to mam właśnie:

for j in range(99):
        k = 399955 + j
        x = 557877 + j
        counter = 1999 + j

Za każdym razem zmieniam seedy - dodaję j (te 3 zmienne można traktować jako 3 seedy). To jest generator wymyślony na poczekaniu, ale jest wystarczająco chaotyczny, żeby dawać zupełnie inne wyniki, nawet, gdy zwiększamy seedy o jeden. Za każdym razem powinienem więc dostawać inną listę a.

Dodam jeszcze, że tutaj:

import matplotlib.pyplot as plt
import numpy as np

a = [0] * 100
lst = []

def myarray():
    for j in range(99):
        k = 399955 + j
        x = 557877 + j
        counter = 1999 + j

        for i in range(100):
            k = np.uint32(k ^ x)
            x = np.uint32(x * (k >> 1))
            counter = np.uint32(counter + 1)
            x = np.uint32(x ^ counter)
            result = np.uint32(k^x)
            a[i] = (result/2**32)-0.5
        lst.append(a)
        #wyszysc liste a

Być może powinienem jakoś czyścić listę a. Tak, żeby przy kolejnej pętli lista a była od nowa uzupełniana. Może coś jest nie tak z nadpisywaniem do tej listy. Ale nie umiem tego zrobić, bo np. napisanie tam:

a = [0] * 100

Zwraca jakieś błędy, których proweniencji też nie rozumiem. Niby czemu nie mogę przypisywać do listy a znów samych zer? Podobnie jest problem z metodą clear(). Ale nie wiem, czy to w ogóle właściwa przyczyna problemów.

@Szalony Programista2: zrobiłem to na tę chwilę tak:

import matplotlib.pyplot as plt
import numpy as np

a = [0] * 100
lst = []

def myarray():
    for j in range(99):
        k = 399955 + j
        x = 557877 + j
        counter = 1999 + j
        a = [0] * 100 #wypelnienie a samymi zerami

        for i in range(100):
            k = np.uint32(k ^ x)
            x = np.uint32(x * (k >> 1))
            counter = np.uint32(counter + 1)
            x = np.uint32(x ^ counter)
            result = np.uint32(k^x)
            a[i] = (result/2**32)-0.5
        lst.append(a)
    numpyarray = np.array(lst)
    return numpyarray

def main():
    n = 100
    d = 100
    T = 1.

    times = np.linspace(0., T, n)
    dt = times[1] - times[0]

    a = myarray()

    #dB = np.sqrt(dt) * np.random.normal(size=(n - 1, d))
    dB = np.sqrt(dt) * a
    B0 = np.zeros(shape=(1, d))
    B = np.concatenate((B0, np.cumsum(dB, axis=0)), axis=0)
    plt.plot(times,B)
    plt.show()

if __name__ == '__main__':
    main()

Coś chyba nadal jest nie tak, bo zawsze jedna linia jest prostą, a do tego wszystkie a[0] są takie same.

def main():
    n = 100
    d = 100
    T = 1.

    times = np.linspace(0., T, n)
    dt = times[1] - times[0]

    a = myarray()

    dB = np.sqrt(dt) * a.T
    B = np.cumsum(dB, axis=0)

    plt.plot(times,B)
    plt.show()

Możesz wyrzucić te dodatkowe concatenate, a transpoze na myarray() wyniku wykonać.

Jak na razie, to mi znowu nie działa:

a = [[0] * 10] * 10
lst = []

def myarray():
    for j in range(10):
        k = 39995 + j
        x = 55787 + j
        counter = 1999 + j

        for i in range(10):
            k = np.uint32(k + x)
            x = np.uint32((x * (k >> 1))
            counter = np.uint32(counter + 1)
            x = np.uint32(x ^ counter)
            result = np.uint32(k^x)
            a[i][j] = result
    numpyarray = np.array(a)
    print(numpyarray)
    return numpyarray

Jak to jest możliwe, że to zwraca listę 10 takich samych list? Dlaczego w nadrzędnej pętli ignorowana jest zmiana seedów:

k = 39995 + j
x = 55787 + j
counter = 1999 + j

I program liczy 10 razy dla tych samych seedów? Wszystko w kodzie wskazuje, że seedy są zmieniane. Znowu jest jakiś problem z zapisem zwykłych liczb do listy. Coś jest nie tak z tymi listami w Pythonie.

Ale ten kod co miałeś na początku generuje poprawny Brownian motion.

import matplotlib.pyplot as plt
import numpy as np

lst = []

def myarray():
    for j in range(99):
        k = 399955 + j
        x = 557877 + j
        counter = 1999 + j
        a = [0] * 100 #wypelnienie a samymi zerami

        for i in range(100):
            k = np.uint32(k ^ x)
            x = np.uint32(x * (k >> 1))
            counter = np.uint32(counter + 1)
            x = np.uint32(x ^ counter)
            result = np.uint32(k^x)
            a[i] = (result/2**32)-0.5
        lst.append(a)
    numpyarray = np.array(lst)
    return numpyarray

def main():
    n = 100
    d = 100
    T = 1.

    times = np.linspace(0., T, n)
    dt = times[1] - times[0]

    a = myarray().T

    #dB = np.sqrt(dt) * np.random.normal(size=(n - 1, d))
    dB = np.sqrt(dt) * a
    B0 = np.zeros(shape=(1, d))
    #B = np.concatenate((B0, np.cumsum(dB, axis=0).T), axis=0)
    B = np.cumsum(dB, axis=0)

    plt.plot(times,B)
    plt.show()

if __name__ == '__main__':
    main()