# 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'.
*********************************************************/

#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);

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