Spot detection using Markov field model with UBN (Uniformly Boosted Normal) distribution

model.stan

data {
  int I;
  int J;
  matrix[I,J] Y;
  real<lower=0> U;
}

parameters {
  matrix<lower=-5, upper=5>[I,J] mu;
  real<lower=0> s_mu;
  real<lower=0> s_Y;
}

model {
  for (j in 2:J)
    for (i in 1:I)
      target += log1p_exp(-log(U) + normal_lpdf(mu[i,j] | mu[i,j-1], s_mu));
  for (j in 1:J)
    for (i in 2:I)
      target += log1p_exp(-log(U) + normal_lpdf(mu[i,j] | mu[i-1,j], s_mu));
  to_vector(Y) ~ normal(to_vector(mu), s_Y);
  s_mu ~ normal(0, 1);
  s_Y ~ normal(0, 1);
}

run-model-1.R

library(rstan)

set.seed(123)
I <- 30
J <- 30
mu <- matrix(0, nrow=I, ncol=J)
mu[11:20, 11:20] <- 1.0
noise <- matrix(rnorm(n=I*J, 0, sd=1/3), nrow=I, ncol=J)
Y <- mu + noise

stanmodel <- stan_model(file='model.stan')
data <- list(I=I, J=J, Y=Y, U=0.1)
fit <- vb(stanmodel, data=data, init=function() { list(mu=Y, s_Y=1) }, seed=123)

run-model-2.R

library(rstan)

set.seed(123)
I <- 300
J <- 600
img <- png::readPNG('Enjoy.png')
mu <- 1 - img[,,1]
noise <- matrix(rnorm(n=I*J, mean=0, sd=1/3), nrow=I, ncol=J)
Y <- mu + noise

stanmodel <- stan_model(file='model.stan')
data <- list(I=I, J=J, Y=Y, U=0.1)
fit <- vb(stanmodel, data=data, init=function() { list(mu=Y, s_Y=0.5) }, seed=123, eta=1, adapt_engaged=FALSE)