bin/IDEAS_2018/IDEAS_APP/genomicTensor.cpp
//mle of mean and precision parameters
mle of mean and precision parameters
mle of mean and precision parameters
void genomicTensor::_getMLE_HyperP1(double *logp, double *ak, double &M, int sz, double N, double precision, int maxitern)
{	int i, cn = 0;
	double newak[sz], newM;
	if(N<=0) return;
	while(cn < maxitern)
	{	double s, err, err2, y = 0;
		for(i = 0; i < sz; i++)
		{	
			y += ak[i] * (logp[i] - gsl_sf_psi(M * ak[i]));
		}
		s = 0;
		for(i = 0; i < sz; i++)
		{	
			newak[i] = _solveDigamma(ak[i], logp[i] - y, precision, maxitern);
			s += newak[i];
		}
		err = y = 0;
		for(i = 0; i < sz; i++) 
		{	newak[i] /= s;
			err += (ak[i] - newak[i]) * (ak[i] - newak[i]) / ak[i];
			ak[i] = newak[i];
			y += ak[i] * (logp[i] - gsl_sf_psi(M * ak[i]));
		}
		y = 0;
		for(i = 0; i < sz; i++)
			y += ak[i] * ak[i] * gsl_sf_psi_1(M * ak[i]);	
		newM = M;
		err2 = (M - newM) * (M - newM) / M;
		M = newM;
		if(err < precision && err2 < precision) break;
		cn++;
	}
}
void genomicTensor :: _getMLE_HyperP1 ( double * logp , double * ak , double & M , int sz , double N , double precision , int maxitern ) { int i , cn = 0 ; double newak [ sz ] , newM ; if ( N <= 0 ) return ; while ( cn < maxitern ) { double s , err , err2 , y = 0 ; for ( i = 0 ; i < sz ; i ++ ) { y += ak [ i ] * ( logp [ i ] - gsl_sf_psi ( M * ak [ i ] ) ) ; } s = 0 ; for ( i = 0 ; i < sz ; i ++ ) { newak [ i ] = _solveDigamma ( ak [ i ] , logp [ i ] - y , precision , maxitern ) ; s += newak [ i ] ; } err = y = 0 ; for ( i = 0 ; i < sz ; i ++ ) { newak [ i ] /= s ; err += ( ak [ i ] - newak [ i ] ) * ( ak [ i ] - newak [ i ] ) / ak [ i ] ; ak [ i ] = newak [ i ] ; y += ak [ i ] * ( logp [ i ] - gsl_sf_psi ( M * ak [ i ] ) ) ; } y = 0 ; for ( i = 0 ; i < sz ; i ++ ) y += ak [ i ] * ak [ i ] * gsl_sf_psi_1 ( M * ak [ i ] ) ; newM = M ; err2 = ( M - newM ) * ( M - newM ) / M ; M = newM ; if ( err < precision && err2 < precision ) break ; cn ++ ; } }
mle of mean and precision parameters
mle of mean and precision parameters
//min(100.,max(1.,M - ds / ds2));