solving knapsack problem

model.stan

functions {
  real max_value(int I, int W, vector value, int[] weight) {
    real dp[I+1,W+1];
    for (w in 0:W) dp[1,w+1] = 0;
    for (i in 1:I) {
      for (w in 0:W) {
        if (w < weight[i]) {
          dp[i+1,w+1] = dp[i,w+1];
        } else {
          dp[i+1,w+1] = fmax(dp[i,w+1], dp[i,w-weight[i]+1] + value[i]);
        }
      }
    }
    return dp[I+1,W+1];
  }
}

data {
  int I;
  int I_obs;
  int I_mis;
  int Idx_obs[I_obs];
  int Idx_mis[I_mis];
  int Weight[I];
  vector[I_obs] V_obs;
  int Max_W;
  int N;
  vector[N] Y;
}

parameters {
  vector<lower=0>[I_mis] v_mis;
  real<lower=0> s_Y;
}

transformed parameters {
  real mu;
  {
    vector[I] v;
    for (i in 1:I_obs) v[Idx_obs[i]] = V_obs[i];
    for (i in 1:I_mis) v[Idx_mis[i]] = v_mis[i];
    mu = max_value(I, Max_W, v, Weight);
  }
}

model {
  v_mis ~ normal(3, 5);
  Y ~ normal(mu, s_Y);
}

run-model.R

library(rstan)

stanmodel <- stan_model(file='model.stan')

d <- read.csv(file='z-data-items.csv')
d_obs <- na.omit(d)
I <- nrow(d)
I_obs <- nrow(d_obs)
I_mis <- I - I_obs
Idx_obs <- d_obs$i
Idx_mis <- which(is.na(d$Value))

Max_W <- 40
Y <- c(63.1, 66.8, 72.5, 58.1, 55.2, 67.5, 69.1, 70.2)
N <- length(Y)

data <- list(I=I, I_obs=I_obs, I_mis=I_mis, Idx_obs=Idx_obs, Idx_mis=Idx_mis,
  Weight=d$Weight, V_obs=d_obs$Value, Max_W=Max_W, N=length(Y), Y=Y)
fit <- sampling(stanmodel, data=data, seed=1234, control=list(adapt_delta=0.99))

z-data-items.csv