A new method is proposed to solve systems of linear approximate equations X θ ≈ y where the unknowns θ and the data y are positive and the matrix X consists of nonnegative elements. Writing the i-th near-equality Xi.θ/yi ≈ 1 the assumed model is Xi.θ/yi = ζi with mutually independent positive errors ζi. The loss function is defined by ∑wi(ζi − 1) logζi in which wi is the importance weight for the i-th near-equality. A reparameterization reduces the method to unconstrained minimization of a smooth strictly convex function implying the unique existence of positive solution and the applicability of Newton's method that converges quadratically. The solution stability is controlled by weighting prior guesses of the unknowns θ. The method matches the maximum likelihood estimation if all weights wi are equal and ζi independently follow the probability density function ∝ tw(1-t), 0 < w.