#include
#include
#include
#include
/*
Program "fprior.c" (12/06/97)
Updated for ANSI C; 18/07/00 - to compile type cc fprior.c -lm
Karen L. Ayres
University of Reading
e-mail: k.l.ayres@reading.ac.uk
Generates a sample from the joint prior for (f,p) for "nall"
alleles (where "nall" is specified in the main part of the program)
as used in the program "hwmet.c" (which itself produces a sample
from the marginal posterior density of the inbreeding coefficient f
by means of a Metropolis-Hastings algorithm which allows for
negative values of f, provided that the condition p_ii>=0 is
satisfied, this condition causing the dependence of f and p - this
algorithm is detailed in Ayres & Balding, Heredity, 1998).
Outputs "f" parameter values to file "fprior.out", and "p" vectors
to "fpriorp.out".
The {p_i} have a multivariate uniform prior, and f has a Beta prior
conditional on the vector p, on the interval (A,1),
where A=-pmin/(1-pmin) and pmin is the smallest value in p (assuming
all p[i] are nonzero), with parameters
alpha=1-A.(beta-1)
beta=g_t(pmin)
and where g_t(*) is a polynomial function of pmin, whose form is
determined by the value "t", which represents where the right-hand
tail of the prior density curve is to become approximately 0. This
program currently caters for three values of "t" - 0.1, 0.2, and 0.4,
which are likely to be useful for many human populations since they
give most support to values of f close to zero. The conditional priors
all have mode at zero, which may be reasonable a priori. The values
of "t" in the program are coded as follows:
for t=0.1 : set btail = 1
for t=0.2 : set btail = 2
for t=0.4 : set btail = 4
for a uniform (i.e. non-informative) prior : set btail = 0
The parameter "btail" is defined at the beginning of the program,
and should be set to one of the above integers accordingly.
The prior density curves can be produced within a suitable statistical
computer package by applying density estimation to the output from this
program.
*/
#define maxall 15 /* maximum number of alleles allowed at a locus */
#define btail 0 /* determines right-hand tail of conditional beta prior of f */
typedef double pvec[maxall];
int i, nall;
long ix, iy, iz, j;
double pmin, f, blb, bpa, bpb, gsum;
pvec p;
FILE *inp, *out1, *out2;
double WichHill()
{
/* Wichmann & Hill random number generator (Algorithm AS183 in Applied
Statistics 31 (1982))
*/
ix=(171*ix)%30269;
iy=(172*iy)%30307;
iz=(170*iz)%30323;
return fmod((ix/30269.0)+(iy/30307.0)+(iz/30323.0),1.0);
}
double rstgam(double alpha1)
{
/* Cheng & Feast algorithm (algorithm 3.20 in Ripley, 1988) for
generating a standard Gamma random variable with shape
parameter alpha>1
*/
int acpt;
double c1, c2, c3, c4, c5, u1, u2, w;
c1=alpha1-1.0;
c2=(alpha1-(1.0/(6.0*alpha1)))/c1;
c3=2.0/c1;
c4=c3+2.0;
c5=1.0/sqrt(alpha1);
do{
do{
u1=WichHill();
u2=WichHill();
if (alpha1>2.5){
u1=u2+c5*(1.0-(1.86*u1));
}
}
while (u1>=1.0 || u1<=0);
w=(c2*u2)/u1;
if ((c3*u1+w+(1.0/w))<=c4){
acpt=1;
}
else if ((c3*log(u1)-log(w)+w)>=1.0){ /* repeat outer "do" loop */
acpt=0;
}
else{
acpt=1;
}
}
while (acpt!=1);
return c1*w;
}
double rgam(double pp)
{
/* returns a standard gamma random variable, with shape parameter "pp" */
double rgm;
if (pp>1.0){
rgm=rstgam(pp);
}
else if (pp<1.0){ /* Ex. 3.22 in Ripley(1988) for alpha<1 */
rgm=rstgam(pp+1.0)*pow(WichHill(),1.0/pp);
}
else { /* Exponential */
rgm= (-log(WichHill()));
}
return rgm;
}
main()
{
srand((unsigned) time(NULL));
out1=fopen("fprior.out","w");
out2=fopen("fpriorp.out","w");
/* initialize seeds for random number generator */
ix=(long)(rand()%30000)+1;
iy=(long)(rand()%30000)+1;
iz=(long)(rand()%30000)+1;
/* set number of alleles */
nall=2;
/* generate a random sample from the joint prior of (f,p) */
for (j=0;j<50000L;++j){
gsum=0;
for (i=0;i