nll_rescaling_agency = function(noise.sigmaLow, extraNoise, crit, bufferFactor, nRNotMe, subjData, data_LowNoise, data_HighNoise) {
  library(here)
  library(tidyr)
  
  noise.sigmaHigh = noise.sigmaLow + extraNoise
  
  nRMe = 6-nRNotMe
  # Get the response counts
  # groupedData = group_by(subjData, Alterations, Noise, SoARatings) %>% 
  #   summarize(count = length(SoARatings))
  # data_LowNoise = dplyr::filter(groupedData, Noise==0)
  # data_HighNoise = dplyr::filter(groupedData, Noise==4)
  
  
  alts = c(0, 0.07, 0.1, 0.2)
  
  # Low Noise
  confBounds_NotMe = crit
  confBounds_Me = crit
  buffer = bufferFactor*noise.sigmaLow
  a = alts[4] + buffer 
  b = alts[1] - buffer
  adist = abs(a-crit)
  bdist = abs(b-crit)
  dist = max(c(adist, bdist))
  notMedist = dist/nRNotMe
  medist = dist/nRMe
  
  nextBound = crit
  for (i in 2:nRNotMe){
    nextBound = nextBound + notMedist
    confBounds_NotMe = c(confBounds_NotMe, nextBound)
  }
  confBounds_NotMe = c(confBounds_NotMe, Inf)
  nextBound = crit
  for (i in 2:nRMe){
    nextBound = nextBound - medist
    confBounds_Me = c(confBounds_Me, nextBound)
  }
  confBounds_Me = c(confBounds_Me, -Inf)
  
  # Calculate probabilities of ratings
  prSoA_alt1 = c(0,0,0,0,0,0)
  prSoA_alt2 = c(0,0,0,0,0,0)
  prSoA_alt3 = c(0,0,0,0,0,0)
  prSoA_alt4 = c(0,0,0,0,0,0)
  
  for (i in 1:nRNotMe){
    prSoA_alt1[i] = (pnorm(confBounds_NotMe[(nRNotMe+2)-i], alts[1], noise.sigmaLow) - pnorm(confBounds_NotMe[(nRNotMe+1)-i], alts[1], noise.sigmaLow))#/(1-pnorm(crit, alts[1], noise.sigmaLow))
    prSoA_alt2[i] = (pnorm(confBounds_NotMe[(nRNotMe+2)-i], alts[2], noise.sigmaLow) - pnorm(confBounds_NotMe[(nRNotMe+1)-i], alts[2], noise.sigmaLow))#/(1-pnorm(crit, alts[2], noise.sigmaLow))
    prSoA_alt3[i] = (pnorm(confBounds_NotMe[(nRNotMe+2)-i], alts[3], noise.sigmaLow) - pnorm(confBounds_NotMe[(nRNotMe+1)-i], alts[3], noise.sigmaLow))#/(1-pnorm(crit, alts[3], noise.sigmaLow))
    prSoA_alt4[i] = (pnorm(confBounds_NotMe[(nRNotMe+2)-i], alts[4], noise.sigmaLow) - pnorm(confBounds_NotMe[(nRNotMe+1)-i], alts[4], noise.sigmaLow))#/(1-pnorm(crit, alts[4], noise.sigmaLow))
  }
  for (i in (nRNotMe+1):6) {
    prSoA_alt1[i] = (pnorm(confBounds_Me[i-nRNotMe], alts[1], noise.sigmaLow) - pnorm(confBounds_Me[i-(nRNotMe-1)], alts[1], noise.sigmaLow))#/pnorm(crit, alts[1], noise.sigmaLow)
    prSoA_alt2[i] = (pnorm(confBounds_Me[i-nRNotMe], alts[2], noise.sigmaLow) - pnorm(confBounds_Me[i-(nRNotMe-1)], alts[2], noise.sigmaLow))#/pnorm(crit, alts[2], noise.sigmaLow)
    prSoA_alt3[i] = (pnorm(confBounds_Me[i-nRNotMe], alts[3], noise.sigmaLow) - pnorm(confBounds_Me[i-(nRNotMe-1)], alts[3], noise.sigmaLow))#/pnorm(crit, alts[3], noise.sigmaLow)
    prSoA_alt4[i] = (pnorm(confBounds_Me[i-nRNotMe], alts[4], noise.sigmaLow) - pnorm(confBounds_Me[i-(nRNotMe-1)], alts[4], noise.sigmaLow))#/pnorm(crit, alts[4], noise.sigmaLow)
  }
  
  logL = 0
  for (i in 1:6) {
    soaData = dplyr::filter(data_LowNoise, SoARatings == i)
    for (alt in 1:4){
      if (alt==1){
        pr = prSoA_alt1
      } else if (alt==2){
        pr = prSoA_alt2
      }else if (alt==3){
        pr = prSoA_alt3
      }else{
        pr = prSoA_alt4
      }
      if (length(soaData$count[which(soaData$Alterations==alts[alt])]) > 0) {
        logL = logL + (soaData$count[which(soaData$Alterations==alts[alt])])*log(pr[i])
      }
    }
  }
  # High Noise
  
  confBounds_NotMe_HighNoise = crit
  confBounds_Me_HighNoise = crit
  buffer = bufferFactor*noise.sigmaHigh
  a = alts[4] + buffer
  b = alts[1] - buffer
  adist = abs(a-crit)
  bdist = abs(b-crit)
  dist = max(c(adist, bdist))
  notMedist = dist/nRNotMe
  medist = dist/nRMe
  
  nextBound = crit
  for (i in 2:nRNotMe){
    nextBound = nextBound + notMedist
    confBounds_NotMe_HighNoise = c(confBounds_NotMe_HighNoise, nextBound)
  }
  confBounds_NotMe_HighNoise = c(confBounds_NotMe_HighNoise, Inf)
  nextBound = crit
  for (i in 2:nRMe){
    nextBound = nextBound - medist
    confBounds_Me_HighNoise = c(confBounds_Me_HighNoise, nextBound)
  }
  confBounds_Me_HighNoise = c(confBounds_Me_HighNoise, -Inf)
  
  # Calculate probabilities of ratings
  prSoA_alt1_HighNoise = c(0,0,0,0,0,0)
  prSoA_alt2_HighNoise = c(0,0,0,0,0,0)
  prSoA_alt3_HighNoise = c(0,0,0,0,0,0)
  prSoA_alt4_HighNoise = c(0,0,0,0,0,0)
  
  for (i in 1:nRNotMe){
    prSoA_alt1_HighNoise[i] = (pnorm(confBounds_NotMe_HighNoise[(nRNotMe+2)-i], alts[1], noise.sigmaHigh) - pnorm(confBounds_NotMe_HighNoise[(nRNotMe+1)-i], alts[1], noise.sigmaHigh))#/(1-pnorm(crit, alts[1], noise.sigmaHigh))
    prSoA_alt2_HighNoise[i] = (pnorm(confBounds_NotMe_HighNoise[(nRNotMe+2)-i], alts[2], noise.sigmaHigh) - pnorm(confBounds_NotMe_HighNoise[(nRNotMe+1)-i], alts[2], noise.sigmaHigh))#/(1-pnorm(crit, alts[2], noise.sigmaHigh))
    prSoA_alt3_HighNoise[i] = (pnorm(confBounds_NotMe_HighNoise[(nRNotMe+2)-i], alts[3], noise.sigmaHigh) - pnorm(confBounds_NotMe_HighNoise[(nRNotMe+1)-i], alts[3], noise.sigmaHigh))#/(1-pnorm(crit, alts[3], noise.sigmaHigh))
    prSoA_alt4_HighNoise[i] = (pnorm(confBounds_NotMe_HighNoise[(nRNotMe+2)-i], alts[4], noise.sigmaHigh) - pnorm(confBounds_NotMe_HighNoise[(nRNotMe+1)-i], alts[4], noise.sigmaHigh))#/(1-pnorm(crit, alts[4], noise.sigmaHigh))
  }
  for (i in (nRNotMe+1):6) {
    prSoA_alt1_HighNoise[i] = (pnorm(confBounds_Me_HighNoise[i-nRNotMe], alts[1], noise.sigmaHigh) - pnorm(confBounds_Me_HighNoise[i-(nRNotMe-1)], alts[1], noise.sigmaHigh))#/pnorm(crit, alts[1], noise.sigmaHigh)
    prSoA_alt2_HighNoise[i] = (pnorm(confBounds_Me_HighNoise[i-nRNotMe], alts[2], noise.sigmaHigh) - pnorm(confBounds_Me_HighNoise[i-(nRNotMe-1)], alts[2], noise.sigmaHigh))#/pnorm(crit, alts[2], noise.sigmaHigh)
    prSoA_alt3_HighNoise[i] = (pnorm(confBounds_Me_HighNoise[i-nRNotMe], alts[3], noise.sigmaHigh) - pnorm(confBounds_Me_HighNoise[i-(nRNotMe-1)], alts[3], noise.sigmaHigh))#/pnorm(crit, alts[3], noise.sigmaHigh)
    prSoA_alt4_HighNoise[i] = (pnorm(confBounds_Me_HighNoise[i-nRNotMe], alts[4], noise.sigmaHigh) - pnorm(confBounds_Me_HighNoise[i-(nRNotMe-1)], alts[4], noise.sigmaHigh))#/pnorm(crit, alts[4], noise.sigmaHigh)
  }
  
  for (i in 1:6) {
    soaData = dplyr::filter(data_HighNoise, SoARatings == i)
    for (alt in 1:4){
      if (alt==1){
        pr = prSoA_alt1_HighNoise
      } else if (alt==2){
        pr = prSoA_alt2_HighNoise
      }else if (alt==3){
        pr = prSoA_alt3_HighNoise
      }else{
        pr = prSoA_alt4_HighNoise
      }
      if (length(soaData$count[which(soaData$Alterations==alts[alt])]) > 0) {
        logL = logL + (soaData$count[which(soaData$Alterations==alts[alt])])*log(pr[i])
      }
    }
  }
  
  if (is.nan(logL)) {
    nll = -Inf
  } else {
    nll = -logL
  }
  
  return(nll)  
  
}