A toy model of a neuron with a gaussian receptive field responding to sparse noise, and recovering parameters with Stan

spatial-interactions.R

library(plyr)
library(ggplot2)
library(R.cache)
library(rstan)

#There is a sensor which either emits "yes" or "no." It adds up
#signals over some region of space (The "receptive field", which we
#model as having a Gaussian profile,) and emits a "yes" if the
#summed-up spatially weighted signals, plus unit additive noise, are
#over some threshold.

#We try to characterize this sensor by presenting random arrangements
#of signals in space and seeing whether it signals "yes" on each trial.
#We would like to recover the receptive field location, receptive field size,
#threshold, and sensitivity of the sensor.

#this function simulates a single trial. with a randomly renerated stimulus.
#There may be more than one stimulus source on any given trial.
#(visual neuroscience would refer to this kind of stimulus as "sparse noise")
spatial_interaction_sim <- function(
    trial.id,
    sensitivity = 1,
    threshold = 0,
    receptive.field = 1,
    sensor.position = 0,
    n.sources = rpois(1, lambda=4),
    source.position = runif(n.sources, -20, 20),
    source.strength = rnorm(n.sources, mean=0, sd=1)
    ) {
  if (n.sources == 0) {
    #Some trials have no sources; emit a row in any case
    n.sources <- 1; source.strength <- 0
  }
  signal <- sum(dnorm(source.position, sensor.position, receptive.field)
                * source.strength)
  link <- signal * sensitivity - threshold;
  response <- link > rnorm(1, 0, 1)
  data.frame(source.position, source.strength,
             trial.id = rep(trial.id, n.sources),
             response = rep(response, n.sources),
             link = rep(link), n.sources)
}

#Simulate data for fitting. We provide parameters for sensor position,
#threshold, RF size, and noise, and see if Stan can recover them.
data <- ldply(1:3000, spatial_interaction_sim,
              sensor.position=5, threshold=0.5, receptive.field=2, sensitivity=10)

print(summary(daply(data, "trial.id", with, unique(response))))

#plot the data using white points for stimuli where the response was
#"yes" and black points "no". The region of influence of the sensor is clear.
print(ggplot(data, aes(source.position, source.strength, color=response))
      + geom_point(size=1.5) + scale_color_manual(values=c("black", "white"))
      + theme(panel.background=element_rect(fill="gray50"),
              panel.grid.minor=element_blank(),
              panel.grid.major=element_blank()))

#Another way of looking at the data is color according to stimulus strength
#and facet according to response.
print(ggplot(data, aes(source.position, trial.id, color=source.strength))
      + geom_point(size=1.5)
      + scale_color_gradient(low="#0000FF", high="#FFFF00", space="rgb")
      + facet_grid(response ~.)
      + theme(panel.background=element_rect(fill="gray50"),
              panel.grid.minor=element_blank(),
              panel.grid.major=element_blank()))

#So that should be way more than enough data to recover the sensor parameters.

#function to format data for Stan.
spatial_interaction_stan_format <- function(df) {
  df <- arrange(df, trial.id)
  within(list(), {
    N_trials <- length(unique(df$trial.id))
    N_sources <- length(df$trial.id)
    source_position <- df$source.position
    source_strength <- df$source.strength
    response <- as.numeric(daply(df, "trial.id", with, unique(response)))
    trial_index <- match(df$trial.id, 1:length(unique(df$trial.id)))
  })
}

stan_data <- spatial_interaction_stan_format(data)
str(stan_data)

#The Stan model
model.code <- '
data {
    int N_sources;
    int N_trials;
    real source_position[N_sources];
    real source_strength[N_sources];
    int response[N_trials];
    int trial_index[N_sources];
}
parameters {
    real sensitivity;
    real threshold;
    real<lower=min(source_position), upper=max(source_position)> position;
    real<lower=0, upper=10> field;
}
model {
    vector[N_trials] signal;
    vector[N_trials] link;
    vector[N_trials] prob;
    for (t in 1:N_trials) signal[t] <- 0;
    for (s in 1:N_sources) {
        signal[trial_index[s]] <- signal[trial_index[s]]
            + source_strength[s] * exp(normal_log(source_position[s], position, field));
    }
    link <- signal * sensitivity - threshold;
    for (t in 1:N_trials) {
        prob[t] <- normal_cdf(link[t], 0, 1);
    }
    response ~ bernoulli(prob);
}
'

#R.cache to save compiling time
model <- memoizedCall(stan_model, model_name="receptive_field",
                      model_code=model.code)

samples <- sampling(model, data=stan_data,
                    warmup=500, iter=1000)

opt <- optimizing(model, data=stan_data)