Hi all,
im triyng to implement a bayesian model with R and c++.
I have a strange problem. I can't reproduce the error with a small
script and then i post the original one.
The problem is after the line
for(MCMC_iter2=0;MCMC_iter2<thin;MCMC_iter2++)
For the first 34 iterations all work fine, after, all the simulations
of mu_acc_P return an "nan". If i delete the line
alpha_acc_P[0] = rnorm(app1_3,sqrt(app1_4) );
the simulations of mu_acc_P work fine. For me it is really strange,
Thanks!.
#include <iostream>
#include <string>
using namespace std;
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Linpack.h>
#include <R_ext/Lapack.h>
#include <R_ext/BLAS.h>
#include "funzioni.h"
#define _USE_MATH_DEFINES
extern "C" {
SEXP Model1_imp11Cpp(
SEXP ad_start_r, SEXP ad_end_r, SEXP ad_esp_r,
SEXP burnin_r, SEXP thin_r,SEXP iter_1_r,SEXP iter_2_r,
SEXP n_j_r,SEXP J_r,
SEXP prior_alpha_r,SEXP prior_rho_r,SEXP prior_sigma2_r,SEXP
prior_phi_r,SEXP prior_zeta2_r,
SEXP sdprop_rho_r,SEXP sdprop_sigma2_r,
SEXP start_alpha_r,SEXP start_mu_r,SEXP start_rho_r,SEXP
start_sigma2_r,SEXP start_k_r,SEXP start_phi_r, SEXP start_zeta2_r,
SEXP x_r,SEXP H_r, SEXP acceptratio_r
// SEXP Sigma_adapt_sp_r, SEXP mean_adapt_sp_r, SEXP sim_acc_sp_r,
SEXP ep_sp_r, SEXP acc_index_sp_r,
// SEXP Sigma_adapt_k_r, SEXP mean_adapt_k_r, SEXP sim_acc_k_r, SEXP
ep_k_r, SEXP acc_index_k_r
//
){
/*****************************************
Varie ed eventuali
*****************************************/
// indici
int i,j,k,l,h,t,f,info,MCMC_iter,MCMC_iter2;
int nProtect= 0;
int indice_adapt=0;
double duepi = (2*M_PI);
// costanti
char const *lower = "L";
char const *upper = "U";
char const *ntran = "N";
char const *ytran = "T";
char const *rside = "R";
char const *lside = "L";
const double one = 1.0;
const double negOne = -1.0;
const double zero = 0.0;
const int incOne = 1;
/*****************************************
Set-up
*****************************************/
double *x_P = REAL(x_r);
int n_j = INTEGER(n_j_r)[0];
int J = INTEGER(J_r)[0];
int N = n_j*J;
int iter_1 = INTEGER(iter_1_r)[0];
int iter_2 = INTEGER(iter_2_r)[0];
int burnin = INTEGER(burnin_r)[0];
int thin = INTEGER(thin_r)[0];
int ad_start = INTEGER(ad_start_r)[0];
int ad_end = INTEGER(ad_end_r)[0];
double *ad_esp_P = REAL(ad_esp_r);
double *acceptratio_P = REAL(acceptratio_r);
double *H_P = REAL(H_r);
// int nSamples_save = INTEGER(nSamples_save_r)[0];
// PRIOR
double *prior_alpha = REAL(prior_alpha_r);
double M_alpha = prior_alpha[0]; double sigma2_alpha =
prior_alpha[1];
double *prior_rho = REAL(prior_rho_r);
double a_rho = prior_rho[0]; double b_rho = prior_rho[1];
double *prior_sigma2 = REAL(prior_sigma2_r);
double a_sigma2 = prior_sigma2[0]; double b_sigma2 = prior_sigma2[1];
double *prior_zeta2 = REAL(prior_zeta2_r);
double a_zeta2 = prior_zeta2[0]; double b_zeta2 = prior_zeta2[1];
double *prior_phi = REAL(prior_phi_r);
double M_phi= prior_phi[0]; double sigma2_phi = prior_phi[1];
double lim_down= prior_phi[2]; double lim_up = prior_phi[3];
// PROPOSAL
double *sdprop_rho = REAL(sdprop_rho_r);
double sd_rho = sdprop_rho[0];
double *sdprop_sigma2 = REAL(sdprop_sigma2_r);
double sd_sigma2 = sdprop_sigma2[0];
// STARTING
double *start_alpha_P = REAL(start_alpha_r);
double start_alpha = start_alpha_P[0];
double *start_mu_P = REAL(start_mu_r);
// for(j=0;j<J;j++)
// {
// start_mu[j] = start_mu_P+j) ;
// }
double *start_phi_P = REAL(start_phi_r);
double start_phi = start_phi_P[0];
double *start_rho_P = REAL(start_rho_r);
double start_rho = start_rho_P[0];
double *start_sigma2_P = REAL(start_sigma2_r);
double start_sigma2 = start_sigma2_P[0];
double *start_zeta2_P = REAL(start_zeta2_r);
double start_zeta2 = start_zeta2_P[0];
double *start_k_P = REAL(start_k_r);
double start_k[N-1];
for(j=0;j<N;j++)
{
start_k[j] = *(start_k_P+j) ;
}
/*****************************************
Varie ed eventuali II
*****************************************/
const int incJ = J;
const int incN = N;
const int nn_j = n_j*n_j;
const int inc4 = 4;
const int inc2 = 2;
double app1_1, app1_2, app1_3, app1_4;
double appn_j_1[n_j];
double app_uno_n_j[n_j];
for(i=0;i<n_j;i++)
{
app_uno_n_j[i] = 1;
}
int nSamples_save = iter_2;
// generatore random
GetRNGstate();
// /*****************************************
// // UPDATE PARAMETRI
// *****************************************/
// alpha
double *alpha_acc_P = (double *) R_alloc(1, sizeof(double));
F77_NAME(dcopy)(&incOne, start_alpha_P, &incOne, alpha_acc_P, &incOne);
// mu_acc
double *mu_acc_P = (double *) R_alloc(J, sizeof(double));
F77_NAME(dcopy)(&incJ, start_mu_P, &incOne, mu_acc_P, &incOne);
// phi
double *phi_acc_P = (double *) R_alloc(1, sizeof(double));
F77_NAME(dcopy)(&incOne, start_phi_P, &incOne, phi_acc_P, &incOne);
// zeta2
double *zeta2_acc_P = (double *) R_alloc(1, sizeof(double));
F77_NAME(dcopy)(&incOne, start_zeta2_P, &incOne, zeta2_acc_P, &incOne);
// sigma2
double *sigma2_acc_P = (double *) R_alloc(1, sizeof(double));
F77_NAME(dcopy)(&incOne, start_sigma2_P, &incOne, sigma2_acc_P,
&incOne);
// rho
double *rho_acc_P = (double *) R_alloc(1, sizeof(double));
F77_NAME(dcopy)(&incOne, start_rho_P, &incOne, rho_acc_P, &incOne);
// k
double *k_acc_P = (double *) R_alloc(n_j*J, sizeof(double));
F77_NAME(dcopy)(&incN, start_k_P, &incOne, k_acc_P, &incOne);
// variabili di appoggio per i parametri
double *mu_app_P = (double *) R_alloc(J, sizeof(double));
F77_NAME(dcopy)(&incJ, start_mu_P, &incOne, mu_app_P, &incOne);
/*****************************************
// OUTPUT di R
*****************************************/
SEXP mu_out_r,alpha_out_r;
PROTECT(mu_out_r = allocMatrix(REALSXP, J, nSamples_save));
nProtect++;
double *mu_out_P = REAL(mu_out_r);
PROTECT(alpha_out_r = allocMatrix(REALSXP, 1, nSamples_save));
nProtect++;
double *alpha_out_P = REAL(alpha_out_r);
//
// PROTECT(accept_r = allocMatrix(REALSXP, 1, 1)); nProtect++;
/*****************************************
Varie ed eventuali III
*****************************************/
// MATRICI COVARIANZA E AFFINI
double *C_P = (double *) R_alloc(nn_j, sizeof(double));
double *C_cand_P = (double *) R_alloc(nn_j, sizeof(double));
double Sum_Sigma_inv;
double Sigma_inv_one[n_j-1];
double logdet;
//
// /*****************************************
// // PARAMETRI PARTE ADATTIVA
// *****************************************/
// /*****************************************
// // STARTING MCMC
// *****************************************/
// MATRICE DI COVARIANZA
covexp(rho_acc_P[0], C_P, H_P, nn_j);
F77_NAME(dscal)(&nn_j, &sigma2_acc_P[0], C_P, &incOne);
//invert C and log det cov
F77_NAME(dpotrf)(upper, &n_j, C_P, &n_j, &info); if(info != 0)
{
error("c++ error: Cholesky failed in sp.lm (N1)\n");
}
// adesso C contiene la parte superiore dellafattorizzazione di
cholesky
logdet = 0.0;
for(i = 0; i < n_j; i++) logdet += 2*log(C_P[i*n_j+i]);
F77_NAME(dpotri)(upper, &n_j, C_P, &n_j, &info);
if(info != 0)
{
error("c++ error: Cholesky inverse failed in sp.lm (N2)\n");
}
for(j=0;j<n_j-1;j++)
{
for(h=j+1;h<n_j;h++)
{
Sum_Sigma_inv += C_P[j*n_j+h] ;
}
}
Sum_Sigma_inv = Sum_Sigma_inv*2;
for(j=0;j<n_j;j++){Sum_Sigma_inv += C_P[j*n_j+j];}
for(j=0;j<n_j;j++)
{
Sigma_inv_one[j] = 0;
for(h=0;h<n_j;h++)
{
Sigma_inv_one[j] += C_P[j*n_j+h];
}
}
// /*****************************************
// // MCMC
// *****************************************/
for(MCMC_iter=0;MCMC_iter<iter_1;MCMC_iter++)
{
}
for(MCMC_iter=0;MCMC_iter<iter_2;MCMC_iter++)
{
for(MCMC_iter2=0;MCMC_iter2<thin;MCMC_iter2++)
{
indice_adapt += 1;
/************ SIMULAZIONE MU **********/
F77_NAME(dcopy)(&incJ, mu_acc_P, &incOne, mu_app_P, &incOne);
// mu_1
app1_2 =
1/((1+pow(phi_acc_P[0],2))/zeta2_acc_P[0]+Sum_Sigma_inv);
app1_1 = phi_acc_P[0]*mu_app_P[1]/zeta2_acc_P[0];
for(t=0; t<n_j;t++)
{
app1_1 +=
Sigma_inv_one[t]*(x_P[t]+2*M_PI*k_acc_P[t]-alpha_acc_P[0]);
}
app1_1 = app1_1*app1_2;
mu_acc_P[0] = rnorm(app1_1,sqrt(app1_2));
for(j=1;j<J-1; j++)
{
// mu_j
app1_1 =
(phi_acc_P[0]*mu_app_P[j-1]+phi_acc_P[0]*mu_app_P[j+1])/zeta2_acc_P[0];
for(t=0; t<n_j;t++)
{
app1_1 +=
Sigma_inv_one[t]*(x_P[t+n_j*j]+2*M_PI*k_acc_P[t+n_j*j]-alpha_acc_P[0]);
}
app1_1 = app1_1*app1_2;
mu_acc_P[j] = rnorm(app1_1,sqrt(app1_2));
}
// // mu_J
app1_2 = 1/(1/zeta2_acc_P[0]+Sum_Sigma_inv);
app1_1 = phi_acc_P[0]*mu_app_P[J-2]/zeta2_acc_P[0];
for(t=0; t<n_j;t++)
{
app1_1 +=
Sigma_inv_one[t]*(x_P[t+n_j*(J-1)]+2*M_PI*k_acc_P[t+n_j*(J-1)]-alpha_acc_P[0]);
}
app1_1 = app1_1*app1_2;
mu_acc_P[J-1] = rnorm(app1_1,sqrt(app1_2));
/************ SIMULAZIONE alpha **********/
app1_3 = 1/(J*Sum_Sigma_inv+sigma2_alpha);
app1_4 = M_alpha;
for(t=0;t<n_j;t++)
{
for(j=0;j<J;j++)
{
app1_4 +=
Sigma_inv_one[t]*(x_P[t+n_j*j]+2*M_PI*k_acc_P[t+n_j*j]-mu_app_P[j]);
}
}
app1_4 = app1_4*app1_3;
alpha_acc_P[0] = rnorm(app1_3,sqrt(app1_4) );
}
/************ SALVO LE VARIABILI **********/
F77_NAME(dcopy)(&incJ, mu_acc_P, &incOne, &mu_out_P[MCMC_iter*J],
&incOne);
F77_NAME(dcopy)(&incOne, alpha_acc_P, &incOne,
&alpha_out_P[MCMC_iter], &incOne);
}
//make return object
SEXP result,resultNames;
// oggetti della lista
int nResultListObjs = 2;
PROTECT(result = allocVector(VECSXP, nResultListObjs)); nProtect++;
PROTECT(resultNames = allocVector(VECSXP, nResultListObjs));
nProtect++;
//samples
SET_VECTOR_ELT(result, 0, mu_out_r);
SET_VECTOR_ELT(resultNames, 0, mkChar("mu"));
SET_VECTOR_ELT(result, 1, alpha_out_r);
SET_VECTOR_ELT(resultNames, 1, mkChar("alpha"));
namesgets(result, resultNames);
//unprotect
UNPROTECT(nProtect);
return(result);
//return(R_NilValue);
}
}
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
