function model = GMRtrain(trY, trX, opt);

    % trY = AU values
    %trX = pose values

    % size(trY) % = 1217, 1
    % size(trX) % = 1217, 3; The number of columns may vary in the future
    
	%generate GaussianMixture  pdfs
	Xm = min(trX); % size = 1, 3; gets the minimum for each column (in this case each pose dimension: yaw, pitch, roll)
	Xr = max(trX) - Xm; % size = 1, 3; range of each column
	E =ones(opt.nGM,1); % size = 8, 1
    
    
    % size(Xm)
    % size(Xr)
    % size(E)
    % (1:opt.nGM) = [ 1 2 3 4 5 6 7 8 ... nGM ]
	Gmean = (1:opt.nGM)'./opt.nGM * Xr + E*Xm;
	Gsigma = (Xr./opt.nGM);
	ntrX = size(trX,1);

	%project trX ont GM
	trXp = zeros(ntrX, length(Gmean(:)));
    tmpSigma1 = (2*E*Gsigma.^2);
    tmpSigma2 = (E*Gsigma * sqrt(2*pi));

%    size(trX)
%    size(Gmean)
	for i=1:ntrX
		f=exp(-(E*trX(i,:) - Gmean).^2 ./ (2*E*Gsigma.^2)) ./ (E*Gsigma * sqrt(2*pi));
%		f=exp(-(E*trX(i,:) - Gmean).^2 ./ tmpSigma1) ./ tmpSigma2;
%		f=E*trX(i,:) - Gmean;
%		f=-(E*trX(i,:) - Gmean).^2;
%		f=exp(f ./ tmpSigma1);
%		f = f ./ tmpSigma2;
		trXp(i,:) = f(:)';
	end

	%linear regression
	c = trXp \ trY; % Same as inv(trXp' * trXp) * trXp' * trY
	c1 = inv(trXp' * trXp) * trXp' * trY;
	% c1 = trXp' * trXp;
	% c1 = inv(trXp' * trXp);
	% c1 = inv(trXp' * trXp) * trXp';
	% Same as inv(trXp' * trXp) * trXp' * trY

    % size(trXp) % = 1217, 24
    % size(trY) % = 1217, 1
    % size(Gmean) = 8, 3
    % size(Gsigma) = 1, 3
    % size(c) = 24, 1

	%make model
	model.Gmean = Gmean;
	model.Gsigma = Gsigma;
	model.c = c;
	model.c1 = c1;
	model.trXp = trXp;
