randomlib.c

Fit a multivariate gaussian mixture by a cross-entropy method. Cross-Entropy is a powerfull tool to

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include "randomlib.h"
#include "matrixjpl.h"



#pragma warning(disable: 4244)
#pragma warning(disable: 4305)

#define ABS(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))



void normal_truncc(int n, double *mu, double *sigma2, double *left, double *right, double *out)
{
	int i;
    double *lowerProb, *upperProb, *u, *sqrtsig;
	double *phi1, *phi2;
    
u = dvector(0,n-1);
sqrtsig = dvector(0,n-1);
lowerProb = dvector(0,n-1);
upperProb = dvector(0,n-1);
phi1 = dvector(0,n-1);
phi2 = dvector(0,n-1);

for(i=0; i<n; i++)
sqrtsig[i] = sqrt(sigma2[i]);

  for(i=0; i<n; i++){   
  lowerProb[i] = (left[i]  - mu[i])/sqrtsig[i];
  upperProb[i] = (right[i] - mu[i])/sqrtsig[i];
  }
  
  for(i=0; i<n; i++){
  phi1[i] = lowerProb[i]/sqrt(2.0);
  phi2[i] = upperProb[i]/sqrt(2.0);
  lowerProb[i] = 0.5*(1.0+erf1(&phi1[i]));
  upperProb[i] = 0.5*(1.0+erf1(&phi2[i]));
  }

  for(i=0; i<n; i++)
      u[i] = lowerProb[i] + (upperProb[i] - lowerProb[i])*ranf();

 for(i=0; i<n; i++)
      out[i] = mu[i] + phiinv(u[i])*sqrtsig[i];

  free_dvector(u,0);
  free_dvector(lowerProb,0);
  free_dvector(upperProb,0);
  free_dvector(sqrtsig,0);
  free_dvector(phi1,0);
  free_dvector(phi2,0);

}


void snormal_rndc(int n, double *randvec)
{ 
	int i;
     
	 for(i=0; i<n; i++)
	 randvec[i] = snorm();

}

void normal_rndc(double **vcov, int k, double *randvec)
{ 
	int i, j;
	double *rvec, *xpxd;
	int invt;

	rvec  = dvector(0,k-1);
	xpxd  = dvector(0,k-1);

    invt = choldcmp(vcov, k, xpxd);
     if (invt != 1)
		 mexErrMsgTxt("cholesky error in normal_rnd \n");

	  
	 for(i=0; i<k; i++){
	 for(j=i; j<k; j++)
		 vcov[i][j] = 0.0;
		 vcov[i][i] = xpxd[i];
	}

	 for(i=0; i<k; i++)
	 rvec[i] = snorm();
	 
	 matvec(vcov,k,k,rvec,randvec);

	 free_dvector(rvec,0);
	 free_dvector(xpxd,0);

}

void beta_rndc(int n, double a, double b, double *array)
{
    int i;

    /*
     BETA deviates
    */
    for(i=0; i<n; i++)   
    *(array+i) = genbet(a,b);
    
}

void beta_cdfc(int n, double a, double b, double *x, double *array)
{
    int i;
	double p, q;
	double *y;
	
	y = dvector(0,n-1);

	for(i=0; i<n; i++){
		y[i] = 1.0 - x[i];
	    cumbet(&x[i],&y[i],&a,&b,&p,&q);
	    *(array+i) = p;
		}
    free_dvector(y,0);
}

void bino_rndc(int n, long a, double p, double *array)
{
	int i;

	/*
     Binomial deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = ignbin(a,p);
    
}

void bino_cdfc(int n, double a, double prob, double *x, double *array)
{
    int i;
	double p, q;
	double ompr;
	
	ompr = 1.0 - prob;
	
	for(i=0; i<n; i++){
	    cumbin(&x[i],&a,&prob,&ompr,&p,&q);
	    *(array+i) = p;
		}
}

void chis_rndc(int n, double dof, double *array)
{
	int i;

	/*
     Chi-square deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = genchi(dof);
    
}

void chis_cdfc(int n, double a, double *x, double *array)
{
    int i;
	double p, q;
	
	for(i=0; i<n; i++){
	    cumchi(&x[i],&a,&p,&q);
	    *(array+i) = p;
		}
}

void exp_rndc(int n, double dof, double *array)
{
	int i;

	/*
     exponential deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = genexp(dof);
    
}

void fdis_rndc(int n, double a, double b, double *array)
{
	int i;

	/*
     F- deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = genf(a,b);   
}

void fdis_cdfc(int n, double a, double b, double *x, double *array)
{
    int i;
	double p, q;

	for(i=0; i<n; i++){
	    cumf(&x[i],&a,&b,&p,&q);
	    *(array+i) = p;
		}
}

void gamm_rndc(int n, double a, double b, double *array)
{
	int i;

	/*
     Gamma deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = gengam(a,b);
    
}

void gamm_cdfc(int n, double a, double *x, double *array)
{
    int i;
	double p, q;
	
	for(i=0; i<n; i++){
	    cumgam(&x[i],&a,&p,&q);
	    *(array+i) = p;
		}
}

void mdis_rndc(int n, int nevents, long ncat, double *p, double **array)
{
	int i, j;
	unsigned long *a;

	a = ivector(0,ncat);
	/*
     multinomial deviates
    */
    for(i=0; i<n; i++) {
     genmul(nevents,p,ncat,a);
     mexPrintf(" a = %ld \n",a[0]);
	 for (j=0; j<ncat; j++)
     array[i][j] = (double) a[j];
	}

	free_ivector(a,0);
    
}

void nbino_rndc(int n, long a, double p, double *array)
{
	int i;

	/*
     Negative Binomial deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = ignnbn(a,p);
    
}

void nbino_cdfc(int n, double a, double prob, double *x, double *array)
{
    int i;
	double p, q;
	double ompr;
	
	ompr = 1.0 - prob;
	
	for(i=0; i<n; i++){
	    cumnbn(&x[i],&a,&prob,&ompr,&p,&q);
	    *(array+i) = p;
		}
}

void nchi_rndc(int n, double a, double b, double *array)
{
	int i;

	/*
     non-central chi-squared deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = gennch(a,b);
    
}

void nchi_cdfc(int n, double a, double b, double *x, double *array)
{
    int i;
	double p, q;

	for(i=0; i<n; i++){
	    cumchn(&x[i],&a,&b,&p,&q);
	    *(array+i) = p;
		}
}


void nfdis_rndc(int n, double a, double b, double c, double *array)
{
	int i;

	/*
     non-central F deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = gennf(a,b,c);
    
}

void nfdis_cdfc(int n, double a, double b, double c, double *x, double *array)
{
    int i;
	double p, q;

	for(i=0; i<n; i++){
	    cumfnc(&x[i],&a,&b,&c,&p,&q);
	    *(array+i) = p;
		}
}

void pois_rndc(int n, double mu, double *array)
{
	int i;
	/*
     poisson deviates
    */
    for(i=0; i<n; i++) 
    *(array+i) = ignpoi(mu);
}

void pois_cdfc(int n, double a, double *x, double *array)
{
    int i;
	double p, q;
	
	for(i=0; i<n; i++){
	    cumpoi(&x[i],&a,&p,&q);
	    *(array+i) = p;
		}
}

void tdis_rndc(int n, double dof, double *array)
{
	int i;
	double *ndis, *cdis;

	ndis = dvector(0,n-1);
	cdis = dvector(0,n-1);

	snormal_rndc(n, ndis);
	chis_rndc(n, dof, cdis);

	/*
     t-distribution random deviates from chi-squared
    */
    for(i=0; i<n; i++) 
    *(array+i) = ndis[i]*sqrt(dof)/sqrt(cdis[i]);

	free_dvector(ndis,0);
	free_dvector(cdis,0);
}

void tdis_cdfc(int n, double a, double *x, double *array)
{
    int i;
	double p, q;
	
	for(i=0; i<n; i++){
	    cumt(&x[i],&a,&p,&q);
	    *(array+i) = p;
		}
}




void ftnstop(char*);
float genbet(float aa,float bb)
/* ********************************************************************
     float genbet(float aa,float bb)
               GeNerate BETa random deviate
                              Function
     Returns a single random deviate from the beta distribution with
     parameters A and B.  The density of the beta is
               x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
                              Arguments
     aa --> First parameter of the beta distribution
       
     bb --> Second parameter of the beta distribution
       
                              Method
     R. C. H. Cheng
     Generating Beta Variatew with Nonintegral Shape Parameters
     Communications of the ACM, 21:317-322  (1978)
     (Algorithms BB and BC)
**********************************************************************
*/
{
/* JJV changed expmax (log(1.0E38)==87.49823), and added minlog */
#define expmax 87.49823
#define infnty 1.0E38
#define minlog 1.0E-37
static float olda = -1.0E37;
static float oldb = -1.0E37;
static float genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
static long qsame;

    qsame = olda == aa && oldb == bb;
    if(qsame) goto S20;
    if(!(aa < minlog || bb < minlog)) goto S10;
    mexPrintf(" AA: %16.6E BB %16.6E\n",aa,bb);
    mexErrMsgTxt("problem generating beta variate");
    exit(1);
S10:
    olda = aa;
    oldb = bb;
S20:
    if(!(min(aa,bb) > 1.0)) goto S100;
/*
     Alborithm BB
     Initialize
*/
    if(qsame) goto S30;
    a = min(aa,bb);
    b = max(aa,bb);
    alpha = a+b;
    beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
    gamma = a+1.0/beta;
S30:
S40:
    u1 = ranf();
/*
     Step 1
*/
    u2 = ranf();
    v = beta*log(u1/(1.0-u1));
/* JJV altered this */
    if(v > expmax) goto S55;
/*
 * JJV added checker to see if a*exp(v) will overflow
 * JJV S50 _was_ w = a*exp(v); also note here a > 1.0
 */
    w = exp(v);
    if(w > infnty/a) goto S55;
    w *= a;
    goto S60;
S55:
    w = infnty;
S60:
    z = pow(u1,2.0)*u2;
    r = gamma*v-1.3862944;
    s = a+r-w;
/*
     Step 2
*/