Saya menerapkan Markov Chain Montecarlo dengan metropolis dan menggonggong alfa untuk integrasi numerik. Saya telah membuat kelas bernama MCMCIntegrator(). Saya telah memuatnya dengan beberapa atribut, salah satunya adalah pdf dari fungsi (a lambda) yang kami coba integrasikan yang disebut g.

import numpy as np
import scipy.stats as st


class MCMCIntegrator:

    def __init__(self):

        self.g = lambda x: st.gamma.pdf(x, 0, 1, scale=1 / 1.23452676)*np.abs(np.cos(1.123454156))
        self.size = 10000
        self.std = 0.6
        self.real_int = 0.06496359

Ada metode lain di kelas ini, size adalah ukuran sampel yang harus dihasilkan kelas, std adalah standar deviasi Kernel Normal, yang akan Anda lihat dalam beberapa detik. real_int adalah nilai integral dari 1 hingga 2 dari fungsi yang kita integrasikan. Saya telah membuatnya dengan skrip R. Sekarang, untuk masalah.

 def _chain(self, method=None):

        """
            Markov chain heat-up with burn-in

        :param method: Metrpolis or barker alpha
        :return: np.array containing the sample
        """
        old = 0
        sample = np.zeros(int(self.size * 1.5))
        i = 0

        if method:
            def alpha(a, b): return min(1, self.g(b) / self.g(a))

        else:
            def alpha(a, b): return self.g(b) / (self.g(a) + self.g(b))

        while i != len(sample):
            new = st.norm(loc=old, scale=self.std).rvs()
            new = abs(new)
            al = alpha(old, new)
            u = st.uniform.rvs()

            if al > u:
                sample[i] = new
                old = new
                i += 1

        return np.array(sample)

Di bawah metode ini adalah metode integrate() yang menghitung proporsi angka dalam interval [1, 2]:

    def integrate(self, method=None):
        """
            Integration step

        """

        sample = self._chain(method=method)
        
        # discarding 30% of the sample for the burn-in
        ind = int(len(sample)*0.3)
        sample = sample[ind:]
        setattr(self, "sample", sample)

        sample = [1 if 1 < v < 2 else 0 for v in sample]
        return np.mean(sample)

Ini adalah fungsi utama:

def main():

    print("-- RESULTS --".center(20), end='\n')
    mcmc = MCMCIntegrator()
    print(f"\t{mcmc.integrate()}", end='\n')
    print(f"\t{np.abs(mcmc.integrate() - mcmc.real_int) / mcmc.real_int}")


if __name__ == "__main__":
    main()

Saya terjebak dalam loop while tak terhingga dan saya tidak tahu mengapa ini terjadi.

0
Occhima 15 Mei 2020, 02:28

1 menjawab

Jawaban Terbaik

Beberapa hal... Anda terpaku pada metode rantai karena komputasi alfa mengembalikan NaN, karena g() mengembalikan NaN. Lihatlah pernyataan cetak yang saya masukkan ke dalam kode Anda dan jalankan ...

Tip:

  1. Jadikan g() fungsi kelas seperti chain.

  2. Uji g() pada beberapa nilai pengujian, ada sesuatu yang jelas salah

  3. Jangan mendefinisikan fungsi secara kondisional seperti alpha. Wayyy membingungkan dan rawan kesalahan/sulit untuk dipecahkan. Cukup berikan alfa apa yang dibutuhkan dan Anda juga dapat menjadikannya fungsi kelas alpha(a, b, method=None)
  4. Lihatlah di mana Anda memperbarui i dalam fungsi `_chain'.... Anda menyadari bahwa Anda mempertaruhkan proses perulangan yang panjang karena Anda hanya memperbarui i secara kondisional!
  5. Anda siap menghadapi bencana dengan penggunaan array numpy. Anda mungkin memiliki banyak angka nol setelah data Anda yang sebenarnya karena Anda terlalu banyak menulis daftar nol besar. Anda TIDAK memerlukan array numpy di sini, cukup gunakan daftar kosong python dan tambahkan nilai baru ke dalamnya, baik nol atau satu...berdasarkan apa pun.

Tambahkan beberapa pernyataan cetak saat Anda memecahkan masalah (atau uji unit fungsi Anda). Coba tambahkan saya ke fungsi Anda di bawah ini ... itulah yang saya gunakan untuk mencari tahu apa yang sedang terjadi

def _chain(self, method=None, verbose=True):

    """
        Markov chain heat-up with burn-in

    :param method: Metrpolis or barker alpha
    :return: np.array containing the sample
    """
    old = 0
    sample = np.zeros(int(self.size * 1.5))
    i = 0

    if method:
        def alpha(a, b): return min(1, self.g(b) / self.g(a))

    else:
        def alpha(a, b): 
            if verbose: print(f'g(a): {self.g(a)}, g(b): {self.g(b)}')
            return self.g(b) / (self.g(a) + self.g(b))

    while i != len(sample):
        new = st.norm(loc=old, scale=self.std).rvs()
        new = abs(new)
        al = alpha(old, new)
        u = st.uniform.rvs()
        if verbose: print(f'old: {old:.3f} new: {new:.3f} alpha: {al:.3f} u: {u:.3f}')
        if al > u:
            sample[i] = new
            old = new
            i += 1              # do you really want to conditionally update i?
        sys.exit(-1)            # to breakout of infinite loop...


    return np.array(sample)
1
AirSquid 15 Mei 2020, 00:51