『Pythonではじめる数理最適化』の7章「商品推薦のための興味のスコアリング」をStanで解く

model.stan

data {
  int I; // number of data
  int R; // number of Rcen
  int F; // number of Freq
  array[I] int Rcen; // value of Rcen
  array[I] int Freq; // value of Freq
  array[I] int N;
  array[I] int PV;
}

parameters {
  matrix<lower=0>[R, F] dx;
}

transformed parameters {
  matrix[R, F] x;
  matrix[R, F] q;
  x[1,1] = 5 - dx[1,1];
  for (r in 2:R) {
    x[r,1] = x[r-1,1] - dx[r,1];
  }
  for (f in 2:F) {
    x[1,f] = x[1,f-1] - dx[1,f];
  }
  for (r in 2:R) {
    for (f in 2:F) {
      x[r,f] = min([x[r-1,f], x[r,f-1]]) - dx[r,f];
    }
  }
  q[1:R,1:F] = inv_logit(x[1:R,1:F]);
}

model {
  for (i in 1:I) {
    PV[i] ~ binomial(N[i], q[Rcen[i], Freq[i]]);
  }
}

plot-result.py

import pandas
import cmdstanpy
import numpy as np
import matplotlib.pyplot as plt
from plotnine import ggplot, aes, theme, element_text, geom_ribbon, geom_line, geom_point, labs

rf_df = pandas.read_csv('input/rf_df.csv')
fit = cmdstanpy.from_csv('output/result')
q_ms = fit.stan_variable(var='q')

Freq = rf_df.freq.unique().tolist()
Rcen = rf_df.rcen.unique().tolist()
Z = [np.median(q_ms[:, rcen-1, 7-freq]) for freq in Freq for rcen in Rcen]
Z = np.array(Z).reshape((len(Freq), len(Rcen)))
X, Y = np.meshgrid(Rcen, Freq)  
fig = plt.figure(dpi=250)
ax = fig.add_subplot(111, projection='3d', xlabel='rcen', ylabel='freq', zlabel='pred_prob')
ax.plot_wireframe(X, Y, Z)
plt.savefig('output/q_median.png') 

plot_f = 5
qua = np.quantile(q_ms[:, :, 7-plot_f], [0.025, 0.25, 0.50, 0.75, 0.975], axis=0)
d_est = pandas.DataFrame(np.column_stack([np.arange(1,8), qua.T]), \
    columns=['rcen', '2.5%', '25%', '50%', '75%', '97.5%'])
rf_df_at_plot_f = rf_df[rf_df.freq == plot_f]
p = (ggplot()
    + theme(text=element_text(size=18))
    + geom_ribbon(d_est, aes(x='rcen', ymin='2.5%', ymax='97.5%'), fill='blue', alpha=1/6)
    + geom_ribbon(d_est, aes(x='rcen', ymin='25%', ymax='75%'), fill='blue', alpha=2/6)
    + geom_line(d_est, aes(x='rcen', y='50%'), color='blue', size=1)
    + geom_point(rf_df_at_plot_f, aes(x='rcen', y='prob'), size=1)
    + labs(y='prob')
) 
p.save(filename=f'output/q_at_f_{plot_f}.png', dpi=300, width=5, height=4)

rf_df.csv

run-model.py

import pandas
import cmdstanpy

d = pandas.read_csv('input/rf_df.csv')
data = {'I':len(d), 'R':7, 'F':7, 
        'Rcen':d.rcen, 'Freq':8 - d.freq, 'N':d.N, 'PV':d.pv}
model = cmdstanpy.CmdStanModel(stan_file='model/model.stan')
fit = model.sample(data=data, seed=123)

fit.save_csvfiles('output/result')
fit.summary().to_csv('output/fit-summary.csv')