https://github.com/cran/aster
Tip revision: 63c8e8a453ea587001e2438a8ce51cf0e1e1675c authored by Charles J. Geyer on 23 March 2009, 00:00:00 UTC
version 0.7-7
version 0.7-7
Tip revision: 63c8e8a
astfam.c
#include <math.h>
#include <stddef.h>
#include <string.h>
#include "aster.h"
#include "raster.h"
/* Bernoulli */
#ifndef __GNUC__
static int bernoulli_parval(double xpred, double hyper1, double hyper2)
#else
static int bernoulli_parval(double xpred,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* must be nonnegative integer */
int foo = 1;
foo = foo && (xpred == ceil(xpred));
foo = foo && (xpred >= 0.0);
return foo;
}
static int bernoulli_validate(double x, double xpred, double hyper1,
double hyper2)
{
/* xpred must be valid */
/* x must be integer between 0 and xpred */
int foo = 1;
foo = foo && bernoulli_parval(xpred, hyper1, hyper2);
foo = foo && (x == ceil(x));
foo = foo && (x >= 0.0);
foo = foo && (x <= xpred);
return foo;
}
#ifndef __GNUC__
static int bernoulli_hypval(double hyper1, double hyper2)
#else
static int bernoulli_hypval(double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
return 1;
}
#ifndef __GNUC__
static double bernoulli(int deriv, double theta, double hyper1, double hyper2)
#else
static double bernoulli(int deriv, double theta,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
double foo, bar;
switch (deriv) {
case 0:
if (theta <= 0)
return my_log1p(exp(theta));
else
return theta + my_log1p(exp(- theta));
case 1:
return 1.0 / (1.0 + exp(- theta));
case 2:
theta = fabs(theta);
foo = exp(- theta);
bar = 1.0 + foo;
return (foo / bar) / bar;
default:
die("deriv %d not valid", deriv);
}
}
#ifndef __GNUC__
static double bernoulli_simulate(double xpred, double theta, double hyper1,
double hyper2)
#else
static double bernoulli_simulate(double xpred, double theta,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
if (xpred == 0.0) {
return 0.0;
} else {
double p = 1.0 / (1.0 + exp(- theta));
return my_rbinom(xpred, p);
}
}
/* Poisson */
#ifndef __GNUC__
static int poisson_parval(double xpred, double hyper1, double hyper2)
#else
static int poisson_parval(double xpred,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* must be nonnegative real */
return (xpred >= 0.0);
}
static int poisson_validate(double x, double xpred, double hyper1,
double hyper2)
{
/* xpred must be valid */
/* x must be nonnegative integer */
/* xpred == 0 implies x == 0 */
int foo = 1;
foo = foo && poisson_parval(xpred, hyper1, hyper2);
foo = foo && (x == ceil(x));
foo = foo && (x >= 0.0);
foo = foo && (xpred > 0.0 || x == 0.0);
return foo;
}
#ifndef __GNUC__
static int poisson_hypval(double hyper1, double hyper2)
#else
static int poisson_hypval(double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
return 1;
}
#ifndef __GNUC__
static double poisson(int deriv, double theta, double hyper1, double hyper2)
#else
static double poisson(int deriv, double theta,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
switch (deriv) {
case 0:
case 1:
case 2:
return exp(theta);
default:
die("deriv %d not valid", deriv);
}
}
#ifndef __GNUC__
static double poisson_simulate(double xpred, double theta, double hyper1,
double hyper2)
#else
static double poisson_simulate(double xpred, double theta,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
double mu = xpred * exp(theta);
if (mu == 0.0)
return 0.0;
return my_rpois(mu);
}
#ifdef ASTER_OLD_STUFF
/* zero-truncated Poisson (now obsolete -- replace by k-truncated) */
static int non_zero_poisson_validate(double x, double xpred,
double hyper1, double hyper2)
{
int foo = 1;
foo = foo && (xpred == my_round(xpred));
foo = foo && (x == my_round(x));
foo = foo && (xpred >= 0.0);
foo = foo && (x >= 0.0);
foo = foo && (xpred > 0.0 || x == 0.0);
foo = foo && (xpred == 0.0 || x > 0.0);
return foo;
}
static double non_zero_poisson(int deriv, double theta, double hyper1,
double hyper2)
{
double mu = exp(theta);
double tau;
if (1.0 - mu == 1.0)
tau = 1.0;
else
tau = mu / (- my_expm1(- mu));
switch (deriv) {
case 0:
return mu + my_log1p(- exp(- mu));
case 1:
return tau;
case 2:
return tau * (1.0 - tau * exp(- mu));
default:
die("deriv %d not valid", deriv);
}
}
static double non_zero_poisson_simulate(double xpred, double theta,
double hyper1, double hyper2)
{
double mu = exp(theta);
double result = 0.0;
int i;
for (i = 0; i < xpred; ++i)
result += my_rnzp(mu);
return result;
}
/* two-truncated Poisson (now obsolete -- replace by k-truncated) */
static int two_trunc_poisson_validate(double x, double xpred,
double hyper1, double hyper2)
{
int foo = 1;
foo = foo && (xpred == my_round(xpred));
foo = foo && (x == my_round(x));
foo = foo && (xpred >= 0.0);
foo = foo && (x >= 0.0);
foo = foo && (xpred > 0.0 || x == 0.0);
foo = foo && (xpred == 0.0 || x > 2.0);
return foo;
}
static double two_trunc_poisson(int deriv, double theta, double hyper1,
double hyper2)
{
int k = 2;
double mu = exp(theta);
double foo, bar, baz, qux, alpha;
switch (deriv) {
case 0:
return mu + my_ppois(k, mu, 0, 1);
case 1:
foo = my_ppois(k + 1, mu, 0, 0);
if (foo == 0.0) {
return mu + k + 1;
} else {
bar = my_dpois(k + 1, mu, 0);
return mu + (k + 1) / (1.0 + foo / bar);
}
case 2:
foo = my_ppois(k + 1, mu, 0, 0);
if (foo == 0.0) {
qux = 0.0;
} else {
bar = my_dpois(k + 1, mu, 0);
qux = foo / bar;
}
alpha = (k + 1) / (1.0 + qux);
if (qux < 1.0) {
baz = qux / (1.0 + qux);
} else {
baz = 1.0 / (1.0 / qux + 1.0);
}
return mu * (1.0 - alpha * (1.0 - (k + 1) * baz / mu));
default:
die("deriv %d not valid", deriv);
}
}
static double two_trunc_poisson_simulate(double xpred, double theta,
double hyper1, double hyper2)
{
double mu = exp(theta);
double result = 0.0;
int i;
for (i = 0; i < xpred; ++i)
result += my_rktp(2, mu);
return result;
}
#endif /* ASTER_OLD_STUFF */
/* k-truncated Poisson */
#ifndef __GNUC__
static int trunc_poisson_parval(double xpred, double hyper1, double hyper2)
#else
static int trunc_poisson_parval(double xpred,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* must be nonnegative integer */
int foo = 1;
foo = foo && (xpred == ceil(xpred));
foo = foo && (xpred >= 0.0);
return foo;
}
static int trunc_poisson_validate(double x, double xpred,
double hyper1, double hyper2)
{
double k = hyper1;
/* xpred must be valid */
/* x must be integer */
/* xpred == 0 implies x == 0 */
/* xpred > 0 implies x > k */
int foo = 1;
foo = foo && trunc_poisson_parval(xpred, hyper1, hyper2);
foo = foo && (x == ceil(x));
foo = foo && (xpred > 0.0 || x == 0.0);
foo = foo && (xpred == 0.0 || x > k);
return foo;
}
#ifndef __GNUC__
static int trunc_poisson_hypval(double hyper1, double hyper2)
#else
static int trunc_poisson_hypval(double hyper1,
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* hyper1 (truncation) must be nonnegative integer */
int foo = 1;
foo = foo && (hyper1 == ceil(hyper1));
foo = foo && (hyper1 >= 0.0);
return foo;
}
#ifndef __GNUC__
static double trunc_poisson(int deriv, double theta,
double hyper1, double hyper2)
#else
static double trunc_poisson(int deriv, double theta, double hyper1,
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
int k = hyper1;
double mu = exp(theta);
double foo, bar, baz, qux, alpha;
switch (deriv) {
case 0:
return mu + my_ppois(k, mu, 0, 1);
case 1:
foo = my_ppois(k + 1, mu, 0, 0);
if (foo == 0.0) {
return mu + k + 1;
} else {
bar = my_dpois(k + 1, mu, 0);
return mu + (k + 1) / (1.0 + foo / bar);
}
case 2:
foo = my_ppois(k + 1, mu, 0, 0);
if (foo == 0.0) {
qux = 0.0;
} else {
bar = my_dpois(k + 1, mu, 0);
qux = foo / bar;
}
alpha = (k + 1) / (1.0 + qux);
if (qux < 1.0) {
baz = qux / (1.0 + qux);
} else {
baz = 1.0 / (1.0 / qux + 1.0);
}
return mu * (1.0 - alpha * (1.0 - (k + 1) * baz / mu));
default:
die("deriv %d not valid", deriv);
}
}
#ifndef __GNUC__
static double trunc_poisson_simulate(double xpred, double theta,
double hyper1, double hyper2)
#else
static double trunc_poisson_simulate(double xpred, double theta,
double hyper1, double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
int k = hyper1;
double mu = exp(theta);
double result = 0.0;
int i;
for (i = 0; i < xpred; ++i)
result += my_rktp(k, mu);
return result;
}
/* (untruncated) negative binomial */
#ifndef __GNUC__
static int neg_bin_parval(double xpred, double hyper1, double hyper2)
#else
static int neg_bin_parval(double xpred,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* must be nonnegative real */
return (xpred >= 0.0);
}
static int neg_bin_validate(double x, double xpred,
double hyper1, double hyper2)
{
/* xpred must be valid */
/* x must be integer */
/* xpred == 0 implies x == 0 */
int foo = 1;
foo = foo && neg_bin_parval(xpred, hyper1, hyper2);
foo = foo && (x == ceil(x));
foo = foo && (xpred > 0.0 || x == 0.0);
return foo;
}
#ifndef __GNUC__
static int neg_bin_hypval(double hyper1, double hyper2)
#else
static int neg_bin_hypval(double hyper1,
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* hyper1 (size) must be nonnegative real */
return (hyper1 >= 0.0);
}
#ifndef __GNUC__
static double neg_bin(int deriv, double theta,
double hyper1, double hyper2)
#else
static double neg_bin(int deriv, double theta, double hyper1,
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
double alpha = hyper1;
double p, q, foo;
if (theta >= 0.0) {
switch (deriv) {
case 0:
return my_posinf();
case 1:
case 2:
return my_nan();
default:
die("deriv %d not valid", deriv);
}
}
switch (deriv) {
case 0:
return alpha * my_log1p(1.0 / my_expm1(- theta));
case 1:
case 2:
q = exp(theta);
p = (- my_expm1(theta));
foo = alpha * q / p;
if (deriv == 1)
return foo;
else
return foo / p;
default:
die("deriv %d not valid", deriv);
}
}
#ifndef __GNUC__
static double neg_bin_simulate(double xpred, double theta,
double hyper1, double hyper2)
#else
static double neg_bin_simulate(double xpred, double theta,
double hyper1, double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
double alpha = hyper1;
double p = - my_expm1(theta);
return my_rnbinom(xpred * alpha, p);
}
/* truncated negative binomial */
#ifndef __GNUC__
static int trunc_neg_bin_parval(double xpred, double hyper1, double hyper2)
#else
static int trunc_neg_bin_parval(double xpred,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* must be nonnegative integer */
int foo = 1;
foo = foo && (xpred == ceil(xpred));
foo = foo && (xpred >= 0.0);
return foo;
}
static int trunc_neg_bin_validate(double x, double xpred,
double hyper1, double hyper2)
{
double k = hyper2;
/* xpred must be valid */
/* x must be integer */
/* xpred == 0 implies x == 0 */
/* xpred > 0 implies x > k */
int foo = 1;
foo = foo && trunc_neg_bin_parval(xpred, hyper1, hyper2);
foo = foo && (x == ceil(x));
foo = foo && (xpred > 0.0 || x == 0.0);
foo = foo && (xpred == 0.0 || x > k);
return foo;
}
static int trunc_neg_bin_hypval(double hyper1, double hyper2)
{
/* hyper1 (size) must be nonnegative real */
/* hyper2 (truncation) must be nonnegative integer */
int foo = 1;
foo = foo && (hyper1 >= 0.0);
foo = foo && (hyper2 == ceil(hyper2));
foo = foo && (hyper2 >= 0.0);
return foo;
}
static double trunc_neg_bin(int deriv, double theta,
double hyper1, double hyper2)
{
double alpha = hyper1;
double k = hyper2;
double p, q, foo, mu, numer, beta, beep;
if (theta >= 0.0) {
switch (deriv) {
case 0:
return my_posinf();
case 1:
case 2:
return my_nan();
default:
die("deriv %d not valid", deriv);
}
}
switch (deriv) {
case 0:
/* psi(theta) = - log(1 - exp(theta)) + log Pr(Y > k) */
p = - my_expm1(theta);
return alpha * my_log1p(1.0 / my_expm1(- theta)) +
my_pnbinom(k, alpha, p, 0, 1);
case 1:
case 2:
q = exp(theta);
p = - my_expm1(theta);
mu = alpha * q / p;
numer = my_pnbinom(k + 1, alpha, p, 0, 0);
if (numer == 0.0)
beta = 0.0;
else
beta = numer / my_dnbinom(k + 1, alpha, p, 0);
if (deriv == 1)
return mu + (k + 1) / (1.0 + beta) / p;
if (beta < 1.0)
beep = beta / (1.0 + beta);
else
beep = 1.0 / (1.0 / beta + 1.0);
foo = k + 1.0 + alpha;
foo = (- q) + foo * q / (1.0 + beta) + (alpha - p * foo) * beep;
foo *= (k + 1.0) / (1.0 + beta) / p;
foo = mu - foo;
foo /= p;
return foo;
default:
die("deriv %d not valid", deriv);
}
}
static double trunc_neg_bin_simulate(double xpred, double theta,
double hyper1, double hyper2)
{
double alpha = hyper1;
double k = hyper2;
double p = (- my_expm1(theta));
double q = exp(theta);
double mu = alpha * q / p;
double result = 0.0;
int i;
for (i = 0; i < xpred; ++i)
result += my_rktnb(alpha, k, mu);
return result;
}
/* normal location (known scale) */
#ifndef __GNUC__
static int norm_loc_parval(double xpred, double hyper1, double hyper2)
#else
static int norm_loc_parval(double xpred,
double hyper1 __attribute__ ((unused)),
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* must be nonnegative real */
return (xpred >= 0.0);
}
static int norm_loc_validate(double x, double xpred,
double hyper1, double hyper2)
{
/* xpred must be valid */
/* x must be real */
/* xpred == 0 implies x == 0 */
int foo = 1;
foo = foo && norm_loc_parval(xpred, hyper1, hyper2);
foo = foo && (xpred > 0.0 || x == 0.0);
return foo;
}
#ifndef __GNUC__
static int norm_loc_hypval(double hyper1, double hyper2)
#else
static int norm_loc_hypval(double hyper1,
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
/* hyper1 (size) must be positive real */
return (hyper1 > 0.0);
}
#ifndef __GNUC__
static double norm_loc(int deriv, double theta,
double hyper1, double hyper2)
#else
static double norm_loc(int deriv, double theta, double hyper1,
double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
double sigmasq = hyper1 * hyper1;
switch (deriv) {
case 0:
return sigmasq * theta * theta / 2.0;
case 1:
return sigmasq * theta;
case 2:
return sigmasq;
default:
die("deriv %d not valid", deriv);
}
}
#ifndef __GNUC__
static double norm_loc_simulate(double xpred, double theta,
double hyper1, double hyper2)
#else
static double norm_loc_simulate(double xpred, double theta,
double hyper1, double hyper2 __attribute__ ((unused)))
#endif /* __GNUC__ */
{
double mu = theta * hyper1 * hyper1 * xpred;
double sigma = hyper1 * sqrt(xpred);
return my_rnorm(mu, sigma);
}
typedef double (*famfun_ptr)(int deriv, double theta, double hyper1,
double hyper2);
typedef int (*famval_ptr)(double x, double xpred, double hyper1,
double hyper2);
typedef int (*fampval_ptr)(double xpred, double hyper1, double hyper2);
typedef int (*famhval_ptr)(double hyper1, double hyper2);
typedef double (*famsim_ptr)(double xpred, double theta, double hyper1,
double hyper2);
struct superfamtab {
char * const name;
famfun_ptr const psi;
famval_ptr const validate;
fampval_ptr const validparent;
famhval_ptr const validhyper;
famsim_ptr const simulate;
int const mincard;
int const maxcard;
int const nhyper;
char * const name_hyper1;
char * const name_hyper2;
double const origin;
};
static struct superfamtab const mysuperfamtab[] =
{
{"bernoulli", bernoulli, bernoulli_validate, bernoulli_parval,
bernoulli_hypval, bernoulli_simulate, 1, 1, 0, "", "", 0.0},
{"poisson", poisson, poisson_validate, poisson_parval,
poisson_hypval, poisson_simulate, 1, 1, 0, "", "", 0.0},
{"truncated.poisson", trunc_poisson, trunc_poisson_validate,
trunc_poisson_parval, trunc_poisson_hypval, trunc_poisson_simulate,
1, 1, 1, "truncation", "", 0.0},
{"negative.binomial", neg_bin, neg_bin_validate,
neg_bin_parval, neg_bin_hypval, neg_bin_simulate,
1, 1, 1, "size", "", -1.0},
{"truncated.negative.binomial", trunc_neg_bin, trunc_neg_bin_validate,
trunc_neg_bin_parval, trunc_neg_bin_hypval, trunc_neg_bin_simulate,
1, 1, 2, "size", "truncation", -1.0},
{"normal.location", norm_loc, norm_loc_validate,
norm_loc_parval, norm_loc_hypval, norm_loc_simulate,
1, 1, 1, "sd", "", 0.0},
{NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0.0},
};
struct famtab {
char *name;
famfun_ptr psi;
famval_ptr validate;
fampval_ptr validparent;
famsim_ptr simulate;
int mincard;
int maxcard;
int nhyper;
double hyper1;
double hyper2;
char *name_hyper1;
char *name_hyper2;
double origin;
};
static struct famtab myfamtab[MAXNFAM];
static int nfam = 0;
void aster_clear_families(void)
{
nfam = 0;
}
void aster_add_family(char **name, double *hyper, int *nhyper)
{
int ifam;
double hyper1 = 0.0;
double hyper2 = 0.0;
if (nfam == MAXNFAM)
die("number of families exceeds family table size");
for (ifam = 0; ; ++ifam) {
if (mysuperfamtab[ifam].name == NULL)
die("family \"%s\" not found", name[0]);
if (strcmp(mysuperfamtab[ifam].name, name[0]) == 0)
break;
}
if (mysuperfamtab[ifam].nhyper != nhyper[0])
die("family \"%s\" has %d hyperparameters, %d specified",
name[0], mysuperfamtab[ifam].nhyper, nhyper[0]);
if (nhyper[0] >= 1)
hyper1 = hyper[0];
if (nhyper[0] >= 2)
hyper2 = hyper[1];
if (! mysuperfamtab[ifam].validhyper(hyper1, hyper2))
die("family \"%s\" specified with invalid hyperparameters", name[0]);
myfamtab[nfam].name = mysuperfamtab[ifam].name;
myfamtab[nfam].psi = mysuperfamtab[ifam].psi;
myfamtab[nfam].validate = mysuperfamtab[ifam].validate;
myfamtab[nfam].validparent = mysuperfamtab[ifam].validparent;
myfamtab[nfam].simulate = mysuperfamtab[ifam].simulate;
myfamtab[nfam].mincard = mysuperfamtab[ifam].mincard;
myfamtab[nfam].maxcard = mysuperfamtab[ifam].maxcard;
myfamtab[nfam].nhyper = mysuperfamtab[ifam].nhyper;
myfamtab[nfam].hyper1 = hyper1;
myfamtab[nfam].hyper2 = hyper2;
myfamtab[nfam].name_hyper1 = mysuperfamtab[ifam].name_hyper1;
myfamtab[nfam].name_hyper2 = mysuperfamtab[ifam].name_hyper2;
myfamtab[nfam].origin = mysuperfamtab[ifam].origin;
++nfam;
}
void aster_get_family(int *famin, char **name, double *hyper, int *nhyper,
char **hypername, double *origin)
{
int fam = famin[0];
if (fam <= 0 || fam > nfam) {
name[0] = "";
return;
}
int ifam = fam - 1; /* fam is 1-origin index */
name[0] = myfamtab[ifam].name;
nhyper[0] = myfamtab[ifam].nhyper;
origin[0] = myfamtab[ifam].origin;
if (nhyper[0] >= 1) {
hyper[0] = myfamtab[ifam].hyper1;
hypername[0] = myfamtab[ifam].name_hyper1;
}
if (nhyper[0] >= 2) {
hyper[1] = myfamtab[ifam].hyper2;
hypername[1] = myfamtab[ifam].name_hyper2;
}
return;
}
void aster_get_superfamily(int *famin, char **name, int *nhyper,
char **hypername)
{
int fam = famin[0];
if (fam <= 0) {
name[0] = "";
return;
}
int ifam = fam - 1; /* fam is 1-origin index */
for (int i = 0; i <= ifam; ++i)
if (mysuperfamtab[ifam].name == NULL) {
name[0] = "";
return;
}
name[0] = mysuperfamtab[ifam].name;
nhyper[0] = mysuperfamtab[ifam].nhyper;
if (nhyper[0] >= 1)
hypername[0] = mysuperfamtab[ifam].name_hyper1;
if (nhyper[0] >= 2)
hypername[1] = mysuperfamtab[ifam].name_hyper2;
return;
}
void aster_byname_superfamily(char **name, int *nhyper, char **hypername)
{
int ifam;
for (ifam = 0; ; ++ifam) {
if (mysuperfamtab[ifam].name == NULL)
die("family \"%s\" not found", name[0]);
if (strcmp(mysuperfamtab[ifam].name, name[0]) == 0)
break;
}
nhyper[0] = mysuperfamtab[ifam].nhyper;
if (nhyper[0] >= 1)
hypername[0] = mysuperfamtab[ifam].name_hyper1;
if (nhyper[0] >= 2)
hypername[1] = mysuperfamtab[ifam].name_hyper2;
return;
}
void aster_family(int *famin, int *derivin, double *thetain, double *value)
{
int fam = famin[0];
int deriv = derivin[0];
double theta = thetain[0];
if (fam <= 0 || fam > nfam)
die("family %d not valid", fam);
int ifam = fam - 1; /* fam is 1-origin index */
double hyper1 = myfamtab[ifam].hyper1;
double hyper2 = myfamtab[ifam].hyper2;
value[0] = myfamtab[ifam].psi(deriv, theta, hyper1, hyper2);
}
int aster_family_validate(int fam, double x, double xpred)
{
if (fam <= 0 || fam > nfam)
die("family %d not valid", fam);
int ifam = fam - 1; /* fam is 1-origin index */
double hyper1 = myfamtab[ifam].hyper1;
double hyper2 = myfamtab[ifam].hyper2;
return myfamtab[ifam].validate(x, xpred, hyper1, hyper2);
}
int aster_family_validate_parent(int fam, double xpred)
{
if (fam <= 0 || fam > nfam)
die("family %d not valid", fam);
int ifam = fam - 1; /* fam is 1-origin index */
double hyper1 = myfamtab[ifam].hyper1;
double hyper2 = myfamtab[ifam].hyper2;
return myfamtab[ifam].validparent(xpred, hyper1, hyper2);
}
int aster_family_number_validate(int fam)
{
return fam > 0 && fam <= nfam;
}
double aster_family_origin(int fam)
{
if (fam <= 0 || fam > nfam)
die("family %d not valid", fam);
int ifam = fam - 1; /* fam is 1-origin index */
return myfamtab[ifam].origin;
}
double aster_family_simulate(int fam, double xpred, double theta)
{
if (fam <= 0 || fam > nfam)
die("family %d not valid", fam);
int ifam = fam - 1; /* fam is 1-origin index */
double hyper1 = myfamtab[ifam].hyper1;
double hyper2 = myfamtab[ifam].hyper2;
return myfamtab[ifam].simulate(xpred, theta, hyper1, hyper2);
}
#ifdef ASTER_OLD_STUFF
struct funtab {
char *name;
famfun_ptr fun;
famval_ptr validate;
famsim_ptr simulate;
};
static struct funtab myfuntab[] =
{
{"bernoulli", bernoulli, bernoulli_validate, bernoulli_simulate},
{"poisson", poisson, poisson_validate, poisson_simulate},
{"non.zero.poisson", non_zero_poisson, non_zero_poisson_validate,
non_zero_poisson_simulate},
{"two.trunc.poisson", two_trunc_poisson, two_trunc_poisson_validate,
two_trunc_poisson_simulate},
{NULL, NULL, NULL, NULL},
};
void aster_family(int *famin, int *derivin, double *thetain, double *value)
{
int fam = famin[0];
int deriv = derivin[0];
double theta = thetain[0];
int i;
for (i = 0; myfuntab[i].fun != NULL; ++i)
if (i == (fam - 1)) /* fam is 1-origin index */ {
value[0] = myfuntab[i].fun(deriv, theta, 0.0, 0.0);
return;
}
die("family %d not valid", fam);
}
int aster_family_validate(int fam, double x, double xpred)
{
int i;
for (i = 0; myfuntab[i].validate != NULL; ++i)
if (i == (fam - 1)) /* fam is 1-origin index */
return myfuntab[i].validate(x, xpred, 0.0, 0.0);
die("family %d not valid", fam);
}
char *aster_family_name(int fam)
{
int i;
for (i = 0; myfuntab[i].name != NULL; ++i)
if (i == (fam - 1)) /* fam is 1-origin index */
return myfuntab[i].name;
return NULL;
}
double aster_family_simulate(int fam, double xpred, double theta)
{
int i;
for (i = 0; myfuntab[i].simulate != NULL; ++i)
if (i == (fam - 1)) /* fam is 1-origin index */
return myfuntab[i].simulate(xpred, theta, 0.0, 0.0);
die("family %d not valid", fam);
}
#endif /* ASTER_OLD_STUFF */