thesis/VectorizedSlopeModel.stan.R at master · crowding/thesis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
splits <- c("subject", "content",
            "displacement",
            "target_number_all",
            "target_number_shown",
            "spacing", "eccentricity", "bias")

relevant <- splits %v% c("n_cw", "n_obs")

lapse_limit <- 0.05

stan_format <- mkchain(
    subset(select=relevant),
    as.list,
    factorify,
    put(names(.), gsub('\\.', '_', names(.))),
    within({
      N <- length(subject_ix)
      lapse_limit <- lapse_limit
    }))

model_code <- '
data {
  int<lower=0> N;
  vector[N] target_number_all;
  vector[N] target_number_shown;
  vector[N] displacement;
  vector[N] content;
  vector[N] spacing;
  vector[N] eccentricity;
  int n_cw[N];
  int n_obs[N];
  real lapse_limit;
}

transformed data {
  vector[N] frac_spacing;
  vector[N] ones;
  vector[N] content_nonlinear;
  frac_spacing <- (2 * pi() * target_number_all);
  for (n in 1:N) {
    ones[n] <- 1.0;
    content_nonlinear[n] <- content[n] * fabs(content[n]);
  }
}

parameters {
  real beta_dx;
  real<lower=0, upper=lapse_limit> lapse;
  real intercept;

  real<lower=0,upper=2*pi()> cs;
  real summation;
  real repulsion;
  real nonlinearity;
}

model {
  vector[N] crowdedness;
  vector[N] link_displacement;
  vector[N] link_repulsion;
  vector[N] link_summation;
  vector[N] link;
  vector[N] response;

  crowdedness <- 2 - (ones*2) ./ (1+exp((ones*cs) ./ frac_spacing));
  link_displacement <- beta_dx * displacement .* crowdedness;
  link_repulsion <- (repulsion * content + nonlinearity * content_nonlinear);
  link_summation <- summation * content .* target_number_shown;
  link <- intercept + link_displacement + link_repulsion + link_summation;
  response <- (ones ./ (1+ exp(-link))) * (1-lapse) + lapse/2;
  n_cw ~ binomial( n_obs, response );
}'

 stan_predict <- mkchain[., coefs](
     mutate(frac_spacing = 2*pi/target_number_all)
   , with(coefs, summarize(
       .
     , link_displacement = (beta_dx * displacement
                            * (2 - 2/(1+exp(-cs/frac_spacing))))
     , link_repulsion = (repulsion * content
                         + nonlinearity * (content * abs(content)))
     , link_summation = (target_number_all * content * summation)
     , link = intercept + link_displacement + link_repulsion + link_summation
     , response = plogis(link) * (1-lapse) + lapse/2
     )))