betasimu.c 1/1

/*************************************************************
 * Simulates noisy data from a beta response function with
 * various alternatives of error distribution. Writes a GLIM
 * file. Parametrization is left as a responsibility of the
 * user...
 *
 * J.Oksanen, June 1997.
 **************************************************************/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <string.h>

/* Returns the expected value from beta response function */

float beta(float k, float a, float b, float alpha, float gamma,
	   float x)
{
  float t2, t3;

  /* Return zero if x is not in (a,b) */
  if (x <= a || x >= b)
    return 0;

  /* Otherwise evaluate the beta-function at x */
  t2 = pow(x-a,alpha);
  t3 = pow(b-x,gamma);
  return k*t2*t3;
}

/* The program asks the maximum height of the response function (to
   the benefit of the user). Function ksol returns the value of k.
 */

float ksol(float a, float b, float alpha, float gamma, float height)
{
  float t1, t4, t6, t11;
  t1 = b-a;
  t4 = t1/(alpha+gamma);
  t6 = pow(alpha*t4, alpha);
  t11= pow(gamma*t4, gamma);
  return height/t6/t11;
}

/* For Beta-Binomial error model: Estimates the parameters a,b of
   beta distribution from expected proportion (pi), binomial
   denominator (m), and shape parameter (tau2). Solution (hopefully
   correct) of Exercise 4.17 of McCullagh & Nelder 1989, helped by
   Moore, Appl Stat 36, 8-14; 1987.
 */

void betapara(float pi, float m, float tau2, float *a, float *b)
{
  float t1,t2,t3,t4;
  t1 = tau2*m;
  t2 = t1-m-tau2+1;
  t3 = 1/(1+t1-tau2);
  t4 = t2*t3;
  *a = -t4*pi;
  *b = t4*(pi-1);
}

/* Ranlib.c (Netlib depository) used for random number generation.
   To re-compile the program, you must either get these routines or
   replace them with your favourite rng-routines.
*/

extern void setall (long, long);          /* initialize rng */
extern void phrtsd (char*, long*, long*); /* find seeds from a string */
extern float ranf ();                     /* uniform [0,1) */
extern long ignbin (long, float);         /* Binomial */
extern long ignpoi (float);               /* Poisson */
extern float genbet (float, float);       /* beta distribution */
extern float gengam (float, float);       /* gamma */


/*********************************************************
 * main() starts.
 * The only commandline parameter is a phrase (i.e., some text) to
 * seed the rng. If not given, rng will be initialized from the clock.
 * Start from the DOS prompt with, e.g., `betapois ThisIsSeed'.
 * Asks for other parameters.
 *********************************************************/

#define BERNOULLI 0
#define BINOM 1
#define POISSON 2
#define BETABIN 3
#define NEGBIN 4

main(int argc, char *argv[])
{
  float x, k, a, b, alpha, gamma, hi, mu, enda, endb, span, temp,
    shape, aa, bb, muorig;
  int laji=0, nsimu, i, resp, error;
  long seed1,seed2,bd;
  char phrase[81]="";
  char *distr[]={"Bernoulli","Binomial, total bd = ","Poisson",
		 "Beta-Binomial, total bd = ","NegBin"};
  FILE *out;

  out = fopen("betasimu.dat","w");

  /* Initialize RNG */

  if (argc>1)
    strcpy(phrase,argv[1]);
  else
    sprintf(phrase,"%x",(long) time(NULL));
  phrtsd(phrase, &seed1, &seed2);
  setall(seed1,seed2);

  fprintf(out,"\n$COM Noisy responses simulated from beta function");
  fprintf(out,"\n$COM phrase used in seeding random numbers: %s",phrase);

  /* Ask for the error distribution */

  printf ("\nSelect the error distribution from the following alternatives:");
  printf ("\n 0 - Bernoulli (Binomial total 1)");
  printf ("\n 1 - Binomial, total asked separately");
  printf ("\n 2 - Poisson");
  printf ("\n 3 - Beta Binomial, shape parameter asked for each species");
  printf ("\n 4 - NegBin, shape parameter asked for each species");
  printf ("\n ");
  scanf ("%d", &error);
  if (error == BINOM || error == BETABIN) {
    printf("\nGive the binomial denominator, bd (integer)  ");
    scanf ("%ld",&bd);
  }

  /* Print info on error distribution used */

  fprintf(out,"\n$COM error distribution %s ",distr[error]);
  if (error == BINOM || error == BETABIN)
    fprintf(out," %d",bd);
  
  /* Define the gradient */
  
  printf ("Gradient endpoints (2 needed)      ");
  scanf ("%f %f", &enda, &endb);
  if (endb < enda) {
    temp = enda;
    enda = endb;
    endb = temp;
  }
  span = endb-enda;
  fprintf(out,"\n$COM Gradient %.4f .. %.4f (span %.4f)",enda,endb,span);
  printf ("Number of simulated observations   ");
  scanf ("%d", &nsimu);
  fprintf(out, "\n$units %d", nsimu);
  if (error == BINOM || error == BETABIN)
    fprintf(out,"\n$cal bd = %d ! Binomial Denominator", bd);
  if (error == BERNOULLI)
    fprintf(out,"\n$cal bd = 1 ! Binomial Denominator");

  printf ("\nThe parameters of beta response function will be asked.");
  printf ("\nTo end: Give expected max abundance of zero");

  /* Asking beta parameters for each species.
     To end: Give expected maximum abundance 0.
     Asks expected maximum abundance instead of parameter k, since
     I think this is more convenient to the user. The value of k is
     solved by ksol, and will be reported in the output.
   */
  for (;;) {
    printf(  "\n\nGive expected max abundance (0 to end) ");
    if (error == BINOM || error == BETABIN || error == BERNOULLI)
      printf  ("\nValue should be probability <1         ");
    scanf ("%f", &hi);
    if (hi <= 0) break; /* Break from the for(;;) loop */
    printf("\nGive range endpoints a, b              ");
    scanf("%f %f", &a, &b);
    if (b < a) {
      temp = a;
      a = b;
      b = temp;
    }
    printf(  "\nGive shape parameters alpha, gamma     ");
    scanf ("%f %f", &alpha, &gamma);
    k = ksol(a, b, alpha, gamma, hi);

    ++laji;
    fprintf (out,"\n$com y%d: a=%.2f, b=%.2f, alpha=%.2f, gamma=%.2f, hi=%.2f (k=%.2f)",
	     laji,a,b,alpha,gamma,hi,k);
    if (error == BETABIN || error == NEGBIN) {
      printf(  "\nGive the overdispersion parameter for species %d ",laji);
      if (error == BETABIN)
	printf("\ntau2 in Var(y) = m*pi*(1-pi)*(1+(m-1)*tau2)      ");
      else
	printf("\nsigma2 in Var(y) = mu + sigma2*mu^2              ");
      scanf ("%f",&shape);
      fprintf(out,"\n$com Overdispersion parameter %.4f", shape);
    }

    /* Proper simulation of observations starts here. Ugly `if (mu > 0.0)'
       are needed since some (!) values of zero seem to cause FP errors in
       ranlib.c. Please note1 the fall-through in BETABIN and NEGBIN.
     */
    
    fprintf(out,"\n$data x mu%d y%d\n$read \n", laji,laji);
    for (i=0; i<nsimu; i++) {
      x = (float)i/(float)(nsimu-1)*span+enda;
      muorig = mu = beta(k,a,b,alpha,gamma,x);
      if (error == BINOM || error == BETABIN)
	muorig *= bd;
      switch (error) {
	case BERNOULLI:
	  if (mu >= ranf()) resp=1;
	  else resp=0;
	  break;
	case BETABIN:
	  if (mu > 0.0) {
	    betapara(mu,bd,shape,&aa,&bb);
	    mu = genbet(aa,bb);
	  }
	case BINOM:
	  if (mu > 0.0) resp = (int) ignbin(bd,mu);
	  else resp = 0;
	  break;
	case NEGBIN:
	  if (mu > 0.0) mu = gengam(1/(shape*mu),1/shape);
	case POISSON:
	  if (mu > 0.0) resp = (int) ignpoi(mu);
	  else resp = 0;
	  break;
      }
      fprintf(out,"%7.4f %7.4f %3d\n", x, muorig, resp);
    }
  }
  fprintf (out,"\n$return");
  printf ("\nWrote GLIM file BETASIMU.DAT");
}


Generated by GNU enscript 1.5.0.