R/internal.R

In nicolemunoz99/multisensory: Predicting the conditions under which multisensory integration is favored

#' Difference between stimulus distribution overlap and uncertainty
#'
#' Internal function for calculating means of distributions centered around 0 based on user-specified uncertainties, u1 or u2,
#' by solving for the root in voi. Distribution in state "a" has mean, mu. Distribution in state "b"
#' has mean -mu.
#'
#' @param mu mean of normal distribution in state "a"
#' @param uu uncertainty - proportion of overlap of distributions in states "a" and "b"
#'
#' @return None
#'
#' @examples
#' uniroot(f.mu, c(-10,10), 0.2)
#'
#' @export

f.mu = function(mu, uu){
  # internal function; difference between the proportion of overlap and u1 (or u2).
  # assummes means of distributions are symmetrical about s1 = 0 (i.e., the distributions overlap at 0).
  # this function will be set to zero and solved using uniroot in order to calculate mu.a and mu.b from u1 (or u2).
  area.b = pnorm(0, mu, 1) # area under b distribution to the left of s1 = 0
  area.a = 1 - pnorm(0, -mu, 1) # area under a distribution to the right of s1 = 0
  overlap = area.a + area.b
  overlap-uu # want to know the mu where overlap = uu
}




#' Value of information of Multisensory Stimuli
#'
#' Calculates the value of information for two sequential stimuli based on signal-detection theory and Bayesian inference.
#' See Munoz & Blumstein's "Optimal Integration" (under review at Behavioral Ecology), for derivation and list of assumptions.
#' In Munoz & BLumstein, "b" = predator (PRED) state, "a" = non-threat (NONE) state.
#'
#'
#' @param t.prior.b prior probability of the world in state "b".
#' @param t.bb.b benefit of correct behavior when world is in state "b".
#' @param t.bb.a benefit of correct behavior when world is in state "a".
#' @param t.s1 magnitude of the first stimulus.
#' @param t.u1 uncertainty of the 1st stimulus (proportion of overlap of distributions of 1st stimulus).
#' @param t.u2 uncertainty of the 2nd stimulus (proportion of overlap of distributions of 2nd stimulus).
#' @param t.k1 cost of attending to the 1st stimulus.
#' @param t.k2 cost of attending to the 2nd stimulus.
#'
#' @return List having components int1, vi1, int2, vi2, st.use
#' @return     int1: T/F (was the 1st stimulus used?).
#' @return     vi1:  Value of Information of 1st stimulus.
#' @return     int2: T/F (was the 2nd stimulus used?).
#' @return     vi2:  Value of Information of the 2nd stimulus.
#' @return     st.use: 0 - neither stimulus used; 1 - only 1st stimulus used; 2 - only 2nd stimulus used; 3 - both stimuli used.
#'
#' @examples
#' voi(0.3, 2, 1, -1, 0.2, 0.8, 0, 1)
#'
#' @export


voi <- function(t.prior.b, t.bb.b, t.bb.a, t.s1, t.u1, t.u2, t.k1, t.k2){
  # calculates the Value of Information (VI) for two sequential stimuli and outputs a list of VI for each stimulus,
  # and whether each stimulus was used
  p.b = t.prior.b

  if(!isTRUE(all.equal(0, t.u1))){
    mu.result = uniroot(f.mu, c(-10,10), uu = t.u1)$root
    mu.a1 = -mu.result
    mu.b1 = mu.result
  } else {mu.a1 = -999; mu.b1 = 999}


  if(!isTRUE(all.equal(0, t.u2))){
    mu.result2 = uniroot(f.mu, c(-10,10), uu = t.u2)$root
    mu.a2 = -mu.result2
    mu.b2 = mu.result2
  }else {mu.a2 = -999; mu.b2 = 999}


  if(t.bb.b == 0 && t.bb.a == 0){
    beta = (1-p.b)/p.b
    s1c = (2*(log(p.b/(1-p.b))) - mu.b1^2 + mu.a1^2) / (2*(mu.a1 - mu.b1))
  }	else {
    beta = (1-p.b)*t.bb.a/(p.b*t.bb.b)
    s1c = (2*(log(p.b*t.bb.b/((1-p.b)*t.bb.a))) - mu.b1^2 + mu.a1^2) / (2*(mu.a1 - mu.b1))
  }

  p.hit = 1 - pnorm(s1c, mu.b1, sd = 1)
  p.fa  = 1 - pnorm(s1c, mu.a1, sd =1)


  f.p.1c = function(x){
    # function for solving for cut-off probability p.1c in the absence of stimuli
    x*(1-p.b)/(p.b*(1-x)) - beta
  }

  # solve for p.1c, the cutoff prob in the absence of stimuli and account for the limits of beta
  p.1c = ifelse(beta == Inf, 1,
                ifelse(beta == 0, 0, uniroot(f.p.1c, c(0, 1))$root))

  vi.1.actual = ifelse(p.b >= p.1c,
                       p.b*p.hit*t.bb.b - (1-p.b)*p.fa*t.bb.a + p.b*t.bb.b*(1-beta) - t.k1,
                       p.b*p.hit*t.bb.b - (1-p.b)*p.fa*t.bb.a - t.k1
  )
  vi.1 = vi.1.actual
  int1 = vi.1 >= 0 # is the first stimulus used?


  ########
  # update prior and new cutoff stimulus if first stimulus is used
  if(int1 == T ){
    p.s1 = p.b*dnorm(t.s1, mu.b1, sd = 1) + (1-p.b)*dnorm(t.s1, mu.a1, sd = 1)
    p.b2 = p.b*dnorm(t.s1, mu.b1, sd = 1) / p.s1

  }

  # don't update prior if first stimulus is ignored
  if(int1 == F){
    p.b2 = p.b # prior isn't updated when 1st stimulus is ignored
  }


  ########
  # 2nd stimulus

  if(t.bb.b == 0 && t.bb.a == 0){
    beta2 = (1-p.b2)/p.b2
    s2c = (2*( log( p.b2/(1-p.b2) )) - mu.b2^2 + mu.a2^2) / (2*(mu.a2 - mu.b2))
  }	else {
    beta2 = (1-p.b2)*t.bb.a/(p.b2*t.bb.b)
    s2c = (2*( log( p.b2*t.bb.b/((1-p.b2)*t.bb.a) )) - mu.b2^2 + mu.a2^2) / (2*(mu.a2 - mu.b2))
  }

  p.hit2 = 1 - pnorm(s2c, mu.b2, sd = 1)
  p.fa2 = 1 - pnorm(s2c, mu.a2, sd = 1)


  f.p.2c = function(x){
    x*(1-p.b2)/(p.b2*(1-x)) - beta2
  }

  p.2c = ifelse(beta2 == Inf, 1,
                ifelse(isTRUE(all.equal(0, beta2)), 0, uniroot(f.p.2c, c(0, 1))$root))

  vi.2.actual = ifelse(p.b2 >= p.2c,
                       p.b2*p.hit2*t.bb.b - (1-p.b2)*p.fa2*t.bb.a + p.b2*t.bb.b*(1-beta2) - t.k2,
                       p.b2*p.hit2*t.bb.b - (1-p.b2)*p.fa2*t.bb.a - t.k2
  )

  vi.2 = ifelse(isTRUE(all.equal(0, vi.2.actual)), 0, vi.2.actual)

  int2 = vi.2 >= 0


  ####### which stimulus/stimuli was/were used? (0 = neither, 1 = 1st only, 2 = 2nd only, 3 = both)
  st.use = ifelse (int1 == T && int2 == T, 3, # integration
                   ifelse (int1 == T && int2 == F, 1, # only 1st stimulus used
                           ifelse (int1 == F && int2 == T, 2, # only 2nd stimulus used
                                   0))) # neither stimulus was used

  list(int1 = int1, vi.1 = vi.1, int2 = int2, vi.2 = vi.2, st.use = st.use)

}