Function runEM

src/myFunc.cpp:1463–1534

}

float runEM (const vector<float>& x,const vector<float>& y,double & a0,double & a1, double & a2, double & a3, int maximalNumberOfIterations,int ploidy, int maximalNumberOfCopies) {

	float rmserror = -1;

	if (x.size() != y.size()) {
		cerr << "Error: object size is different";
		return 0;
	}
	vector <int> cluster (x.size());
	vector <float> res (maximalNumberOfCopies+1);

	int count = 0;

	for (count = 0; count < maximalNumberOfIterations; count++) {

		for (int i = 0; i <(int)x.size(); i++) {
			for (int j = 0; j <= maximalNumberOfCopies; j++)
				res[j] = fabs(polynomial(x[i],j*a0/ploidy,j*a1/ploidy,j*a2/ploidy,j*a3/ploidy)-y[i]);
			cluster[i] = get_min_index(res);
		}
		int npoints = 0;
		for (int i = 0; i <(int)x.size(); i++)
			if (cluster[i] == ploidy)
				npoints++;
		int nvars = 3; //fit by cubic polynomial
		ap::real_2d_array xy;
		xy.setlength(npoints,nvars+1);
		int pos = 0;
		for (int i = 0; i <(int)x.size(); i++) {
			if (cluster[i] == ploidy){
				xy(pos,2) = x[i];
				xy(pos,0) = x[i]*x[i]*x[i];
				xy(pos,1) = x[i]*x[i];
				xy(pos,3) = y[i];
				pos++;
			}
		}

		linearmodel lm;
		int info;
		lrreport ar;
		lrbuild(xy,npoints,nvars,info,lm,ar);
		if (info != 1) {
			cerr << "Error in linear regression, code: " << info <<"\n";
			break;
		}
		rmserror = float(ar.rmserror);
		ap::real_1d_array v;
		v.setlength(nvars+1);
		lrunpack(lm,v,nvars);
		//cout << v(0);  //x
		//cout << v(1);  //x^3
		//cout << v(2);  //x^2
		//cout << v(3);  //intercept

		cout << a0 - v(0) << "\t"; //x^3
		cout << a1 - v(1) << "\t"; //x^2
		cout << a2 - v(2) << "\t"; //x