[Rd] Error in simulation. NAN
Martyn Byng
Martyn.Byng at nag.co.uk
Tue Aug 27 15:14:12 CEST 2013
Hi,
You are going to have to track down exactly where the NaN is being
produced. These usually occur due to performing invalid numerical
operations, like trying to take the square-root of a negative number,
trying to calculate 0/0 etc.
Try adding a series of if / then blocks prior to performing any
operation that could result in a NaN, so something like:
if (app1_4 <= 0.0) {
error("app1_4 is -ve");
} else {
alpha_acc_P[0] = rnorm(app1_3,sqrt(app1_4) );
}
this should then give you a better idea of what is going wrong.
Martyn
-----Original Message-----
From: r-devel-bounces at r-project.org
[mailto:r-devel-bounces at r-project.org] On Behalf Of
gianluca.mastrantonio at yahoo.it
Sent: 27 August 2013 13:15
To: r-devel at r-project.org
Subject: [Rd] Error in simulation. NAN
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 at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
________________________________________________________________________
This e-mail has been scanned for all viruses by Star.\ _...{{dropped:12}}
More information about the R-devel
mailing list