#include <math.h>
#include <stdio.h>
/* These versions of the routines allow for recombination between the
loci.
The function get.2locexp.recomb() takes as input the 2-locus
sperm-typing data counts (length 16) in order
-,b,a,ab,B,Bb,Ba,Bab,A,Ab,Aa,Aab,AB,ABb,ABa,ABab
and the current parameter values (length 11) in order
s (segregation probability for allele 1, assoc. with AB in parent)
p=(p1,p2) probabilities of 1 and 2 sperm
deposited with p0 assumed = 1-p1-p2, am vector of 4 amplification
probabilities for alleles A,B,a, and b; co vector of 4 contamination
probabilities for alleles A,B,a, and b. Also takes as input the
multilocus recombination probabilities: Pr(no rec., i.e. 1AB or 2ab),
Pr(rec btw disease and marker haplotype, i.e. 1ab or 2AB), Pr(1Ab or
2aB), Pr(1aB or 2Ab), where parental haplotypes are 1AB and 2ab.
The output of get.2locexp.recomb() is the expected values of the
following
quantities (length 21), conditional on the data and the current
parameter values:
number of segregations of allele 1, associated in parent with AB,
number of segregations of allele 2, associated in parent with ab,
number of times 1 sperm is deposited, number of times 2 sperm are
deposited, number of times 0 sperm are deposited,
number of times
allele A is amplified, number of times allele A is present but not
amplified, number of times allele B is amplified, number of times B is
present but not amplified, similarly for a and b, number of times
allele
A is a contaminant, number of times allele A is not a contaminant,
number of times allele B is/is not a contaminant, a is/is not a
contaminant, b is/is not a contaminant.
The intention is that this function will be run on each two-locus data
set, the function get.1locexp.recomb() will be run on each one-locus
data
set, and then the results of all these runs will be used in the next
step of the E-M algorithm to estimate the parameters (under selected
constraints) by maximum likelihood.*/
/*x is 16x4, ith row of x2is i-1 base 2. Represents observed phenotype*/
/* y is 8x2, z is 64x4, code is 4x4, prpheno0 is 16, prpheno1 is 16x8,
prpheno2 is 16x64, mpp is 4, denom is 16, tmp is 16, tmpt is 16, ltmp is
16 */
/* The values of x, y, z, and code are unchanging and are prespecified
outside this subroutine. Can think of the prpheno's, mpp, denom, tmp,
tmpt, ltmp as temporary storage passed to this routine.*/
void get2locexprecomb(double *data,double *s, double *p, double *am,
double *co, double *mp, double *output, int **x, int **y, int **z,
int **code, double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom, double *tmp, double *tmpt, double *ltmp)
{
int i,j,k;
/*prpheno0[i] is prob. of (A,B,a,b)^(x[i,]) and no sperm present*/
/*printf("prpheno0\n");*/
for(i=0;i<16;i++){
prpheno0[i]=(1.-p[0]-p[1]);
for(j=0;j<4;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
pow(1.-co[j],1.-(double)x[i][j]));
/*printf("%e ",prpheno0[i]);*/
}
mpp[0]=mp[1]; mpp[1]=mp[0]; mpp[2]=mp[3]; mpp[3]=mp[2];
/*printf("\nmpp\n");
for(i=0;i<4;i++) printf("%lf ",mpp[i]);*/
/* Now find prob. of each phenotype event and each genotype involving
1 sperm present. prpheno1 is 16x8 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix). */
/*printf("\nprpheno1\n");*/
for(k=0;k<4;k++) for(i=0;i<16;i++){
prpheno1[i][k]=p[0]**s*mp[k];
prpheno1[i][4+k]=p[0]*(1.-*s)*mpp[k];
}
for(k=0;k<8;k++) for(i=0;i<16;i++) for(j=0;j<4;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));
/*for(i=0;i<16;i++){
for(k=0;k<8;k++) printf("%e ",prpheno1[i][k]);
printf("\n");
}*/
/* Now find prob. of each phenotype event and each genotype involving
2 sperm present. prpheno2 is 16x64 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix), 4th digit is 2nd sperm: same */
for(k=0;k<64;k++) for(i=0;i<16;i++){
prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(mp[z[k][2]],
1.-(double)z[k][0])*pow(mpp[z[k][2]],(double)z[k][0])*
pow(mp[z[k][3]],1.-(double)z[k][1])*
pow(mpp[z[k][3]],(double)z[k][1]);
for(j=0;j<4;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
(double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}
/*printf("\nprpheno2\n");
for(i=0;i<16;i++){
for(j=0;j<64;j++) printf("%e ",prpheno2[i][j]);
printf("\n");
}*/
/*printf("denom\n");*/
for(i=0;i<16;i++){
denom[i]=prpheno0[i];
for(j=0;j<8;j++) denom[i]+=prpheno1[i][j];
for(j=0;j<64;j++) denom[i]+=prpheno2[i][j];
/*printf("%e ",denom[i]);*/
}
/* number of segregations of 1 (mutant), of 2 (normal) */
output[0]=0.;
for(i=0;i<16;i++){
for(j=0;j<4;j++) if(denom[i]>0.) output[0]+=(data[i]/denom[i]*
prpheno1[i][j]);
for(j=0;j<16;j++) if(denom[i]>0.) output[0]+=(data[i]/denom[i]*2.*
prpheno2[i][j]);
for(j=16;j<48;j++) if(denom[i]>0.) output[0]+=(data[i]/denom[i]*
prpheno2[i][j]);
}
output[1]=0.;
for(i=0;i<16;i++){
for(j=4;j<8;j++) if(denom[i]>0.) output[1]+=(data[i]/denom[i]*
prpheno1[i][j]);
for(j=48;j<64;j++) if(denom[i]>0.) output[1]+=(data[i]/denom[i]*2.*
prpheno2[i][j]);
for(j=16;j<48;j++) if(denom[i]>0.) output[1]+=(data[i]/denom[i]*
prpheno2[i][j]);
}
/* number of times 1 sperm is deposited, 2 sperm, 0 sperm */
output[2]=0.;
for(i=0;i<16;i++) for(j=0;j<8;j++) if(denom[i]>0.)
output[2]+=(data[i]/denom[i]*
prpheno1[i][j]);
output[3]=0.;
for(i=0;i<16;i++) for(j=0;j<64;j++) if(denom[i]>0.)
output[3]+=(data[i]/denom[i]*
prpheno2[i][j]);
output[4]=-output[2]-output[3];
for(i=0;i<16;i++) output[4]+=data[i];
/* number of times each allele is amplified/ is present but is not
amplified */
for(j=0;j<4;j++){
for(i=0;i<16;i++){
tmp[i]=0.;
tmpt[i]=0.;
ltmp[i]=0.;
}
for(k=0;k<8;k++){
for(i=0;i<16;i++){
tmp[i]+=((double)code[y[k][1]][j]*
prpheno1[i][k]*am[j]/(am[j]+(1.-am[j])*
co[j]));
tmpt[i]+=((double)code[y[k][1]][j]*
prpheno1[i][k]);
}
}
for(k=0;k<64;k++){
for(i=0;i<16;i++){
if(pow(1.-pow(1.-am[j],(double)(
code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j]),(double)x[i][j])*pow(pow(
1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-
(double)x[i][j])==0.) ltmp[i]=1.;
tmp[i]+=((1.-ltmp[i])*(double)(
code[z[k][2]][j]+code[z[k][3]][j])*
prpheno2[i][k]*am[j]/(
pow(1.-pow(1.-am[j],(double)(
code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j]),(double)x[i][j])*pow(pow(
1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-
(double)x[i][j])+ltmp[i]));
tmpt[i]+=((1.-ltmp[i])*(double)(
code[z[k][2]][j]+code[z[k][3]][j])*
prpheno2[i][k]);
}
}
output[3+2*(j+1)]=0.;
output[4+2*(j+1)]=0.;
for(i=0;i<16;i++){
if(x[i][j]!=1) tmp[i]=0.;
if (denom[i]>0.)
output[3+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
if(denom[i]>0.)
output[4+2*(j+1)]+=(data[i]/denom[i]*(tmpt[i]-tmp[i]));
}
}
/* number of times each allele is/is not a contaminant */
for(j=0;j<4;j++){
for(i=0;i<16;i++){
tmp[i]=prpheno0[i];
ltmp[i]=0.;
}
for(k=0;k<8;k++){
for(i=0;i<16;i++) tmp[i]+=(prpheno1[i][k]*pow((1.-
am[j])*co[j]/(am[j]+(1.-am[j])*co[j]),
(double)code[y[k][1]][j]));
}
for(k=0;k<64;k++){
for(i=0;i<16;i++){
if(1.-pow(1.-am[j],(double)(
code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j])==0.) ltmp[i]=1.;
tmp[i]+=((1.-ltmp[i])*prpheno2[i][k]*pow(1.-
am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*co[j]/(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j])+ltmp[i]));
}
}
output[11+2*(j+1)]=0.;
output[12+2*(j+1)]=0.;
for(i=0;i<16;i++){
if(x[i][j]!=1) tmp[i]=0.;
if(denom[i]>0.)
output[11+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
output[12+2*(j+1)]+=(data[i]*
(double)(1-x[i][j]));
}
}
}
void makexyzcode(int **x,int **y,int **z,int **code)
/* x is 16x4, y is 8x2, z is 64x4, code is 4x4 */
{
int i,j;
for(i=0;i<16;i++){
x[i][3]=(int)fmod((double)i,2.);
x[i][2]=(int)fmod((double)((i-x[i][3])/2),2.);
x[i][1]=(int)fmod((double)(((i-x[i][3])/2-x[i][2])/2),2.);
x[i][0]=(int)fmod((double)((((i-x[i][3])/2-x[i][2])/2-x[i][1])/
2),2.);
}
for(j=0;j<4;j++){
y[j][0]=0;
y[j+4][0]=1;
y[j][1]=j;
y[j+4][1]=j;
}
code[0][0]=1;
code[0][1]=1;
code[0][2]=0;
code[0][3]=0;
code[1][0]=0;
code[1][1]=0;
code[1][2]=1;
code[1][3]=1;
code[2][0]=1;
code[2][1]=0;
code[2][2]=0;
code[2][3]=1;
code[3][0]=0;
code[3][1]=1;
code[3][2]=1;
code[3][3]=0;
for(i=0;i<32;i++){
z[i][0]=0;
z[i+32][0]=1;
}
for(i=0;i<16;i++){
z[i][1]=0;
z[i+16][1]=1;
z[i+32][1]=0;
z[i+48][1]=1;
}
for(i=0;i<4;i++){
for(j=0;j<4;j++){
z[i+j*16][2]=0;
z[i+4+j*16][2]=1;
z[i+8+j*16][2]=2;
z[i+12+j*16][2]=3;
}
}
for(i=0;i<4;i++) for(j=0;j<16;j++) z[j*4+i][3]=i;
}
void makeonexyzcode(int **x1,int **y1,int **z1,int **code1)
/* x1 is 4x2, y1 is 4x2, z1 is 16x4 code1 is 2x2 */
{
int i,j;
for(i=0;i<4;i++){
x1[i][1]=(int)fmod((double)i,2.);
x1[i][0]=(int)fmod((double)((i-x1[i][1])/2),2.);
}
for(i=0;i<2;i++){
y1[i][0]=0;
y1[i+2][0]=1;
y1[i][1]=i;
y1[i+2][1]=i;
}
code1[0][0]=1;
code1[0][1]=0;
code1[1][0]=0;
code1[1][1]=1;
for(i=0;i<8;i++){
z1[i][0]=0;
z1[i+8][0]=1;
}
for(i=0;i<4;i++){
z1[i][1]=0;
z1[i+4][1]=1;
z1[i+8][1]=0;
z1[i+12][1]=1;
}
for(i=0;i<2;i++){
for(j=0;j<4;j++){
z1[i+j*4][2]=0;
z1[i+2+j*4][2]=1;
}
}
for(i=0;i<2;i++) for(j=0;j<8;j++) z1[j*2+i][3]=i;
}
/* The function get.1locexp.recomb() takes as input the 1-locus sperm
typing
data counts (length 4) in order -,a,A,Aa
and the current parameter values (length 7) in order s (segregation
probability for allele 1, assoc. with AB in parent) p=(p1,p2)
probabilities of 1 and 2 sperm
deposited with p0 assumed = 1-p1-p2, am vector of 2 amplification
probabilities for alleles A and a; co vector of 2 contamination
probabilities for alleles A and a. Also takes as input the
recombination fraction between the gene and the marker.
The output of get.1locexp.recomb() is the expected values of the
following
quantities (length 13), condional on the data and the current
parameter values:
number of segregations of allele 1, associated in parent with A,
number of segregations of allele 2, associated in parent with a,
number of times 1 sperm is deposited, number of times 2 sperm are
deposited, number of times 0 sperm are deposited,
number of times
allele A is amplified, number of times allele A is present but not
amplified, similarly for a, number of times allele
A is a contaminant, number of times allele A is not a contaminant,
number of times a is/is not a contaminant
For the purpose, see under get.2locexp.recomb() above.
x is 4x2, ith row is i-1 base 2. y is 4x2, z is 16x4, code is 2x2,
prpheno0 is 4, prpheno1 is 4x4, prpheno2 is 4x16, denom is 4,
tmp is 4, tmpt is 4, ltmp is 4*/
void get1locexprecomb(double *data,double *s,double *p,double *am,double
*co,double *rf,double *output,int **x,int **y,int **z, int **code,
double *prpheno0, double **prpheno1,double **prpheno2,double *denom,
double *tmp,double *tmpt,double *ltmp)
{
int i,j,k;
/* prpheno0[i] is prob of (A,a)^x[i,] and no sperm present */
/*printf("prpheno0\n");*/
for(i=0;i<4;i++){
prpheno0[i]=(1.-p[0]-p[1]);
for(j=0;j<2;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
pow(1.-co[j],1.-(double)x[i][j]));
/*printf("%e ",prpheno0[i]);
fflush(stdout);*/
}
/* Now find prob. of each phenotype event and each genotype involving
1 sperm present. prpheno1 is 4x4 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=A, 2=a
(given by code matrix).*/
/*printf("\nprpheno1\n");*/
for(i=0;i<4;i++){
prpheno1[i][0]=p[0]**s*(1.-*rf);
prpheno1[i][1]=p[0]**s**rf;
prpheno1[i][2]=p[0]*(1.-*s)**rf;
prpheno1[i][3]=p[0]*(1.-*s)*(1.-*rf);
}
for(k=0;k<4;k++) for(i=0;i<4;i++) for(j=0;j<2;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));
/*for(i=0;i<4;i++){
for(k=0;k<4;k++) printf("%e ",prpheno1[i][k]);
printf("\n");
fflush(stdout);
}*/
/* Now find prob. of each phenotype event and each genotype involving
2 sperm present. prpheno2 is 4x16 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=A, 2=a
(given by code matrix), 4th digit is 2nd sperm: same */
for(k=0;k<16;k++) for(i=0;i<4;i++){
prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(*rf,fmod(fabs((double)
(z[k][0]-z[k][2])),2.)+fmod(fabs((double)(z[k][1]-z[k][3])),2.))*
pow(1.-*rf,fmod(fabs((double)(z[k][0]-z[k][2]+1)),2.)+
fmod(fabs((double)(z[k][1]-z[k][3]+1)),2.));
for(j=0;j<2;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
(double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}
/*printf("\nprpheno2\n");
for(i=0;i<4;i++){
for(j=0;j<16;j++) printf("%e ",prpheno2[i][j]);
printf("\n");
fflush(stdout);
}*/
/*printf("denom\n");*/
for(i=0;i<4;i++){
denom[i]=prpheno0[i];
for(j=0;j<4;j++) denom[i]+=prpheno1[i][j];
for(j=0;j<16;j++) denom[i]+=prpheno2[i][j];
/*printf("%e ",denom[i]);
fflush(stdout);*/
}
/* number of segregations of 1 (mutant), of 2 (normal) */
output[0]=0.;
for(i=0;i<4;i++){
for(j=0;j<2;j++) if(denom[i]>0.)
output[0]+=(data[i]/denom[i]*
prpheno1[i][j]);
for(j=0;j<4;j++) if(denom[i]>0.)
output[0]+=(data[i]/denom[i]*2.*
prpheno2[i][j]);
for(j=4;j<12;j++) if(denom[i]>0.)
output[0]+=(data[i]/denom[i]*
prpheno2[i][j]);
}
/*printf("output0 %lf\n",output[0]);
fflush(stdout);*/
output[1]=0.;
for(i=0;i<4;i++){
for(j=2;j<4;j++) if(denom[i]>0.)
output[1]+=(data[i]/denom[i]*
prpheno1[i][j]);
for(j=12;j<16;j++) if(denom[i]>0.)
output[1]+=(data[i]/denom[i]*2.*
prpheno2[i][j]);
for(j=4;j<12;j++) if(denom[i]>0.)
output[1]+=(data[i]/denom[i]*
prpheno2[i][j]);
}
/*printf("output1 %lf\n",output[1]);
fflush(stdout);*/
/* number of times 1 sperm is deposited, 2 sperm, 0 sperm */
output[2]=0.;
for(i=0;i<4;i++) for(j=0;j<4;j++) if(denom[i]>0.)
output[2]+=(data[i]/denom[i]*
prpheno1[i][j]);
/*printf("output2 %lf\n",output[2]);
fflush(stdout);*/
output[3]=0.;
for(i=0;i<4;i++) for(j=0;j<16;j++) if(denom[i]>0.)
output[3]+=(data[i]/denom[i]*
prpheno2[i][j]);
/*printf("output3 %lf\n",output[3]);
fflush(stdout);*/
output[4]=-output[2]-output[3];
for(i=0;i<4;i++) output[4]+=data[i];
/*printf("output4 %lf\n",output[4]);
fflush(stdout);*/
/* number of times each allele is amplified/ is present but is not
amplified */
for(j=0;j<2;j++){
for(i=0;i<4;i++){
tmp[i]=0.;
tmpt[i]=0.;
ltmp[i]=0.;
}
for(k=0;k<4;k++){
if(am[j]+(1.-am[j])*co[j]==0.) *ltmp=1.; else *ltmp=0.;
for(i=0;i<4;i++){
tmp[i]+=((1.-*ltmp)*(double)code[y[k][1]][j]*
prpheno1[i][k]*am[j]/(am[j]+(1.-am[j])*
co[j]+*ltmp));
tmpt[i]+=((double)code[y[k][1]][j]*
prpheno1[i][k]);
}
}
ltmp[0]=0.;
for(k=0;k<16;k++){
for(i=0;i<4;i++){
if(pow(1.-pow(1.-am[j],(double)(
code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j]),(double)x[i][j])*pow(pow(
1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-
(double)x[i][j])==0.) ltmp[i]=1.;
else ltmp[i]=0.;
tmp[i]+=((1.-ltmp[i])*(double)(
code[z[k][2]][j]+code[z[k][3]][j])*
prpheno2[i][k]*am[j]/(
pow(1.-pow(1.-am[j],(double)(
code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j]),(double)x[i][j])*pow(pow(
1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-
(double)x[i][j])+ltmp[i]));
tmpt[i]+=((double)(
code[z[k][2]][j]+code[z[k][3]][j])*
prpheno2[i][k]);
}
}
output[3+2*(j+1)]=0.;
output[4+2*(j+1)]=0.;
for(i=0;i<4;i++){
if(x[i][j]!=1) tmp[i]=0.;
if(denom[i]>0.)
output[3+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
if(denom[i]>0.)
output[4+2*(j+1)]+=(data[i]/denom[i]*(tmpt[i]-tmp[i]));
}
}
/*printf("output5 %lf\n",output[5]);
printf("output6 %lf\n",output[6]);
printf("output7 %lf\n",output[7]);
printf("output8 %lf\n",output[8]);
fflush(stdout);*/
/* number of times each allele is/is not a contaminant */
for(j=0;j<2;j++){
for(i=0;i<4;i++){
tmp[i]=prpheno0[i];
ltmp[i]=0.;
}
for(k=0;k<4;k++){
if(pow(am[j]+(1.-am[j])*co[j],
(double)code[y[k][1]][j])==0.) *ltmp=1.;
else *ltmp=0.;
for(i=0;i<4;i++) tmp[i]+=((1.-*ltmp)*
prpheno1[i][k]*pow((1.-
am[j])*co[j]/(am[j]+(1.-am[j])*co[j]+*ltmp),
(double)code[y[k][1]][j]));
}
*ltmp=0.;
for(k=0;k<16;k++){
for(i=0;i<4;i++){
if(1.-pow(1.-am[j],(double)(
code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j])==0.) ltmp[i]=1.; else ltmp[i]=0.;
tmp[i]+=((1.-ltmp[i])*prpheno2[i][k]*pow(1.-
am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*co[j]/(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*
(1.-co[j])+ltmp[i]));
}
}
output[7+2*(j+1)]=0.;
output[8+2*(j+1)]=0.;
for(i=0;i<4;i++){
if(x[i][j]!=1) tmp[i]=0.;
if(denom[i]>0.)
output[7+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
output[8+2*(j+1)]+=(data[i]*
(double)(1-x[i][j]));
}
}
}
/* The function get.2loglike.recomb() takes as input the 2-locus
sperm-typing
data counts (length 16) in order given above and the current
parameter values (length 11) in order given above. Also takes as input
the
multilocus recombination probabilities: Pr(no rec., i.e. 1AB or 2ab),
Pr(rec btw disease and marker haplotype, i.e. 1ab or 2AB), Pr(1Ab or
2aB), Pr(1aB or 2Ab), where parental haplotypes are 1AB and 2ab.
The output of get.2loglike.recomb() is the log-likelihood of the data
given
the parameter values.*/
/* x is 16x4, ith row of x is i-1 base 2.*/
void get2loglikerecomb(double *data,double *s, double *p,double *am,double
*co,double *mp,double *loglike,int **x, int **y, int **z,
int **code, double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom)
{
int i,j,k;
/*prpheno0[i] is prob. of (A,B,a,b)^(x[i,]) and no sperm present*/
/*printf("prpheno0\n");*/
for(i=0;i<16;i++){
prpheno0[i]=(1.-p[0]-p[1]);
for(j=0;j<4;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
pow(1.-co[j],1.-(double)x[i][j]));
/*printf("%e ",prpheno0[i]);*/
}
mpp[0]=mp[1]; mpp[1]=mp[0]; mpp[2]=mp[3]; mpp[3]=mp[2];
/*printf("\nmpp\n");
for(i=0;i<4;i++) printf("%lf ",mpp[i]);*/
/* Now find prob. of each phenotype event and each genotype involving
1 sperm present. prpheno1 is 16x8 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix). */
/*printf("\nprpheno1\n");*/
for(k=0;k<4;k++) for(i=0;i<16;i++){
prpheno1[i][k]=p[0]**s*mp[k];
prpheno1[i][4+k]=p[0]*(1.-*s)*mpp[k];
}
for(k=0;k<8;k++) for(i=0;i<16;i++) for(j=0;j<4;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));
/*for(i=0;i<16;i++){
for(k=0;k<8;k++) printf("%e ",prpheno1[i][k]);
printf("\n");
}*/
/* Now find prob. of each phenotype event and each genotype involving
2 sperm present. prpheno2 is 16x64 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix), 4th digit is 2nd sperm: same */
for(k=0;k<64;k++) for(i=0;i<16;i++){
prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(mp[z[k][2]],
1.-(double)z[k][0])*pow(mpp[z[k][2]],(double)z[k][0])*
pow(mp[z[k][3]],1.-(double)z[k][1])*
pow(mpp[z[k][3]],(double)z[k][1]);
for(j=0;j<4;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
(double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}
/*printf("\nprpheno2\n");
for(i=0;i<16;i++){
for(j=0;j<64;j++) printf("%e ",prpheno2[i][j]);
printf("\n");
}*/
/*printf("denom\n");*/
for(i=0;i<16;i++){
denom[i]=prpheno0[i];
for(j=0;j<8;j++) denom[i]+=prpheno1[i][j];
for(j=0;j<64;j++) denom[i]+=prpheno2[i][j];
/*printf("%e ",denom[i]);*/
}
*loglike=0.;
for(i=0;i<16;i++) if(denom[i]>0.) *loglike+=(data[i]*log(denom[i]));
}
/*The function get.1loglike.recomb() takes as input the 1-locus sperm
typing
data counts (length 4) in order given above and the current
parameter values (length 7) in order given above. Also takes as input
the
recombination fraction between the gene and the marker.
The output of get.1loglike.recomb() is the log-likelihood of the data
given
the parameter values.*/
void get1loglikerecomb(double *data,double *s,double *p,double *am,double
*co,double *rf,double *loglike,int **x,int **y,int **z, int **code,
double *prpheno0, double **prpheno1,double **prpheno2,double *denom)
{
int i,j,k;
/* prpheno0[i] is prob of (A,a)^x[i,] and no sperm present */
/*printf("prpheno0\n");*/
for(i=0;i<4;i++){
prpheno0[i]=(1.-p[0]-p[1]);
for(j=0;j<2;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
pow(1.-co[j],1.-(double)x[i][j]));
/*printf("%e ",prpheno0[i]);
fflush(stdout);*/
}
/* Now find prob. of each phenotype event and each genotype involving
1 sperm present. prpheno1 is 4x4 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=A, 2=a
(given by code matrix).*/
/*printf("\nprpheno1\n");*/
for(i=0;i<4;i++){
prpheno1[i][0]=p[0]**s*(1.-*rf);
prpheno1[i][1]=p[0]**s**rf;
prpheno1[i][2]=p[0]*(1.-*s)**rf;
prpheno1[i][3]=p[0]*(1.-*s)*(1.-*rf);
}
for(k=0;k<4;k++) for(i=0;i<4;i++) for(j=0;j<2;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));
/*for(i=0;i<4;i++){
for(k=0;k<4;k++) printf("%e ",prpheno1[i][k]);
printf("\n");
fflush(stdout);
}*/
/* Now find prob. of each phenotype event and each genotype involving
2 sperm present. prpheno2 is 4x16 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=A, 2=a
(given by code matrix), 4th digit is 2nd sperm: same */
for(k=0;k<16;k++) for(i=0;i<4;i++){
prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(*rf,fmod(fabs((double)
(z[k][0]-z[k][2])),2.)+fmod(fabs((double)(z[k][1]-z[k][3])),2.))*
pow(1.-*rf,fmod(fabs((double)(z[k][0]-z[k][2]+1)),2.)+
fmod(fabs((double)(z[k][1]-z[k][3]+1)),2.));
for(j=0;j<2;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
(double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}
/*printf("\nprpheno2\n");
for(i=0;i<4;i++){
for(j=0;j<16;j++) printf("%e ",prpheno2[i][j]);
printf("\n");
fflush(stdout);
}*/
/*printf("denom\n");*/
for(i=0;i<4;i++){
denom[i]=prpheno0[i];
for(j=0;j<4;j++) denom[i]+=prpheno1[i][j];
for(j=0;j<16;j++) denom[i]+=prpheno2[i][j];
/*printf("%e ",denom[i]);
fflush(stdout);*/
}
*loglike=0.;
for(i=0;i<4;i++) if(denom[i]>0.) *loglike+=(data[i]*log(denom[i]));
}
/* The function get2expcellprobs() takes as input the current
parameter values (length 11) in order given above. Also takes as input
the
multilocus recombination probabilities: Pr(no rec., i.e. 1AB or 2ab),
Pr(rec btw disease and marker haplotype, i.e. 1ab or 2AB), Pr(1Ab or
2aB), Pr(1aB or 2Ab), where parental haplotypes are 1AB and 2ab.
The output of get2expcellprobs() is the cell probabilities for
the data given the parameter values.*/
/* x is 16x4, ith row of x is i-1 base 2.*/
void get2expcellprobs(double *s, double *p,double *am,double
*co,double *mp,int **x, int **y, int **z, int **code,
double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom,double *expcellprobs)
{
int i,j,k;
/*prpheno0[i] is prob. of (A,B,a,b)^(x[i,]) and no sperm present*/
/*printf("prpheno0\n");*/
for(i=0;i<16;i++){
prpheno0[i]=(1.-p[0]-p[1]);
for(j=0;j<4;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
pow(1.-co[j],1.-(double)x[i][j]));
/*printf("%e ",prpheno0[i]);*/
}
mpp[0]=mp[1]; mpp[1]=mp[0]; mpp[2]=mp[3]; mpp[3]=mp[2];
/*printf("\nmpp\n");
for(i=0;i<4;i++) printf("%lf ",mpp[i]);*/
/* Now find prob. of each phenotype event and each genotype involving
1 sperm present. prpheno1 is 16x8 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix). */
/*printf("\nprpheno1\n");*/
for(k=0;k<4;k++) for(i=0;i<16;i++){
prpheno1[i][k]=p[0]**s*mp[k];
prpheno1[i][4+k]=p[0]*(1.-*s)*mpp[k];
}
for(k=0;k<8;k++) for(i=0;i<16;i++) for(j=0;j<4;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));
/*for(i=0;i<16;i++){
for(k=0;k<8;k++) printf("%e ",prpheno1[i][k]);
printf("\n");
}*/
/* Now find prob. of each phenotype event and each genotype involving
2 sperm present. prpheno2 is 16x64 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix), 4th digit is 2nd sperm: same */
for(k=0;k<64;k++) for(i=0;i<16;i++){
prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(mp[z[k][2]],
1.-(double)z[k][0])*pow(mpp[z[k][2]],(double)z[k][0])*
pow(mp[z[k][3]],1.-(double)z[k][1])*
pow(mpp[z[k][3]],(double)z[k][1]);
for(j=0;j<4;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
(double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}
/*printf("\nprpheno2\n");
for(i=0;i<16;i++){
for(j=0;j<64;j++) printf("%e ",prpheno2[i][j]);
printf("\n");
}*/
/*printf("denom\n");*/
for(i=0;i<16;i++){
denom[i]=prpheno0[i];
for(j=0;j<8;j++) denom[i]+=prpheno1[i][j];
for(j=0;j<64;j++) denom[i]+=prpheno2[i][j];
/*printf("%e ",denom[i]);*/
}
for(i=0;i<16;i++) expcellprobs[i]=denom[i];
}
/* The function get1expcellprobs() takes as input the current
parameter values (length 7) in order given above. Also takes as input
the recombination fraction between the gene and the marker.
The output of get1expcellprobs() is the cell probabilities for
the data given the parameter values.*/
void get1expcellprobs(double *s, double *p,double *am,double
*co,double *rf,int **x, int **y, int **z, int **code,
double *prpheno0, double **prpheno1, double **prpheno2,
double *denom,double *expcellprobs)
{
int i,j,k;
/* prpheno0[i] is prob of (A,a)^x[i,] and no sperm present */
/*printf("prpheno0\n");*/
for(i=0;i<4;i++){
prpheno0[i]=(1.-p[0]-p[1]);
for(j=0;j<2;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
pow(1.-co[j],1.-(double)x[i][j]));
/*printf("%e ",prpheno0[i]);
fflush(stdout);*/
}
/* Now find prob. of each phenotype event and each genotype involving
1 sperm present. prpheno1 is 4x4 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=A, 2=a
(given by code matrix).*/
/*printf("\nprpheno1\n");*/
for(i=0;i<4;i++){
prpheno1[i][0]=p[0]**s*(1.-*rf);
prpheno1[i][1]=p[0]**s**rf;
prpheno1[i][2]=p[0]*(1.-*s)**rf;
prpheno1[i][3]=p[0]*(1.-*s)*(1.-*rf);
}
for(k=0;k<4;k++) for(i=0;i<4;i++) for(j=0;j<2;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));
/*for(i=0;i<4;i++){
for(k=0;k<4;k++) printf("%e ",prpheno1[i][k]);
printf("\n");
fflush(stdout);
}*/
/* Now find prob. of each phenotype event and each genotype involving
2 sperm present. prpheno2 is 4x16 matrix, ith row corresponds to
phenotype x[i,]. kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=A, 2=a
(given by code matrix), 4th digit is 2nd sperm: same */
for(k=0;k<16;k++) for(i=0;i<4;i++){
prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(*rf,fmod(fabs((double)
(z[k][0]-z[k][2])),2.)+fmod(fabs((double)(z[k][1]-z[k][3])),2.))*
pow(1.-*rf,fmod(fabs((double)(z[k][0]-z[k][2]+1)),2.)+
fmod(fabs((double)(z[k][1]-z[k][3]+1)),2.));
for(j=0;j<2;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
(double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
(double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}
/*printf("\nprpheno2\n");
for(i=0;i<4;i++){
for(j=0;j<16;j++) printf("%e ",prpheno2[i][j]);
printf("\n");
fflush(stdout);
}*/
/*printf("denom\n");*/
for(i=0;i<4;i++){
denom[i]=prpheno0[i];
for(j=0;j<4;j++) denom[i]+=prpheno1[i][j];
for(j=0;j<16;j++) denom[i]+=prpheno2[i][j];
/*printf("%e ",denom[i]);
fflush(stdout);*/
}
for(i=0;i<4;i++) expcellprobs[i]=denom[i];
}
#include <malloc.h>
#include <stdio.h>
void nrerror(error_text)
char error_text[];
{
void exit();
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
double *vector(nl,nh)
int nl,nh;
{
double *v;
v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
if (!v) nrerror("allocation failure in vector()");
return v-nl;
}
int *ivector(nl,nh)
int nl,nh;
{
int *v;
v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
if (!v) nrerror("allocation failure in ivector()");
return v-nl;
}
double *dvector(nl,nh)
int nl,nh;
{
double *v;
v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
if (!v) nrerror("allocation failure in dvector()");
return v-nl;
}
double **matrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
int i;
double **m;
m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
if (!m) nrerror("allocation failure 1 in matrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
if (!m[i]) nrerror("allocation failure 2 in matrix()");
m[i] -= ncl;
}
return m;
}
double **dmatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
int i;
double **m;
m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
if (!m) nrerror("allocation failure 1 in dmatrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
m[i] -= ncl;
}
return m;
}
int **imatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
int i,**m;
m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
if (!m) nrerror("allocation failure 1 in imatrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
if (!m[i]) nrerror("allocation failure 2 in imatrix()");
m[i] -= ncl;
}
return m;
}
double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
double **a;
int oldrl,oldrh,oldcl,oldch,newrl,newcl;
{
int i,j;
double **m;
m=(double **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(double*));
if (!m) nrerror("allocation failure in submatrix()");
m -= newrl;
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;
return m;
}
void free_vector(v,nl,nh)
double *v;
int nl,nh;
{
free((char*) (v+nl));
}
void free_ivector(v,nl,nh)
int *v,nl,nh;
{
free((char*) (v+nl));
}
void free_dvector(v,nl,nh)
double *v;
int nl,nh;
{
free((char*) (v+nl));
}
void free_matrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
{
int i;
for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
free((char*) (m+nrl));
}
void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
{
int i;
for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
free((char*) (m+nrl));
}
void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
int nrl,nrh,ncl,nch;
{
int i;
for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
free((char*) (m+nrl));
}
void free_submatrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
free((char*) (b+nrl));
}
double **convert_matrix(a,nrl,nrh,ncl,nch)
double *a;
int nrl,nrh,ncl,nch;
{
int i,j,nrow,ncol;
double **m;
nrow=nrh-nrl+1;
ncol=nch-ncl+1;
m = (double **) malloc((unsigned) (nrow)*sizeof(double*));
if (!m) nrerror("allocation failure in convert_matrix()");
m -= nrl;
for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl;
return m;
}
void free_convert_matrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
free((char*) (b+nrl));
}
#include <math.h>
#include <stdio.h>
/* sfix is a vector of length ndon indicating which of the s variables
are
to be fixed (1 means fix it, 0 means let it vary). The
corresponding
values are given in sval (vector of length ndon). The entry for a
variable will be ignored if that variable is not fixed. pfix is a
ndonx2 matrix, similarly for pval. amfix, amval,
cofix, coval are vectors of length nal. sstart, pstar, amstart,
costart are of the same type as s, p, am, and co, and contain the
starting values for these parameters. Now the probabilities of 1 sperm
deposited and the probabilities of 2 sperm deposited can have equality
relationships that are separately specified. Also have changed the way
the EM algorithm decides when it's done.*/
void spermanalysisrfssdem(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *sstart,double **pstart,double
*amstart,double *costart,double *dconv,int *niter,int simind,int
**ptype,double *dtmp)
{
int i,j,k,l, iconv;
void makexyzcode(int **x,int **y,int **z,int **code);
void makeonexyzcode(int **x1,int **y1,int **z1,int **code1);
void get2locexprecomb(double *data,double *s, double *p, double *am,
double *co, double *mp, double *output, int **x, int **y, int **z,
int **code, double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom, double *tmp, double *tmpt, double *ltmp);
void get1locexprecomb(double *data,double *s,double *p,double *am,double
*co,double *rf,double *output,int **x,int **y,int **z, int **code,
double *prpheno0, double **prpheno1,double **prpheno2,double *denom,
double *tmp,double *tmpt,double *ltmp);
void get2loglikerecomb(double *data,double *s, double *p, double *am,
double *co,double *mp,double *loglike,int **x, int **y, int **z,
int **code, double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom);
void get1loglikerecomb(double *data,double *s,double *p,double *am,
double *co,double *rf,double *loglike,int **x,int **y,int **z,
int **code, double *prpheno0, double **prpheno1,double **prpheno2,
double *denom);
makexyzcode(x,y,z,code);
makeonexyzcode(x1,y1,z1,code1);
for(i=0;i<ndon;i++){
if(sfix[i]==1) s[i]=sval[i]; else s[i]=sstart[i];
if(pfix[i][0]==1) p[i][0]=pval[i][0];
else p[i][0]=pstart[i][0];
if(pfix[i][1]==1) p[i][1]=pval[i][1];
else p[i][1]=pstart[i][1];
}
/* Amplification an contamination params are in order A,B,a,b */
for(i=0;i<nal;i++){
if(amfix[i]==1) am[i]=amval[i]; else am[i]=amstart[i];
if(cofix[i]==1) co[i]=coval[i]; else co[i]=costart[i];
}
iconv=0;
*niter=0;
while(iconv==0){
*niter+=300;
for(k=0;k<300;k++){
/*printf("d2 1 %f\n",p[1][1]);
fflush(stdout);*/
for(i=0;i<n2loc;i++){
for(j=0;j<4;j++) amin[j]=am[amal2[i][j]];
for(j=0;j<4;j++) coin[j]=co[coal2[i][j]];
for(j=0;j<2;j++) pin[j]=p[pind2[i][j]][j];
/*if(*niter==2400){
printf("iteration %d\n",k);
printf("twoloc\n");
for(j=0;j<16;j++) printf("%f ",twoloc[i][j]);
printf("\ns p1 p2 a1 a2 a3 a4 c1 c2 c3 c4\n");*/
/*printf("%f %f %f %f %f %f %f %f %f %f %f\n",s[sind2[i]],
pin[0],pin[1],amin[0],amin[1],amin[2],amin[3],coin[0],coin[1],
coin[2],coin[3]);*/
/*printf("twolocmp %f %f %f %f\n",twolocmp[i][0],twolocmp[i][1],
twolocmp[i][2],twolocmp[i][3]);
}*/
get2locexprecomb(twoloc[i],s+sind2[i],pin,
amin,coin,twolocmp[i],output,x,y,z,code,prpheno0,
prpheno1,prpheno2,mpp,denom,tmp,tmpt,ltmp);
/*if(*niter==2400){
printf("Output of get2locexprecomb\n");
for(j=0;j<21;j++) printf("%f ",output[j]);
printf("\n");
fflush(stdout);
}*/
for(j=0;j<21;j++) exp2[i][j]=output[j];
}
for(i=0;i<n1loc;i++){
for(j=0;j<2;j++) amin[j]=am[amal1[i][j]];
for(j=0;j<2;j++) coin[j]=co[coal1[i][j]];
for(j=0;j<2;j++) pin[j]=p[pind1[i][j]][j];
get1locexprecomb(oneloc[i],s+sind1[i],pin,
amin,coin,onelocrf+i,output,x1,y1,z1,code1,prpheno0,
prpheno11,prpheno21,denom,tmp,tmpt,ltmp);
for(j=0;j<13;j++) exp1[i][j]=output[j];
}
for(i=0;i<ndon;i++){
s[i]=0.;
sden[i]=0.;
p[i][0]=0.;
pden[i][0]=0.;
p[i][1]=0.;
pden[i][1]=0.;
}
for(i=0;i<nal;i++){
am[i]=0.;
amden[i]=0.;
co[i]=0.;
coden[i]=0.;
}
for(j=0;j<n2loc;j++){
s[sind2[j]]+=exp2[j][0];
sden[sind2[j]]+=(exp2[j][0]+exp2[j][1]);
p[pind2[j][0]][0]+=exp2[j][2];
p[pind2[j][1]][1]+=exp2[j][3];
for(i=0;i<2;i++) for(l=2;l<5;l++)
pden[pind2[j][i]][i]+=exp2[j][l];
am[amal2[j][0]]+=exp2[j][5];
am[amal2[j][1]]+=exp2[j][7];
am[amal2[j][2]]+=exp2[j][9];
am[amal2[j][3]]+=exp2[j][11];
amden[amal2[j][0]]+=(exp2[j][5]+exp2[j][6]);
amden[amal2[j][1]]+=(exp2[j][7]+exp2[j][8]);
amden[amal2[j][2]]+=(exp2[j][9]+exp2[j][10]);
amden[amal2[j][3]]+=(exp2[j][11]+exp2[j][12]);
co[coal2[j][0]]+=exp2[j][13];
co[coal2[j][1]]+=exp2[j][15];
co[coal2[j][2]]+=exp2[j][17];
co[coal2[j][3]]+=exp2[j][19];
coden[coal2[j][0]]+=(exp2[j][13]+exp2[j][14]);
coden[coal2[j][1]]+=(exp2[j][15]+exp2[j][16]);
coden[coal2[j][2]]+=(exp2[j][17]+exp2[j][18]);
coden[coal2[j][3]]+=(exp2[j][19]+exp2[j][20]);
}
for(j=0;j<n1loc;j++){
s[sind1[j]]+=exp1[j][0];
sden[sind1[j]]+=(exp1[j][0]+exp1[j][1]);
p[pind1[j][0]][0]+=exp1[j][2];
p[pind1[j][1]][1]+=exp1[j][3];
for(i=0;i<2;i++) for(l=2;l<5;l++)
pden[pind1[j][i]][i]+=exp1[j][l];
am[amal1[j][0]]+=exp1[j][5];
am[amal1[j][1]]+=exp1[j][7];
amden[amal1[j][0]]+=(exp1[j][5]+exp1[j][6]);
amden[amal1[j][1]]+=(exp1[j][7]+exp1[j][8]);
/*printf("%lf %lf\n",co[coal1[j][0]],co[coal1[j][1]]);*/
co[coal1[j][0]]+=exp1[j][9];
co[coal1[j][1]]+=exp1[j][11];
/*printf("%lf %lf\n",co[coal1[j][0]],co[coal1[j][1]]);*/
/*printf("%lf %lf\n",coden[coal1[j][0]],coden[coal1[j][1]]);*/
coden[coal1[j][0]]+=(exp1[j][9]+exp1[j][10]);
coden[coal1[j][1]]+=(exp1[j][11]+exp1[j][12]);
/*printf("%lf %lf\n",coden[coal1[j][0]],coden[coal1[j][1]]);*/
}
for(j=0;j<ndon;j++){
if(sden[j]>0.) s[j]=s[j]/sden[j];
if(sfix[j]==1) s[j]=sval[j];
}
for(i=0;i<ndon;i++) for(j=0;j<2;j++){
if(ptype[i][4*j+2]==1){
/*printf("i=%d, ptype[i][4*j+2]==1\n",i);*/
p[i][j]=pval[i][j];
}
else if(ptype[i][4*j+1]==1 && ptype[i][4*j+3]==0){
/*printf("i=%d, ptype[i][4*j+1]==1 && ptype[i][4*j+3]==0\n",i);*/
if(pden[i][j]>0.) p[i][j]=p[i][j]/pden[i][j];
}
else if(ptype[i][4*j+1]==1 && ptype[i][4*j+3]==1){
dtmp[0]=pden[i][j];
dtmp[2]=p[i][j];
dtmp[3]=pden[i][j];
dtmp[4]=p[i][j];
for(l=0;l<n2loc && !(pind2[l][j]==i &&
pfix[pind2[l][1-j]][1-j]==1);l++);
if(l<n2loc){
dtmp[5]=pval[pind2[l][1-j]][1-j];
dtmp[2]*=(1.-dtmp[5]);
}
else{
for(l=0;l<n1loc && !(pind1[l][j]==i &&
pfix[pind1[l][1-j]][1-j]==1);l++);
if(l<n1loc){
dtmp[5]=pval[pind1[l][1-j]][1-j];
dtmp[2]*=(1.-dtmp[5]);
}
}
for(l=0;l<n2loc;l++) if(pind2[l][j]==i &&
pfix[pind2[l][1-j]][1-j]==1){
dtmp[0]-=(exp2[l][3-j]);
dtmp[3]-=(exp2[l][3-j]+exp2[l][4]);
dtmp[4]+=(exp2[l][4]);
}
for(l=0;l<n1loc;l++) if(pind1[l][j]==i &&
pfix[pind1[l][1-j]][1-j]==1){
dtmp[0]-=(exp1[l][3-j]);
dtmp[3]-=(exp1[l][3-j]+exp1[l][4]);
dtmp[4]+=(exp1[l][4]);
}
dtmp[1]=-(1.-dtmp[5])*dtmp[3]-dtmp[4];
if(dtmp[0]<=0.) p[i][j]=0.;
else if((dtmp[3]=pow(dtmp[1],2.)-4.*dtmp[0]*
dtmp[2])<0.) p[i][j]=1.-dtmp[5];
else p[i][j]=(-dtmp[1]-pow(dtmp[3],.5))/
(2.*dtmp[0]);
/*printf("i=%d, ptype[i][4*j+1]==1 && ptype[i][4*j+3]==1\n",i);*/
}
}
for(i=0;i<ndon;i++) for(j=0;j<2;j++){
if(ptype[i][4*j+1]==0 && ptype[i][4*j]==1 &&
ptype[i][4*j+2]==0){
/*printf("i=%d, ptype[i][4*j+1]==0 && ptype[i][4*j]==1 && ptype[i][4*j+2]==0\n",
i);*/
p[i][j]=0.;
for(l=0;l<n2loc;l++) if(pind2[l][j]==i){
p[i][j]+=(exp2[l][j+2]*(1.-
p[pind2[l][1-j]][1-j]));
pden[i][j]-=(exp2[l][3-j]);
/*printf("p[i][j] %f\n",p[i][j]);
printf("pden[i][j] %f\n",pden[i][j]);*/
}
for(l=0;l<n1loc;l++) if(pind1[l][j]==i){
p[i][j]+=(exp1[l][j+2]*(1.-
p[pind1[l][1-j]][1-j]));
pden[i][j]-=(exp1[l][3-j]);
/*printf("p[i][j] %f\n",p[i][j]);
printf("pden[i][j] %f\n",pden[i][j]);*/
}
if(pden[i][j]>0.) p[i][j]=p[i][j]/pden[i][j];
}
}
for(j=0;j<nal;j++){
if(amden[j]>0.) am[j]=am[j]/amden[j];
if(coden[j]>0.) co[j]=co[j]/coden[j];
if(amfix[j]==1) am[j]=amval[j];
if(cofix[j]==1) co[j]=coval[j];
}
/*for(i=0;i<ndon;i++) printf("%7e ",s[i]);
printf("\n");
for(i=0;i<ndon;i++){
for(j=0;j<2;j++) printf("%7e ",p[i][j]);
printf("\n");
}
for(i=0;i<nal;i++) printf("%7e ",am[i]);
printf("\n");
for(i=0;i<nal;i++) printf("%7e ",co[i]);
printf("\n\n");*/
/* compare the parameters after 250th and 300th iterations and look
for a change. */
if(k==249){
j=0;
for(i=0;i<ndon;i++){
dconv[j]=s[i];
j++;
}
for(i=0;i<ndon;i++) for(l=0;l<2;l++){
dconv[j]=p[i][l];
j++;
}
for(i=0;i<nal;i++){
dconv[j]=am[i];
j++;
}
for(i=0;i<nal;i++){
dconv[j]=co[i];
j++;
}
*maxloglike=0.;
for(l=0;l<n2loc;l++){
for(i=0;i<4;i++) amin[i]=am[amal2[l][i]];
for(i=0;i<4;i++) coin[i]=co[coal2[l][i]];
for(i=0;i<2;i++) pin[i]=p[pind2[l][i]][i];
get2loglikerecomb(twoloc[l],s+sind2[l],pin,amin,
coin,twolocmp[l],loglike,x,y,z,code,prpheno0,prpheno1,
prpheno2,mpp,denom);
*maxloglike+=*loglike;
maxloglike[l+1]=*loglike;
}
for(i=0;i<n1loc;i++){
for(l=0;l<2;l++) amin[l]=am[amal1[i][l]];
for(l=0;l<2;l++) coin[l]=co[coal1[i][l]];
for(l=0;l<2;l++) pin[l]=p[pind1[i][l]][l];
get1loglikerecomb(oneloc[i],s+sind1[i],pin,amin,
coin,onelocrf+i,loglike,x1,y1,z1,code1,prpheno0,prpheno11,
prpheno21,denom);
*maxloglike+=*loglike;
maxloglike[n2loc+i+1]=*loglike;
}
dconv[3*ndon+2*nal]=*maxloglike;
}
} /* End of loop over number of iterations */
iconv=1;
j=0;
for(i=0;i<ndon;i++){
if(!(s[i]>pow(10.,-10.) || s[i]>dconv[j])
|| !(fabs(dconv[j]-s[i])/s[i]>pow(10.,-6.)) ){
/*||
!(fabs(dconv[j]-s[i])/s[i]>pow(10.,-3.) || fabs(dconv[j]-s[i])>
pow(10.,-6.))){*/
iconv*=1;
/*printf("%f %f\n",dconv[j],s[i]);*/
}
else{
iconv*=0;
/*printf("%d %f %f\n",*niter,dconv[j],s[i]);*/
}
j++;
}
for(i=0;i<ndon;i++) for(l=0;l<2;l++){
if(!(p[i][l]>pow(10.,-10.) || p[i][l]>dconv[j])
|| !(fabs(dconv[j]-p[i][l])/p[i][l]>pow(10.,-6.))){
/* ||
!(fabs(dconv[j]-p[i][l])/p[i][l]>pow(10.,-3.) || fabs(dconv[j]-p[i][l])>
pow(10.,-6.))){*/
iconv*=1;
/*printf("%f %f\n",dconv[j],p[i][l]);*/
}
else{
iconv*=0;
/*printf("%d %e %e\n",*niter,dconv[j],p[i][l]);*/
}
j++;
}
for(i=0;i<nal;i++){
if(!(am[i]>pow(10.,-10.) || am[i]>dconv[j])
|| !(fabs(dconv[j]-am[i])/am[i]>pow(10.,-6.))){
/* ||
!(fabs(dconv[j]-am[i])/am[i]>pow(10.,-3.) || fabs(dconv[j]-am[i])>
pow(10.,-6.))){*/
iconv*=1;
/*printf("%f %f\n",dconv[j],am[i]);*/
}
else{
iconv*=0;
/*printf("%d %f %f\n",*niter,dconv[j],am[i]);*/
}
j++;
}
for(i=0;i<nal;i++){
if(!(co[i]>pow(10.,-10.) || co[i]>dconv[j])
|| !(fabs(dconv[j]-co[i])/co[i]>pow(10.,-6.))){
/* ||
!(fabs(dconv[j]-co[i])/co[i]>pow(10.,-3.) || fabs(dconv[j]-co[i])>
pow(10.,-6.))){*/
iconv*=1;
/*printf("%f %f\n",dconv[j],co[i]);*/
}
else{
iconv*=0;
/*printf("%d %e %e\n",*niter,dconv[j],co[i]);*/
}
j++;
}
fflush(stdout);
if(iconv==0 && simind==1){
*maxloglike=0.;
for(i=0;i<n2loc;i++){
for(j=0;j<4;j++) amin[j]=am[amal2[i][j]];
for(j=0;j<4;j++) coin[j]=co[coal2[i][j]];
for(j=0;j<2;j++) pin[j]=p[pind2[i][j]][j];
get2loglikerecomb(twoloc[i],s+sind2[i],pin,amin,
coin,twolocmp[i],loglike,x,y,z,code,prpheno0,prpheno1,
prpheno2,mpp,denom);
*maxloglike+=*loglike;
maxloglike[i+1]=*loglike;
}
for(i=0;i<n1loc;i++){
for(j=0;j<2;j++) amin[j]=am[amal1[i][j]];
for(j=0;j<2;j++) coin[j]=co[coal1[i][j]];
for(j=0;j<2;j++) pin[j]=p[pind1[i][j]][j];
get1loglikerecomb(oneloc[i],s+sind1[i],pin,amin,
coin,onelocrf+i,loglike,x1,y1,z1,code1,prpheno0,prpheno11,
prpheno21,denom);
*maxloglike+=*loglike;
maxloglike[n2loc+i+1]=*loglike;
}
if(fabs(dconv[ndon*3+2*nal]-*maxloglike)<pow(10.,-4.)) iconv=1;
}
if(*niter>90000){
printf("Warning: exiting EM algorithm --- too many iterations\n");
for(i=0;i<n2loc;i++){
for(j=0;j<16;j++) printf("%f ",twoloc[i][j]);
printf("\n");
}
for(i=0;i<n1loc;i++){
for(j=0;j<4;j++) printf("%f ",oneloc[i][j]);
printf("\n");
}
for(i=0;i<ndon;i++) printf("%f ",sstart[i]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][0]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][1]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",amstart[i]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",costart[i]);
printf("\n");
break;
}
}/* end of while statement on iconv */
/* Now calculate the likelihood */
*maxloglike=0.;
for(i=0;i<n2loc;i++){
for(j=0;j<4;j++) amin[j]=am[amal2[i][j]];
for(j=0;j<4;j++) coin[j]=co[coal2[i][j]];
for(j=0;j<2;j++) pin[j]=p[pind2[i][j]][j];
get2loglikerecomb(twoloc[i],s+sind2[i],pin,amin,
coin,twolocmp[i],loglike,x,y,z,code,prpheno0,prpheno1,
prpheno2,mpp,denom);
*maxloglike+=*loglike;
maxloglike[i+1]=*loglike;
}
for(i=0;i<n1loc;i++){
for(j=0;j<2;j++) amin[j]=am[amal1[i][j]];
for(j=0;j<2;j++) coin[j]=co[coal1[i][j]];
for(j=0;j<2;j++) pin[j]=p[pind1[i][j]][j];
get1loglikerecomb(oneloc[i],s+sind1[i],pin,amin,
coin,onelocrf+i,loglike,x1,y1,z1,code1,prpheno0,prpheno11,
prpheno21,denom);
*maxloglike+=*loglike;
maxloglike[n2loc+i+1]=*loglike;
}
if(!(*maxloglike<0.)){
for(i=0;i<n2loc;i++){
for(j=0;j<16;j++) printf("%f ",twoloc[i][j]);
printf("\n");
}
for(i=0;i<n1loc;i++){
for(j=0;j<4;j++) printf("%f ",oneloc[i][j]);
printf("\n");
}
for(i=0;i<ndon;i++) printf("%f ",sstart[i]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][0]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][1]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",amstart[i]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",costart[i]);
printf("\n");
}
}
#include <math.h>
#include <stdio.h>
/* maxvals[3*nd(n+2*nal+1] is the vector of parameter values that
maximize the
likelihood, with last entry the maximized likelihood, quan is
qchisq(.95,1) for a 95% confidence interval, qchisq(.68,1) to get
standard errors */
void spermanalysisfssdioemcis(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *maxvals,double *quan,
int *sfix2,int **pfix2,int *amfix2,int *cofix2,double *sval2,
double **pval2,double *amval2,double *coval2,double *sstart,double
**pstart,double *amstart,double *costart,double *dconv,int *niter, int
nstep1,int nstep2,int nstep3,double *sstart2,double
**pstart2,double *amstart2,double *costart2,int **ptype,int **ptype2,
double *dtmp)
{
int i,j,k,l,ll,m,nstep,nt;
double eps,min,max,a,b,upp,low;
void spermanalysisrfssdem(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *sstart,double **pstart,double
*amstart,double *costart,double *dconv,int *niter,int simind,int
**ptype,double *dtmp);
eps=pow(10.,-10.);
min=eps;
max=1-eps;
for(i=0;i<ndon;i++){
sfix2[i]=sfix[i];
sval2[i]=sval[i];
pfix2[i][0]=pfix[i][0];
pval2[i][0]=pval[i][0];
pfix2[i][1]=pfix[i][1];
pval2[i][1]=pval[i][1];
}
for(i=0;i<nal;i++){
amfix2[i]=amfix[i];
amval2[i]=amval[i];
cofix2[i]=cofix[i];
coval2[i]=coval[i];
}
for(i=0;i<ndon;i++) for(j=0;j<8;j++) ptype2[i][j]=ptype[i][j];
for(i=0;i<ndon;i++){
sstart2[i]=maxvals[i];
for(j=0;j<2;j++) pstart2[i][j]=maxvals[ndon+2*i+j];
}
for(i=0;i<nal;i++){
amstart2[i]=maxvals[3*ndon+i];
costart2[i]=maxvals[3*ndon+nal+i];
}
printf("Confidence intervals for estimated parameters\n\n");
printf("Segregation parameters\n\n");
printf("Parameter Estimate Interval\n");
for(i=0;i<0;i++){
/*for(i=0;i<ndon;i++){*/
for(j=0;j<ndon;j++){
sfix2[j]=sfix[j];
sval2[j]=sval[j];
}
nt=1;
for(j=0;j<n2loc;j++) if(sind2[j]==i && sfix[i]==0) nt=0;
for(j=0;j<n1loc;j++) if(sind1[j]==i && sfix[i]==0) nt=0;
if(nt==0){
sfix2[i]=1;
a=maxvals[i];
b=max;
if(a>b) upp=a; else{
for(j=0;j<nstep1;j++){
sval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[k]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) a=(a+b)/2.;
else b=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
upp=b;
}
/*printf("%7e\n",upp);*/
a=min;
b=maxvals[i];
if(a>b) low=b; else{
for(j=0;j<nstep1;j++){
sval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) b=(a+b)/2.;
else a=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
low=a;
}
/*printf("%7e\n",low);*/
printf("s %4d %15f %9f %7f\n",i,maxvals[i],low,upp);
fflush(stdout);
}
}
for(j=0;j<ndon;j++){
sfix2[j]=sfix[j];
sval2[j]=sval[j];
}
printf("\nDeposit parameters\n\n");
printf("Parameter Estimate Interval\n");
for(i=0;i<ndon;i++){
/*for(j=0;j<2;j++) if(pfix[i][j]==0){*/
for(j=0;j<1;j++) if(pfix[i][j]==0){
nt=1;
b=max;
for(k=0;k<n2loc;k++) if(pind2[k][j]==i){
nt=0;
if(pfix[pind2[k][1-j]][1-j]==1)
if(b>1.-pval[pind2[k][1-j]][1-j]) b=1.-pval[pind2[k][1-j]][1-j];
}
for(k=0;k<n1loc;k++) if(pind1[k][j]==i){
nt=0;
if(pfix[pind1[k][1-j]][1-j]==1)
if(b>1.-pval[pind1[k][1-j]][1-j]) b=1.-pval[pind1[k][1-j]][1-j];
}
if(nt==0){
if(j==0) nstep=nstep1;
if(j==1) nstep=nstep2;
for(k=0;k<ndon;k++) for(l=0;l<2;l++){
pfix2[k][l]=pfix[k][l];
pval2[k][l]=pval[k][l];
ptype2[k][l*4+2]=ptype[k][l*4+2];
ptype2[k][l*4+3]=ptype[k][l*4+3];
}
pfix2[i][j]=1;
ptype2[i][j*4+2]=1;
for(k=0;k<n2loc;k++) if(pind2[k][j]==i){
if(ptype[pind2[k][1-j]][4*(1-j)+3]==1){
printf("Sorry, confidence interval calculation for this parameter\n");
printf("not yet implemented (deposit parameter d%d %d paired with\n",
j+1,i);
printf("d%d %d which is paired with another d%d parameter that is\n",
2-j,pind2[k][1-j],j+1);
printf("set to a fixed value.\n");
nt=1;
}
ptype2[pind2[k][1-j]][4*(1-j)+3]=1;
}
for(k=0;k<n1loc;k++) if(pind1[k][j]==i){
if(ptype[pind1[k][1-j]][4*(1-j)+3]==1){
printf("Sorry, confidence interval calculation for this parameter\n");
printf("not yet implemented (deposit parameter d%d %d paired with\n",
j+1,i);
printf("d%d %d which is paired with another d%d parameter that is\n",
2-j,pind1[k][1-j],j+1);
printf("set to a fixed value.\n");
nt=1;
}
ptype2[pind1[k][1-j]][4*(1-j)+3]=1;
}
if(nt==0){
a=maxvals[ndon+2*i+j];
if(a>b) upp=a; else{
for(k=0;k<nstep;k++){
pval2[i][j]=(a+b)/2.;
for(l=0;l<ndon;l++)
for(ll=0;ll<2;ll++)
pstart2[l][ll]=maxvals[ndon+2*l+ll];
for(l=0;l<n2loc;l++) if(pind2[l][j]==i &&
pstart2[pind2[l][1-j]][1-j]+pval2[i][j]>1.){
if(j==0) pstart2[pind2[l][1-j]][1-j]=.2*(1.-pval2[i][j]);
if(j==1) pstart2[pind2[l][1-j]][1-j]=.95*(1.-pval2[i][j]);
}
for(l=0;l<n1loc;l++) if(pind1[l][j]==i &&
pstart2[pind1[l][1-j]][1-j]+pval2[i][j]>1.){
if(j==0) pstart2[pind1[l][1-j]][1-j]=.2*(1.-pval2[i][j]);
if(j==1) pstart2[pind1[l][1-j]][1-j]=.95*(1.-pval2[i][j]);
}
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype2,
dtmp);
/*for(l=0;l<ndon;l++)
printf("%7e ",s[l]);
printf("\n");
for(m=0;m<ndon;m++){
for(l=0;l<2;l++)
printf("%7e ",p[m][l]);
printf("\n");
}
for(l=0;l<nal;l++)
printf("%7e ",am[l]);
printf("\n");
for(l=0;l<nal;l++)
printf("%7e ",co[l]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>
maxvals[3*ndon+2*nal]- *quan/2.)
a=(a+b)/2.;
else b=(a+b)/2.;
/*printf("%7e %7e\n",a,b);*/
fflush(stdout);
}
upp=b;
}
/*printf("%7e\n",upp);*/
a=min;
b=maxvals[ndon+2*i+j];
if(a>b) low=b; else{
if(j==1) nstep=nstep3;
for(k=0;k<nstep;k++){
pval2[i][j]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype2,
dtmp);
/*for(l=0;l<ndon;l++)
printf("%7e ",s[l]);
printf("\n");
for(m=0;m<ndon;m++){
for(l=0;l<2;l++)
printf("%7e ",p[m][l]);
printf("\n");
}
for(l=0;l<nal;l++)
printf("%7e ",am[l]);
printf("\n");
for(l=0;l<nal;l++)
printf("%7e ",co[l]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>
maxvals[3*ndon+2*nal]- *quan/2.)
b=(a+b)/2.;
else a=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
low=a;
}
/*printf("%7e\n",low);*/
printf("d%d %3d %15f %9f %7f\n",j+1,i,
maxvals[ndon+2*i+j],low,upp);
fflush(stdout);
}
}
}
}
for(k=0;k<ndon;k++) for(l=0;l<2;l++){
pfix2[k][l]=pfix[k][l];
pval2[k][l]=pval[k][l];
}
printf("\nAmplification parameters\n\n");
printf("Parameter Estimate Interval\n");
/*for(i=0;i<nal;i++){*/
for(i=0;i<0;i++){
for(j=0;j<nal;j++){
amfix2[j]=amfix[j];
amval2[j]=amval[j];
}
nt=1;
for(j=0;j<n2loc;j++) for(l=0;l<4;l++) if(amal2[j][l]==i &&
amfix[i]==0) nt=0;
for(j=0;j<n1loc;j++) for(l=0;l<2;l++) if(amal1[j][l]==i &&
amfix[i]==0) nt=0;
if(nt==0){
amfix2[i]=1;
a=maxvals[3*ndon+i];
b=max;
if(a>b) upp=a; else{
for(j=0;j<nstep1;j++){
amval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) a=(a+b)/2.;
else b=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
upp=b;
}
/*printf("%7e\n",upp);*/
a=min;
b=maxvals[3*ndon+i];
if(a>b) low=b; else{
for(j=0;j<nstep1;j++){
amval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike); */
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) b=(a+b)/2.;
else a=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
low=a;
}
/*printf("%7e\n",low);*/
printf("a %4d %15f %9f %7f\n",i,maxvals[3*ndon+i],low,upp);
fflush(stdout);
}
}
for(j=0;j<nal;j++){
amfix2[j]=amfix[j];
amval2[j]=amval[j];
}
printf("\nContamination parameters\n\n");
printf("Parameter Estimate Interval\n");
/*for(i=0;i<nal;i++){*/
for(i=0;i<0;i++){
for(j=0;j<nal;j++){
cofix2[j]=cofix[j];
coval2[j]=coval[j];
}
nt=1;
for(j=0;j<n2loc;j++) for(l=0;l<4;l++) if(coal2[j][l]==i &&
cofix[i]==0) nt=0;
for(j=0;j<n1loc;j++) for(l=0;l<2;l++) if(coal1[j][l]==i &&
cofix[i]==0) nt=0;
if(nt==0){
cofix2[i]=1;
a=maxvals[3*ndon+nal+i];
b=max;
if(a>b) upp=a; else{
for(j=0;j<nstep2;j++){
coval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) a=(a+b)/2.;
else b=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
upp=b;
}
/*printf("%7e\n",upp);*/
a=min;
b=maxvals[3*ndon+nal+i];
if(a>b) low=b; else{
for(j=0;j<nstep3;j++){
coval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) b=(a+b)/2.;
else a=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
low=a;
}
/*printf("%7e\n",low);*/
printf("c %4d %15f %9f %7f\n",i,maxvals[3*ndon+nal+i],low,upp);
fflush(stdout);
}
}
}
#include <math.h>
#include <stdio.h>
/* maxvals[3*nd(n+2*nal+1] is the vector of parameter values that
maximize the
likelihood, with last entry the maximized likelihood, quan is
qchisq(.95,1) for a 95% confidence interval, qchisq(.68,1) to get
standard errors */
void spermanalysisfssdioemcisonly(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *maxvals,double *quan,
int *sfix2,int **pfix2,int *amfix2,int *cofix2,double *sval2,
double **pval2,double *amval2,double *coval2,double *sstart,double
**pstart,double *amstart,double *costart,double *dconv,int *niter,int
nstep1,double *sstart2,double **pstart2,double *amstart2,
double *costart2,int **ptype,double *dtmp)
{
int i,j,nt;
double eps,min,max,a,b,upp,low;
void spermanalysisrfssdem(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *sstart,double **pstart,double
*amstart,double *costart,double *dconv,int *niter,int simind,int
**ptype,double *dtmp);
eps=pow(10.,-10.);
min=eps;
max=1-eps;
for(i=0;i<ndon;i++){
sfix2[i]=sfix[i];
sval2[i]=sval[i];
pfix2[i][0]=pfix[i][0];
pval2[i][0]=pval[i][0];
pfix2[i][1]=pfix[i][1];
pval2[i][1]=pval[i][1];
}
for(i=0;i<nal;i++){
amfix2[i]=amfix[i];
amval2[i]=amval[i];
cofix2[i]=cofix[i];
coval2[i]=coval[i];
}
for(i=0;i<ndon;i++){
sstart2[i]=maxvals[i];
for(j=0;j<2;j++) pstart2[i][j]=maxvals[ndon+2*i+j];
}
for(i=0;i<nal;i++){
amstart2[i]=maxvals[3*ndon+i];
costart2[i]=maxvals[3*ndon+nal+i];
}
printf("Confidence intervals for segregation parameters\n\n");
printf("Parameter Estimate Interval\n");
for(i=0;i<ndon;i++){
for(j=0;j<ndon;j++){
sfix2[j]=sfix[j];
sval2[j]=sval[j];
}
nt=1;
for(j=0;j<n2loc;j++) if(sind2[j]==i && sfix[i]==0) nt=0;
for(j=0;j<n1loc;j++) if(sind1[j]==i && sfix[i]==0) nt=0;
if(nt==0){
sfix2[i]=1;
a=maxvals[i];
b=max;
if(a>b) upp=a; else{
for(j=0;j<nstep1;j++){
sval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[k]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) a=(a+b)/2.;
else b=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
upp=b;
}
/*printf("%7e\n",upp);*/
a=min;
b=maxvals[i];
if(a>b) low=b; else{
for(j=0;j<nstep1;j++){
sval2[i]=(a+b)/2.;
spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);
/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
printf("\n");
for(k=0;k<ndon;k++){
for(l=0;l<2;l++)
printf("%7e ",p[k][l]);
printf("\n");
}
for(k=0;k<nal;k++) printf("%7e ",am[k]);
printf("\n");
for(k=0;k<nal;k++) printf("%7e ",co[k]);
printf("\n");
printf("%7e\n",*maxloglike);*/
if(*maxloglike>maxvals[3*ndon+2*nal]-
*quan/2.) b=(a+b)/2.;
else a=(a+b)/2.;
/*printf("%7e %7e\n",a,b);
fflush(stdout);*/
}
low=a;
}
/*printf("%7e\n",low);*/
printf("s %4d %15f %9f %7f\n",i,maxvals[i],low,upp);
fflush(stdout);
}
}
}
/* This routine calculates the chi-square quantile (1 df) associated
with the probability cilevel, i.e. this is the inverse cdf of a
chi-square on 1 df. */
#include <math.h>
#include <stdio.h>
void chisqquan(double cilevel,double *quan,double *tmp){
int i,j;
double gammacdf(double df, double z,double *tmp);
for(i=1;gammacdf(.5,(double)i/2.,tmp+2)<cilevel;i++);
*tmp=(double)(i-1);
tmp[1]=(double)i;
for(j=0;j<24;j++){
*quan=(*tmp+tmp[1])/2.;
if(gammacdf(.5,*quan/2.,tmp+2)<cilevel) *tmp=*quan;
else tmp[1]=*quan;
}
*quan=tmp[1];
}
/* This routine calculates the log of the gamma function */
double lngamma(double z,double *tmp){
*tmp=(z+.5)*log(z+5.5)-z-5.5;
*tmp+=log(2.5066282746310005*(1.000000000190015+76.18009172947146/
(z+1)-86.50532032941677/(z+2)+24.01409824083091/(z+3)-
1.231739572450155/(z+4)+0.1208650973866179e-2/(z+5)-
0.5395239384953e-5/(z+6))/z);
return *tmp;
}
/* This routine calculates the series representation of the incomplete
gamma function */
void igseries(double *ans,double df,double z,double *tmp){
double lngamma(double z,double *tmp);
int i;
if(z<=0.) *ans=0.;
else{
tmp[2]=1./df;
tmp[1]=df;
*tmp=1./df;
for(i=0;i<200;i++){
tmp[1]++;
*tmp*=(z/tmp[1]);
tmp[2]+=*tmp;
}
*ans=exp(-z+df*log(z)-lngamma(df,tmp))*tmp[2];
}
}
/* This routine calculates the continued fraction representation of the
incomplete gamma function */
void igcontfrac(double *ans,double df,double z,double *tmp){
double lngamma(double z,double *tmp);
int i;
tmp[6]=pow(10.,-30.);
tmp[1]=z+1.-df;
tmp[2]=1./tmp[6];
tmp[3]=1./tmp[1];
tmp[5]=tmp[3];
for(i=1;i<=200;i++){
*tmp=-(double)i*((double)i-df);
tmp[1]+=2.;
tmp[3]=*tmp*tmp[3]+tmp[1];
if(fabs(tmp[3])<tmp[6]) tmp[3]=tmp[6];
tmp[2]=tmp[1]+*tmp/tmp[2];
if(fabs(tmp[2])< tmp[6]) tmp[2]=tmp[6];
tmp[3]=1./tmp[3];
tmp[4]=tmp[3]*tmp[2];
tmp[5]*=tmp[4];
if(fabs(tmp[4]-1.)<pow(10.,-8.)) break;
}
*ans=exp(-z+df*log(z)-lngamma(df,tmp))*tmp[5];
}
/* This routines calculates the cdf of the gamma density */
double gammacdf(double df,double z,double *tmp){
void igcontfrac(double *ans,double df,double z,double *tmp);
void igseries(double *ans,double df,double z,double *tmp);
if(z<(df+1.)){
igseries(tmp+7,df,z,tmp);
return tmp[7];
} else{
igcontfrac(tmp+7,df,z,tmp);
return 1.-tmp[7];
}
}
/*Must use with seed.c*/
#include <math.h>
#include <stdio.h>
double randomnum(long *seed,long *ab,long *cd)
{
int j;
long k;
double temp;
if (*seed <= 0 || !*cd){
if (-(*seed) < 1) *seed=1;
else *seed = -(*seed);
for (j=32+7; j >=0; j--){
k = (*seed)/127773;
*seed = 16807*(*seed-k*127773)-2836*k;
if (*seed < 0) *seed += 2147483647;
if (j<32) ab[j] = *seed;
}
*cd = ab[0];
}
k = (*seed)/127773;
*seed = 16807*(*seed-k*127773)-2836*k;
if (*seed <0) *seed += 2147483647;
j = *cd/(1+1073741823/16);
*cd = ab[j];
ab[j] = *seed;
if ((temp = (1./2147483647)**cd) > 1.-1.2e-7) temp = (1.-1.2e-7);
return temp;
}