@@ -14,7 +14,7 @@ void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* m

 // 	FILE* pFile = fopen("chromStar.log", "w");

 // 	Output2FILE::Stream() = pFile;

  	FILELog::ReportingLevel() = FILELog::FromString("ERROR");

-//  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG1");

+//  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG2");

 	// Parallelization settings

 	omp_set_num_threads(*num_threads);

@@ -337,8 +337,8 @@ void NegativeBinomial::update_constrained(double** weights, int fromState, int t

 	{

 		for (int t=0; t<this->T; t++)

 		{

-			numerator+=weights[i+fromState][t]*this->size*(i+1);

-			denominator+=weights[i+fromState][t]*(this->size*(i+1)+this->obs[t]);

+			numerator += weights[i+fromState][t] * this->size*(i+1);

+			denominator += weights[i+fromState][t] * (this->size*(i+1) + this->obs[t]);

 		}

 	}

 	this->prob = numerator/denominator; // Update of size (r) is now done with old prob

@@ -371,22 +371,22 @@ void NegativeBinomial::update_constrained(double** weights, int fromState, int t

 				{

 					if(this->obs[t]==0)

 					{

-						Fr+=weights[i+fromState][t]*(i+1)*logp;

+						Fr += weights[i+fromState][t] * (i+1) * logp;

 						//dFrdr+=0;

 					}

 					if(this->obs[t]!=0)

 					{

-						Fr+=weights[i+fromState][t]*(i+1)*(logp-DigammaR+DigammaRplusX[(int)obs[t]]);

-						dFrdr+=weights[i+fromState][t]/dr*(i+1)*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX[(int)obs[t]]-DigammaRplusX[(int)obs[t]]);

+						Fr += weights[i+fromState][t] * (i+1) * (logp - DigammaR + DigammaRplusX[(int)obs[t]]);

+						dFrdr += weights[i+fromState][t] / dr * (i+1) * (DigammaR - DigammaRplusDR + DigammaRplusDRplusX[(int)obs[t]] - DigammaRplusX[(int)obs[t]]);

 					}

 				}

 				if(fabs(Fr)<eps)

 				{

 					break;

 				}

-				if(Fr/dFrdr<rhere) rhere=rhere-Fr/dFrdr;

-				if(Fr/dFrdr>rhere) rhere=rhere/2.0;

 			}

+			if(Fr/dFrdr<rhere) rhere=rhere-Fr/dFrdr;

+			if(Fr/dFrdr>rhere) rhere=rhere/2.0;

 		}

 	}

 	else

@@ -406,13 +406,13 @@ void NegativeBinomial::update_constrained(double** weights, int fromState, int t

 					DigammaRplusDRplusX = digamma((i+1)*(rhere+dr)+this->obs[t]); // boost::math::digamma<>(rhere+dr+this->obs[ti]);

 					if(this->obs[t]==0)

 					{

-						Fr+=weights[i+fromState][t]*(i+1)*logp;

+						Fr += weights[i+fromState][t] * (i+1) * logp;

 						//dFrdr+=0;

 					}

 					if(this->obs[t]!=0)

 					{

-						Fr+=weights[i+fromState][t]*(i+1)*(logp-DigammaR+DigammaRplusX);

-						dFrdr+=weights[i+fromState][t]/dr*(i+1)*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX-DigammaRplusX);

+						Fr += weights[i+fromState][t] * (i+1) * (logp - DigammaR + DigammaRplusX);

+						dFrdr += weights[i+fromState][t] / dr * (i+1) * (DigammaR - DigammaRplusDR + DigammaRplusDRplusX - DigammaRplusX);

 					}

 				}

 			}

@@ -262,34 +262,7 @@ void ScaleHMM::baumWelch(int* maxiter, int* maxtime, double* eps)

 // 			this->densityFunctions[iN]->update(this->gamma[iN]);

 // 		}

-// 		// Update distribution of state 'null-mixed' and 'monosomic', set others as multiples of 'monosomic'

-// 		// This loop assumes that the negative binomial states come last and are consecutive

-// 		for (int iN=0; iN<this->N; iN++)

-// 		{

-// 			if (this->densityFunctions[iN]->get_name() == ZERO_INFLATION) {}

-// 			if (this->densityFunctions[iN]->get_name() == GEOMETRIC)

-// 			{

-// 				this->densityFunctions[iN]->update(this->gamma[iN]);

-// 			}

-// 			if (this->densityFunctions[iN]->get_name() == NEGATIVE_BINOMIAL)

-// 			{

-// 				FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

-// 				this->densityFunctions[iN]->update_constrained(this->gamma, iN, this->N);

-// 				double mean1 = this->densityFunctions[iN]->get_mean();

-// 				double variance1 = this->densityFunctions[iN]->get_variance();

-// 				FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

-// 				// Set others as multiples

-// 				for (int jN=iN+1; jN<this->N; jN++)

-// 				{

-// 					this->densityFunctions[jN]->set_mean(mean1 * (jN-iN+1));

-// 					this->densityFunctions[jN]->set_variance(variance1 * (jN-iN+1));

-// 					FILE_LOG(logDEBUG1) << "mean(state="<<jN<<") = " << this->densityFunctions[jN]->get_mean() << ", var(state="<<jN<<") = " << this->densityFunctions[jN]->get_variance();

-// 				}

-// 				break;

-// 			}

-// 		}

-			

-		// Update distribution of state 'null-mixed' and 'disomic', set others as multiples of 'disomic'/2

+		// Update distribution of state 'null-mixed' and 'monosomic', set others as multiples of 'monosomic'

 		// This loop assumes that the negative binomial states come last and are consecutive

 		int xsomy = 1;

 		for (int iN=0; iN<this->N; iN++)

@@ -301,18 +274,18 @@ void ScaleHMM::baumWelch(int* maxiter, int* maxtime, double* eps)

 			}

 			if (this->densityFunctions[iN]->get_name() == NEGATIVE_BINOMIAL)

 			{

-				if (xsomy==2)

+				if (xsomy==1)

 				{

 					FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

-					this->densityFunctions[iN]->update(this->gamma[iN]);

-					double mean1 = this->densityFunctions[iN]->get_mean() / xsomy;

-					double variance1 = this->densityFunctions[iN]->get_variance() / xsomy;

+					this->densityFunctions[iN]->update_constrained(this->gamma, iN, this->N);

+					double mean1 = this->densityFunctions[iN]->get_mean();

+					double variance1 = this->densityFunctions[iN]->get_variance();

 					FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

 					// Set others as multiples

-					for (int jN=2; jN<this->N; jN++)

+					for (int jN=iN+1; jN<this->N; jN++)

 					{

-						this->densityFunctions[jN]->set_mean(mean1 * (jN-iN+2));

-						this->densityFunctions[jN]->set_variance(variance1 * (jN-iN+2));

+						this->densityFunctions[jN]->set_mean(mean1 * (jN-iN+1));

+						this->densityFunctions[jN]->set_variance(variance1 * (jN-iN+1));

 						FILE_LOG(logDEBUG1) << "mean(state="<<jN<<") = " << this->densityFunctions[jN]->get_mean() << ", var(state="<<jN<<") = " << this->densityFunctions[jN]->get_variance();

 					}

 					break;

@@ -321,6 +294,38 @@ void ScaleHMM::baumWelch(int* maxiter, int* maxtime, double* eps)

 			}

 		}

+// 		// Update distribution of state 'null-mixed' and 'disomic', set others as multiples of 'disomic'/2

+// 		// This loop assumes that the negative binomial states come last and are consecutive

+// 		int xsomy = 1;

+// 		for (int iN=0; iN<this->N; iN++)

+// 		{

+// 			if (this->densityFunctions[iN]->get_name() == ZERO_INFLATION) {}

+// 			if (this->densityFunctions[iN]->get_name() == GEOMETRIC)

+// 			{

+// 				this->densityFunctions[iN]->update(this->gamma[iN]);

+// 			}

+// 			if (this->densityFunctions[iN]->get_name() == NEGATIVE_BINOMIAL)

+// 			{

+// 				if (xsomy==2)

+// 				{

+// 					FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

+// 					this->densityFunctions[iN]->update(this->gamma[iN]);

+// 					double mean1 = this->densityFunctions[iN]->get_mean() / xsomy;

+// 					double variance1 = this->densityFunctions[iN]->get_variance() / xsomy;

+// 					FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

+// 					// Set others as multiples

+// 					for (int jN=2; jN<this->N; jN++)

+// 					{

+// 						this->densityFunctions[jN]->set_mean(mean1 * (jN-iN+2));

+// 						this->densityFunctions[jN]->set_variance(variance1 * (jN-iN+2));

+// 						FILE_LOG(logDEBUG1) << "mean(state="<<jN<<") = " << this->densityFunctions[jN]->get_mean() << ", var(state="<<jN<<") = " << this->densityFunctions[jN]->get_variance();

+// 					}

+// 					break;

+// 				}

+// 				xsomy++;

+// 			}

+// 		}

+			

 		dtime = clock() - clocktime;

 	 	FILE_LOG(logDEBUG) << "updating distributions: " << dtime << " clicks";

 		R_CheckUserInterrupt();