Как найти лучшее значение для прыжков MCMC? - PullRequest
Новая
       105

Как найти лучшее значение для прыжков MCMC?

0 голосов
/ 10 июля 2019

Я использую MCMC для заполнения области на основе данного образца.

Вы можете увидеть визуальные эффекты алгоритма по ссылке ниже.

https://imgur.com/a/wa5Shkl

Оранжевые точки - это мои смоделированные точки, а зеленые - настоящие точки.

Как хорошо видно на рисунке, алгоритм MCMC не может перепрыгнуть через пустое пространство.

Как рассчитать наилучшее значение случайного прыжка, решить ситуации, подобные этой?

Примечание: код не симпатичен: ')

РЕДАКТИРОВАТЬ: добавлен код 

public List<MCMCSimulation> MCMCSimulate(int nSimulations)
    {
        float[,] m_PDF = (float[,])Pdf.Values;
        int n1Orginal = m_PDF.GetLength(0);
        int n2Orginal = m_PDF.GetLength(1);
        int n1 = n1Orginal + 2;
        int n2 = n2Orginal + 2;
        float[,] pdfMat = new float[n1, n2];
        float[] x1f = new float[n1];
        float[] x2f = new float[n2];
        for (int i = 0; i < n1Orginal; i++)
        {
            for (int j = 1; j < n2Orginal; j++)
            {
                pdfMat[i + 1, j + 1] = m_PDF[i, j];
            }
        }
        for (int i = 0; i < n1; i++)
            x1f[i] = i + 1;
        for (int i = 0; i < n2; i++)
            x2f[i] = i + 1;

        float sig1 = Computation.Std(x1f) / 2;
        float sig2 = Computation.Std(x2f) / 2;


        //Initialize random start of the random walk
        Random randGen = new Random(); // seed (int) System.DateTime.Now.Ticks
        float randN1 = (float)randGen.NextDouble();
        float randN2 = (float)randGen.NextDouble();
        float n10 = randN1 * n1;
        float n20 = randN2 * n2;
        float[] rVector1 = new float[nSimulations];
        float[] rVector2 = new float[nSimulations];
        float pdfValue0 = 0;
        int count = 0;
        int i1, i2;
        while (pdfValue0 <= 0)
        {
            i1 = Computation.MinimumResidual(n10, x1f);
            i2 = Computation.MinimumResidual(n20, x2f);
            pdfValue0 = pdfMat[i1, i2];
            if (pdfValue0 <= 0)
            {
                randN1 = (float)randGen.NextDouble();
                randN2 = (float)randGen.NextDouble();
                n10 = n1 * randN1;
                n20 = n2 * randN2;
                count++;
            }
            if (count > 1e3)
                break;
        }
        rVector1 = Computation.AddConstantToArray(rVector1, n10);
        rVector2 = Computation.AddConstantToArray(rVector2, n20);

        //Simulate
        float sig1r;
        float sig2r;
        float r1a;
        float r2a;
        float pdfValue2;
        float a;
        int acc = 0;
        for (int iSim = 0; iSim < nSimulations - 1; iSim++)
        {
            i1 = Computation.MinimumResidual(rVector1[iSim], x1f);
            i2 = Computation.MinimumResidual(rVector2[iSim], x2f);
            pdfValue0 = pdfMat[i1, i2];
            randN1 = (float)(randGen.Next(-10000, 10000) / 10000.0);
            randN2 = (float)(randGen.Next(-10000, 10000) / 10000.0);
            sig1r = sig1 * randN1;
            sig2r = sig2 * randN2;
            r1a = sig1r + rVector1[iSim];
            r2a = sig2r + rVector2[iSim];
            i1 = Computation.MinimumResidual(r1a, x1f);
            i2 = Computation.MinimumResidual(r2a, x2f);
            pdfValue2 = pdfMat[i1, i2];
            a = (float)Math.Min((pdfValue2 / (pdfValue0 + 1e-10)), 1);

            randN1 = (float)randGen.NextDouble();
            if (randN1 < a)
            {
                rVector1[iSim + 1] = r1a;
                rVector2[iSim + 1] = r2a;
                acc++;
            }
            else
            {
                rVector1[iSim + 1] = rVector1[iSim];
                rVector2[iSim + 1] = rVector2[iSim];
            }
        }
...