[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] gsl_ran_beta can return 0 or 1
From: |
Toby Johnson |
Subject: |
[Bug-gsl] gsl_ran_beta can return 0 or 1 |
Date: |
Mon, 15 Nov 2004 17:28:36 +0000 |
User-agent: |
Mutt/1.5.6+20040722i |
I'm not sure whether this is a bug.
gsl_ran_beta() can return values that are 0 or 1. This can cause
problems if I then want to look up the pdf for the value just
generated by calling gsl_ran_beta_pdf, because the pdf can be infinite
at 0 or 1. This sequence of events is often used in a
Metropolis--Hastings algorithm. The following code generates this
`bug' on the 10th iteration, using the default generator and seed. I
haven't tried to fix this since I'm not sure whether others will
consider it a bug.
#include <iostream>
using namespace std;
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
int main() {
const gsl_rng_type *rng_type;
long int rng_seed;
gsl_rng *rng;
gsl_rng_env_setup();
rng_type = gsl_rng_default;
rng_seed = gsl_rng_default_seed;
rng = gsl_rng_alloc (rng_type);
gsl_rng_set(rng,rng_seed);
double x,newx,probx;
x = 0.5;
do {
newx = gsl_ran_beta(rng,5.*x,5.*(1.-x));
cout << x << " -> " << newx << endl;
probx = gsl_ran_beta_pdf(newx,5.*x,5.*(1.-x));
cout << "p = " << probx << endl;
x = newx;
} while (1);
return(0);
}
signature.asc
Description: Digital signature
- [Bug-gsl] gsl_ran_beta can return 0 or 1,
Toby Johnson <=