Выборка Гиббса дает малые вероятности

Question

Как часть нашего окончательного дизайн-проекта, мы должны спроектировать сэмплер Gibbs для шумоподавления изображения. Мы решили использовать алгоритм Метрополиса вместо обычного сэмплера Гиббса. Примерный набросок алгоритма выглядит следующим образом, все пиксели имеют 0-255 значений оттенков серого. Также мы используем простое сглаживание предварительного распределения.

main()
    get input image as img
    originalImg = img
    for k = 1 to 1000
        beta = 3/log(1+k)
        initialEnergy = energy(img,originalImg)
        for i = 0 to imageRows
            for j = 0 to imageCols
                img[i][j] = metropolisSample(beta,originalImg,img,initialEnergy,i,j)


energy(img,originalImg)
    for i = 1 to imageRows
        for j = 1 to imageCols
            ans += (img[i][j] - originalImg[i][j])^2 / (255*255)
            ans += (img[i][j] - image[i][j+1])^2 / (255*255)
            ans += (img[i][j] - image[i][j-1])^2 / (255*255)
            ans += (img[i][j] - image[i+1][j])^2 / (255*255)
            ans += (img[i][j] - image[i-1][j])^2 / (255*255)
    return ans


metropolisSample(beta,originalImg,img,initialEnergy,i,j)
    imageCopy = img
    imageCopy[i][j] = random int between 0 and 255
    newEnergy = energy(imageCopy,originalImg)
    if (newEnergy < initialEnergy)
        initialEnergy = newEnergy
        return imageCopy[i][j]
    else
        rand = random float between 0 and 1 
        prob = exp(-(1/beta) * (newEnergy - initialEnergy))
        if rand < prob
            initialEnergy = newEnergy
            return imageCopy[i][j]
        else
            return img[i][j]

Это в значительной степени суть программы. Моя проблема в том, что на этапе, где я вычисляю вероятность

prob = exp(-(1/beta) * (newEnergy - initialEnergy))

Разница в энергии настолько велика, что вероятность почти всегда равна нулю. Как правильно смягчить это? Мы также попробовали подход выборки Гиббса, но столкнулись с аналогичной проблемой. Код сэмплера Гиббса выглядит следующим образом. Вместо использования metropolisSample мы используем вместо этого gibbsSample

gibbsSample(beta,originalImg,img,initialEnergy,i,j)
    imageCopy = img
    sum = 0
    for k = 0 to 255
        imageCopy[i][j] = k
        energies[k] = energy(imageCopy,originalImg)
        prob[k] = exp(-(1/beta) * energies[k])
        sum += prob[k]

    for k = 0 to 255
        prob[k] / sum

    for k = 1 to 255
        prob[k] = prob[k-1] + prob[k] //Convert our PDF to a CDF

    rand = random float between 0 and 1
    k = 0
    while (1)
        if (rand < prob[k])
            break
        k++
    initialEnergy = energy[k]
    return k

У нас были похожие проблемы и с этой реализацией. Когда мы рассчитали

prob[k] = exp(-(1/beta) * energies[k])

наши энергии были настолько велики, что все наши вероятности сводились к нулю. Теоретически, это не должно быть проблемой, потому что мы суммируем их все, а затем делим на сумму, но представление с плавающей запятой просто недостаточно точно. Что было бы хорошим способом исправить это?

Answer 1

Я ничего не знаю о вашей конкретной проблеме, но мой первый ответ - масштабирование энергий. Ваши пиксели находятся в диапазоне 0..255, что является произвольным. Если бы пиксели были долями от нуля до единицы, у вас были бы очень разные результаты.

Если единицы энергии в пикселях ^2, попробуйте разделить энергии на 256^2. Иначе попробуйте разделить на 256.

Кроме того, учитывая, что данные являются полностью случайными, возможно, что есть очень высокие энергии, и на самом деле не должно быть высокой вероятности.

Мое незнание вашей проблемы могло привести к бесполезному ответу. Если так, пожалуйста, игнорируйте это.

Answer 2

Я думаю, что вероятность для выборки Гиббса в модели Изинга должна быть

p = 1 / (1 + np.exp(-2 * beta * Energy(x,y)))