[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