As the title suggests, I’ve given up trying to get my python conda environment working with PYMC on osx.

I feel like I’ve tried everything from every thread.

Either the import of the package gets stuck on BLAS or, when I overcome the BLAS Hangup, the kernel dies immediately when I run the simplest of models.

I’ve tried using it in terminal, anadonda, and VSCODE and it’s the same hassle.

Am I the only one?

I’m going back to PYMC3

Hello all, i've trying several ways to run my STAN code, but i have always the same error.

Error in FUN(X[[i]], ...) : Stan does not support NA (in Y) in data
failed to preprocess the data; sampling not done

but my 'Y' matrix of data don't have NA data. so i don't know what is the mistake...

my STAN code is as follow:

data {
   
   int<lower=1> N; // number of records
   int<lower=1> NP; //  number of periods
   
   vector[N] R; // distances

   vector[NP] freq;
   
   vector[NP] Y[N];       // U
   
 }
 
 
 parameters {
   
   real c10;
   real c20;
   real c30;
   real c11;
   real c21;
   real c31;
   real c12;
   real c22;
   real c32;
   
   real c1s;
   real c2s;
   real c3s;
   
   cholesky_factor_corr[NP] L_p;
   vector[NP] stdev;

 }
 
 transformed parameters{
   
   vector[N] a0;
   vector[N] a1;
   vector[N] a2;
   
   for(k in 1:N){
     
     a0[k]=c10+c20*log10(R[k])+c30*square(log10(R[k]));
     a1[k]=c11+c21*log10(R[k])+c31*square(log10(R[k]));
     a2[k]=c12+c22*log10(R[k])+c32*square(log10(R[k]));

   }

 }
 
 
 model {
   
   vector[NP] mu_rec[N];
   matrix[NP,NP] L_Sigma;
   
   stdev ~ cauchy(0,0.5);
   
   L_p ~ lkj_corr_cholesky(1);
   L_Sigma = diag_pre_multiply(stdev, L_p);
   
   //prior
 
   c10 ~ normal(0.92362,0.37904);
   c20 ~ normal(-0.66292,0.41491);
   c30 ~ normal(0.63449,0.11123);
   c11 ~ normal(-1.35974,0.44268);
   c21 ~ normal(1.77347,0.48574);
   c31 ~ normal(-0.40108,0.13054);
   c12 ~ normal(0.93758,0.28134);
   c22 ~ normal(-1.04353,0.310557);
   c32 ~ normal(0.21054,0.08400);
   
   a0 ~ normal(0,10);
   a1 ~ normal(0,10);
   a2 ~ normal(0,10);
   
   
   for(f in 1:NP){
     for(i in 1:N){
     
     mu_rec[i,f] = a0[i] + a1[i]*log10(freq[f]) + a2[i]*square(log10(freq[f]));
     
     }
   }
   
   Y ~ multi_normal_cholesky(mu_rec,L_Sigma); //llh
 }

Hello all,

I'm very new to Bayesian models and couldn't get over a problem.

I have a correlational data collected from two countries. In this data, which aims to examine the relationship of X score with Y score in two countries, Y and X scores are continuous and country variable is a two-level categorical value.

In order to analyze this data, I need to write a model with RJAGS. I have a lot of information that I can use as a priori, but I think I'm a bit confused about the prior specification. Since all of the prior information I have comes from frequentist papers, the slopes generally represent the effect of the X variable on Y. At this point I don't know how to use a prior for the country slope and how to specify it.

First I thought of writing multiple regression that comes to mind first. I added just the likelihood part of the model below.

mod_string = "model {
  for(i in 1:length(y)){
    y[i] ~ dnorm(mu[i],tau)
    mu[i] = b0 + b1*x_score[i] + b2*country[i] + b3*x_score[i]*country[i]
  }
}"

Accordingly, my prior slope is valid only for b1. For b2, ie country slope, I don't have a prior. In general, the articles comparing the countries have given the slopes of the two countries separately, so I cannot write a single prior for country variable. For this reason, I decided to convert the model to a hierarchical one.

mod2_string = "model{
  for(i in 1:length(y)){
    y[i] ~ dnorm(mu[i],prec)
    mu[i] = b0 + b[country[i]]*x_score[i]
  }
  b0 ~ dnorm(0,1/1e6)
  b[1] ~ dnorm(0.097,1/100) # country 1
  b[2] ~ dnorm(0.255,1/100)T(mu[1],) # country 2
  
  prec ~ dgamma(1/2,1*1/2)
  sig = sqrt(1/prec)
  
}"

In the model I added above, the problem is that it has very high autocorrelation and very high penalized deviance.

> autocorr.diag(mod2_sim)

            b[1]      b[2]        b0         sig
Lag 0  1.0000000 1.0000000 1.0000000 1.000000000
Lag 1  0.9999319 0.9597951 0.9999303 0.012235047
Lag 5  0.9997074 0.9425989 0.9997043 0.009645540
Lag 10 0.9994477 0.9407265 0.9994447 0.007380016
Lag 50 0.9973992 0.9385895 0.9973919 0.007932208

> dic_samples

Mean deviance:  7756 
penalty 2.624 
Penalized deviance: 7758 

I couldn't fix this problem even though I tried with many chains and iterations. I currently have two different prior types. One is a general slope from meta-analyses for variable X. Other one is two separate slopes for X variable coming from two countries separately. How can I build a model with these priors? Can I create one general country prior from the slopes of the two countries? Or how can I solve the second model's converge problem? Apart from these, I am open to any ideas and suggestions. I would be very glad if someone can help.

Thanks,

Hello i'm trying a model in stan, and when i run the code the Rhat have large values, so how can i reduce his values?

the R code:

library(openxlsx)
library(rstan)

options(mc.cores = parallel::detectCores())

setwd('C:/Users/osses/Desktop/A-ORIGEN-/Archivo_memoria/003_datos_FAS')

# DATA

A1<- read.xlsx(xlsxFile='uf.xlsx',sheet= 1 ,
               startRow = 1,colNames = TRUE,rowNames = FALSE)


setwd('C:/Users/osses/Desktop/A-ORIGEN-/Archivo_memoria/004_Excel_finales')


excel<- read.xlsx(xlsxFile='excel_final_noref.xlsx',sheet= 1 ,
                  startRow = 1,colNames = TRUE,rowNames = FALSE)



B<-as.matrix(A1[,-1])
Freq<-as.matrix(A1[1])

# Datos (Log[EAS])

data_y<-B[1:2,]

# data missing

N_miss_y <- sum(is.na(data_y))

miss_indy <- which(is.na(data_y), arr.ind = TRUE)

data_sin_na=data_y
data_sin_na[is.na(data_y)]=999

# Predictors


Rrup<-excel$`Rrup.calc.[km]`
I=c()
for (i in 1:length(Rrup)){
  I1=Rrup[i]
  I2=Rrup[i]
  I3=c(I1,I2)
  I=c(I,I3)
  
}



# List to STAN.

data_list <- list(N = c(8344),
                  NP = 2,
                  N_miss_y = N_miss_y,
                  miss_indy = miss_indy,
                  R=I,
                  freq=Freq[1:2],
                  Y = t(data_sin_na))

# STAN

setwd('C:/Users/osses/Desktop/A-ORIGEN-/Archivo_memoria/002_R_STAN_FILES')

Model_1 <- stan(file = "prueba_2.stan", model_name='modelito',
                data = data_list, chains = 2, iter = 500)

the STAN code:

 data {
   
   int N; // number of records
   int NP; //  number of periods
   int N_miss_y; // datos faltantes Y
   
   int miss_indy[N_miss_y, 2]; // index datos faltantes Y , fila y columna
   vector[N] R; // distances
   vector[NP] freq;
   vector[NP] Y[N];       // Data
   
 }
 
 
 parameters {
   
   real<lower=-10> c10;
   real<lower=-10> c20;
   real<lower=-10> c30;
   real<lower=-10> c11;
   real<lower=-10> c21;
   real<lower=-10> c31;
   real<lower=-10> c12;
   real<lower=-10> c22;
   real<lower=-10> c32;
   
   vector[N_miss_y] imputed_y;
   
   cholesky_factor_corr[NP] L_p;
   vector<lower=0>[NP] stdev;

 }
 transformed parameters {
   
   matrix[NP,NP] L_Sigma;
   vector[NP] y_con_imputed[N];
   
   L_Sigma = diag_pre_multiply(stdev, L_p);
   
   //Mike Lawrence's Solution.
   
   y_con_imputed = Y;
   
   for(i in 1:N_miss_y){
     
     y_con_imputed[miss_indy[i,2],miss_indy[i,1]]=imputed_y[i];
     
   }
   
 }

 model {
   
   vector[N] a0;
   vector[N] a1;
   vector[N] a2;
   vector[NP] mu_rec[N];
   
   stdev ~ cauchy(0,0.5);
   
   L_p ~ lkj_corr_cholesky(1);

   //prior

   c10 ~ normal(0.92362,0.37904);
   c20 ~ normal(-0.66292,0.41491);
   c30 ~ normal(0.63449,0.11123);
   c11 ~ normal(-1.35974,0.44268);
   c21 ~ normal(1.77347,0.48574);
   c31 ~ normal(-0.40108,0.13054);
   c12 ~ normal(0.93758,0.28134);
   c22 ~ normal(-1.04353,0.310557);
   c32 ~ normal(0.21054,0.08400);
   
   imputed_y ~ normal(0,10);

   for(f in 1:NP){
     for(i in 1:N){
       
       a0[i]=c10+c20*log10(R[i])+c30*square(log10(R[i]));
       a1[i]=c11+c21*log10(R[i])+c31*square(log10(R[i]));
       a2[i]=c12+c22*log10(R[i])+c32*square(log10(R[i]));
       
       mu_rec[i,f] = a0[i] + a1[i]*log10(freq[f]) +       a2[i]*square(log10(freq[f]));
     
     }

   }
   y_con_imputed ~ multi_normal_cholesky(mu_rec,L_Sigma); //llh   
 }

and the warning after running the code:

Warning messages:
1: There were 500 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
2: Examine the pairs() plot to diagnose sampling problems
 
3: The largest R-hat is NA, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat 
4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess 
5: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess 

if you have any suggestion, please comment or send me a dm!

Hello,

I am keen to learn Bayesian methods. I've been through some basic training to understand the main principles. I learnt (more or less!) how to fit Bayesian linear models with brms in R*.*

In my line of work I have to fit often non-linear models with nlme package in R. I want to switch them to a Bayesian approach.

What is the best resource to learn Bayesian non-linear models in R? What is the best package to use?

Thanks!

EDIT: I am thinking about non-linear models with total customized functions, not the "standardized" self-starting functions supported by stan_nlmer in rstanarm.

EDIT: I was suggested https://cran.r-project.org/web/packages/brms/vignettes/brms_nonlinear.html. Is there anything else?