Как использовать brms (пакет R) для генерации кода Stan, который мне нужен для воспроизведения оценки моделей в pystan? - PullRequest
Новая
       136

Как использовать brms (пакет R) для генерации кода Stan, который мне нужен для воспроизведения оценки моделей в pystan?

0 голосов
/ 12 марта 2020

Я разработал конвейер для оценки модели с использованием пакета brms R, и теперь мне нужно преобразовать его в python. Я понимаю, что ближе всего к brms в python я могу найти pystan, где я должен написать свою модель с использованием синтаксиса Stan. Я хотел бы знать, есть ли функция brms, которая генерирует код Stan, который можно использовать в качестве аргумента model_code для функции pystan.StanModel в python. Я попытался использовать код, сгенерированный функцией make_stancode, но он не сработал.

Это код, сгенерированный make_stancode: 

life_span_code = """
// generated with brms 2.10.0
functions {

  /* compute monotonic effects
   * Args:
   *   scale: a simplex parameter
   *   i: index to sum over the simplex
   * Returns:
   *   a scalar between 0 and 1
   */
  real mo(vector scale, int i) {
    if (i == 0) {
      return 0;
    } else {
      return rows(scale) * sum(scale[1:i]);
    }
  }
}
data {
  int<lower=1> N;  // number of observations
  vector[N] Y;  // response variable
  int<lower=1> Ksp;  // number of special effects terms
  int<lower=1> Imo;  // number of monotonic variables
  int<lower=2> Jmo[Imo];  // length of simplexes
  // monotonic variables
  int Xmo_1[N];
  // prior concentration of monotonic simplexes
  vector[Jmo[1]] con_simo_1;
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  // temporary intercept for centered predictors
  real Intercept;
  // special effects coefficients
  vector[Ksp] bsp;
  // simplexes of monotonic effects
  simplex[Jmo[1]] simo_1;
  real<lower=0> sigma;  // residual SD
}
transformed parameters {
}
model {
  // initialize linear predictor term
  vector[N] mu = Intercept + rep_vector(0, N);
  for (n in 1:N) {
    // add more terms to the linear predictor
    mu[n] += (bsp[1]) * mo(simo_1, Xmo_1[n]);
  }
  // priors including all constants
  target += student_t_lpdf(Intercept | 3, 65, 12);
  target += dirichlet_lpdf(simo_1 | con_simo_1);
  target += student_t_lpdf(sigma | 3, 0, 12)
    - 1 * student_t_lccdf(0 | 3, 0, 12);
  // likelihood including all constants
  if (!prior_only) {
    target += normal_lpdf(Y | mu, sigma);
  }
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
}
"""

Это код, который я использую в python: 

## Libraries
import pandas as pd
import pystan
import numpy as np
import random as rd

## Build data for life span example with ordenated factors

income_options =  ["below_20", "20_to_40", "40_to_100", "greater_100"]
income_mean = [30, 60, 70, 75]
income_factor = [0, 1, 2, 3]

dict_data = {'income_options' : income_options,
             'income_mean' : income_mean,
             'income_factor' :  income_factor}

map_df = pd.DataFrame(dict_data)

income_rep = rd.sample(income_factor*25, 100)

rand_inc = np.random.normal(loc = 0, scale = 1, size = 100).tolist()


data_df = pd.DataFrame({'income_factor': income_rep,
                        'rand_inc' : rand_inc})

data_df = pd.merge(data_df, map_df, on = 'income_factor')

data_df['ls'] = data_df['income_mean'] + data_df['rand_inc']

N = data_df.shape[0]
Y = data_df['ls'].tolist()
K = 1
X = [1]*N
Ksp = 1
Imo = 1
Xmo_1 = data_df['income_factor'].tolist()
Jmo = len(data_df['income_factor'].unique().tolist())-1
con_simo_1 = [1]*Jmo
prior_only = 0


life_span_data = {'N' : N,
                  'Y' : Y,
                  'K' : K,
                  'X' : X,
                  'Ksp' : Ksp,
                  'Imo' : Imo,
                  'Xmo_1' : Xmo_1,
                  'Jmo' : Jmo,
                  'con_simo_1' : con_simo_1,
                  'prior_only' : prior_only}

life_span_sm = pystan.StanModel(model_code = life_span_code)
life_span_fit = life_span_sm.sampling(data= life_span_data, iter=1000, chains=2)

И эта ошибка, которую я получаю:

"RuntimeError: Исключение: несоответствие в измерениях чисел, объявленных и найденных в контексте; этап обработки = инициализация данных; имя переменной = Jmo; объявлено dims = (1); найдено dims = () (в «неизвестное имя файла» в строке 24) "

Спасибо за помощь

1 Ответ

0 голосов
/ 14 апреля 2020

Оказывается, проблема была не в коде модели, созданном brms, а в том, как я определил аргументы. В частности, Jmo должен быть списком вместо int. 

N = data_df.shape[0]
Y = data_df['ls'].tolist()
K = 1
X = [1]*N
Ksp = 1
Imo = 1
Xmo_1 = data_df['income_factor'].tolist()

## The following two lines have changed
Jmo = [len(data_df['income_factor'].unique().tolist())-1]
con_simo_1 = [1, 1, 1]
## End of changes

prior_only = 0

Остальная часть кода такая же. Я все еще хотел бы получить разъяснения о том, почему некоторые аргументы могут быть объявлены как целые, а другие только как списки.

Еще раз спасибо

...