From Martyn.Byng at nag.co.uk  Tue Aug 27 15:14:12 2013
From: Martyn.Byng at nag.co.uk (Martyn Byng)
Date: Tue, 27 Aug 2013 14:14:12 +0100
Subject: [Rd] Error in simulation. NAN
References: <521C9837.9000701@yahoo.it>
Message-ID: <49E76DF37649DC48A4CE882BC8CE51C9031AFB91@nagmail2.nag.co.uk>

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}}



[Rd] Error in simulation. NAN
Martyn Byng