3

I have a linear system of equations in the form Ax = b. The elements of A and b were experimentally determined and as such have some uncertainty. Within each row, A and b are correlated. Between rows, A and b are uncorrelated. How would I propagate this to the elements of x?

\begin{equation} \begin{bmatrix} a_{11} \pm \sigma_{a_{11}} & a_{12} \pm \sigma_{a_{12}} & a_{13} \pm \sigma_{a_{13}} \\ a_{21} \pm \sigma_{a_{21}} & a_{22} \pm \sigma_{a_{22}} & a_{23} \pm \sigma_{a_{23}} \\ \end{bmatrix} \begin{bmatrix} x_{11} \pm \sigma_{x_{11}} \\ x_{21} \pm \sigma_{x_{21}} \\ x_{31} \pm \sigma_{x_{31}} \\ \end{bmatrix} =\begin{bmatrix} b_{11} \pm \sigma_{b_{11}} \\ b_{21} \pm \sigma_{b_{21}} \end{bmatrix} \end{equation}

, where $x_{11}+x_{21}+x_{31}=1$, and $x_{11}$, $x_{21}$, $x_{21}$ $>0$. I am trying to find $\sigma_{x_{11}}$, $\sigma_{x_{21}}$, and $\sigma_{x_{31}}$.

I have been looking at the threads on propagation of uncertainty through a linear system of equations and variance of dependent variables but I am unable to combine the two ideas.

There has been a similar unanswered question. Will be hugely grateful for any pointers!

Mark Stone has been extremely helpful with a statistical sampling method. If there is a deterministic way, or some understanding of the challenges of a deterministic way, that would also be of interest, possibly not to just myself.

Update: I tested the solution to this problem but it raised a couple of questions on its own.

Community
  • 1
Pete
  • 86
  • 10
  • Do you have 2 linear equations in 3 unknowns, and therefore, leaving aside the uncertainty, an under-determined system, not over-determined? – Mark L. Stone Mar 08 '16 at 01:05
  • Thanks Mark for your comment. Apologies, you may be right! It's very late at night here. It is indeed 2 linear equations. I have 3 unknowns, but the third unknown is a linear combination of the first two. Does it qualify as an independent third unknown? – Pete Mar 08 '16 at 01:21
  • I don't understand what constraint or relation you have on the 3rd unknown, or as some linear combination of the first two. Really need to understand what the problem is first, then can worry about the solution. – Mark L. Stone Mar 08 '16 at 01:23
  • Thank you for asking! B is an expression of the extent of x-ray CT absorption within a volume, as a function of pure components A. Here, x represents volume fractions of the components of A, and must sum to 1 (thus the constraint on the third unknown). In this problem, x represents fractions of dolomite, calcite and air; 'A' represent the experimentally measured absorptions of the three pure components; and 'B' represent the volume-averaged absorption at two different conditions where the pure components A respond differently. Please let me know if I may further clarify. – Pete Mar 08 '16 at 01:33
  • Constraint that the x's sum to 1 is your 3rd linear equation, which has no uncertainty in LHS or RHS elements. Leaving aside the uncertainty, your system is exactly determined, presuming your A instantiation (with 3rd row) is not singular.The most straightforward approach, which allows for any form of dependency, is to perform n stochastic simulation replications, where for each replication, you draw random values for A and B, and solve A_augmented x = B_augmented, where augmented means adding the row specifying x's sum to 1.Then you have n sample values of the x vector, w/ all dependencies. – Mark L. Stone Mar 08 '16 at 01:49
  • Thanks Mark, that sounds like something I could do! If the uncertainties of A and B were from samples of size m, are there any implications on n? Is there a number n that is too many or too little? I recall that the uncertainty of a random parameter decays with sqrt(n). If you copy your comment to an answer I could vote on it. – Pete Mar 08 '16 at 02:12
  • Are A and b measurements obtained in matched "sets", i.e., there are m samples each consisting of one observation of A and one of b, or are A and b measurements independent of each other? Depending on how the A and b data is collected, and how well it can be fit to a distribution, you can either use each of m A and b observation pairs as one simulation replication, or if fitting, number of replications is essentially unrelated to value of m. – Mark L. Stone Mar 08 '16 at 02:33
  • They are obtained independently, and will comprise different number of measurements. – Pete Mar 08 '16 at 02:36
  • Sorry, I had computer problems, please see the completed version of my preceding comment re: how well you can fit a distribution for errors in A and b. Is there any reason to think they are Normally distributed? Are errors in various elements of A correlated? Are errors in the 2 elements of b correlated? – Mark L. Stone Mar 08 '16 at 02:42
  • The errors within B are not correlated. Between rows of A they are not correlated. However, between A and B within the same row, they are correlated. Each row of A and B are greyscale pixels from a single 3D image (positively correlated) and their uncertainty is Gaussian noise. A and B are random and normally distributed. Distributions can be fitted well to A and B in most cases. There is an application (situation) where B may comprise a limited number of measurements and a fitted distribution may some uncertainty. – Pete Mar 08 '16 at 03:07

2 Answers2

2

Basically, the question is how to calculate uncertainties of the derived coefficients in errors-in-variables models.

One way is to use Monte Carlo simulations, as showed in Mark's answer. For more analytical solutions, you may need to read https://en.wikipedia.org/wiki/Errors-in-variables_models and other materials. I haven't used this before. Hope other people can comment on it. You may refer to this post for the difference between uncertainty from MC simulations and uncertainty from analytical derivations: Is Monte Carlo uncertainty estimation equivalent to analytical error propagation?

By the way, you said that 'Within each row, A and b are correlated'. Of course, A1 and b1 is are correlated, otherwise we won't derive X, right? By stating 'Between rows, A and b are uncorrelated', do you mean your observations of (A1, b1) and (A2, b2) are uncorrelated?

Community
  • 1
yuqian
  • 1,234
  • 9
  • 8
  • Thanks yuqian. I had a look at the wikipedia page, and I admit it was a struggle. I am primarily an experimentalist and so relegated to simpler maths. Yes (A1, b1) and (A2, b2) are uncorrelated, they come from two independent measurements. – Pete Mar 09 '16 at 00:22
1

Here is a stochastic simulation approach. I am not ruling out that someone may come up with an approach not requiring simulation. Nevertheless, the simulation run time should be negligible, so if the goal is a practical solution, this should qualify.

You have a system of 2 linear equations in 3 unknowns, $A x = b$, augmented by the additional linear equation that the sum of the x's equals 1. This adds a 3rd row to A consisting of all 1's, and a 3rd element of b equal to 1. There is no uncertainty associated with this 3rd row. Let $A_{augmented}$ be the 2 rows of A augmented by a (3rd row) row of 1's, and $b_{augmented}$ be the 2 elements of b with a 3rd element equal to 1. So your complete system of linear equations is $A_{augmented} x = b_{augmented}$

The first step is to either

1) fit distributions to the values of A and b, accounting for any correlations among the elements, including across A and b, then for each simulation replication, draw a sample of A and b from the fitted distributions (for instance, Multivariate Normal which can allow correlation between A and b) or

2) use raw measurements of A and b "as is", with each measurement (combination between A and b) being used for a simulation replication.

For each of n simulation replications, solve the equation $A_{augmented} x = b_{augmented}$, which provides the simulated value of x for that simulation replication. Once you have completed the simulation replications, you can assess the error distribution of the n sample values of x (including correlation across elements). The accuracy you achieve in assessing the uncertainty of x will ultimately be limited by the number m of A and b measurements (decreasing as its square root), and if fitting distributions to A and b, the validity of the distributions you fit. Nevertheless, if using fitted distributions for A and b, you might as well do a large number of replications n, so that you get maximum mileage out of the data you do have - no need for number of replications n to be limited to number of observations m.

Mark L. Stone
  • 12,546
  • 1
  • 31
  • 51
  • Mark, thanks for this, it is working as far as I can tell. However, I'm getting negative values of x. Since x is a physical volume fraction, how can I add a constraint for only positive values of x? – Pete Mar 08 '16 at 17:34
  • Instead of solving the system of linear equations, you could do bound-constrained linear least squares, imposing the bound constraint x >= 0. You can use linearly-constrained linear least squares, which will include bound constraints as a special case. You'll have to decide if you want to impose the x's summing to 1 as a hard constraint (which would be a linear equality constraint to the linear least squares problem), or if you want to allow some wiggle room on it. If the former, remove it from the least squares part, and impose as a linear constraint, in addition to the x >= 0 constraint. – Mark L. Stone Mar 08 '16 at 17:52
  • Thanks Mark. I simulated x and using your method, attempted to back-derive x and its error or standard deviation. I get the right mean but the standard deviation differs significantly. I simply don't understand what is wrong. If you have any comments, or pointers, I'd be immensely grateful. [The script was written in Matlab](http://pastebin.com/CjrjmbYG). This is similar to [another question](http://stats.stackexchange.com/questions/198445/is-monte-carlo-uncertainty-estimation-equivalent-to-analytical-error-propagation) posted on this site. – Pete Mar 09 '16 at 01:20
  • I don't know what you're doing on the APPROXIMATE analytical error propagation, but surely you're not accounting for the nonnegativity constraint, likely among other things. – Mark L. Stone Mar 09 '16 at 01:42
  • Thanks Mark for this. Do you mean that the constraints affect the propagation of uncertainties of A and x to b in $a_{11} x_{11} + a_{12} x_{21} + a_{13} x_{31} = b_{11}$, and $\sigma_{a11x11} = a_{11} x_{11} \sqrt{\left(\frac{\sigma_{a11}}{a_{11}}\right)^2 + \left(\frac{\sigma_{x11}}{x_{11}}\right)^2} etc$, and $\sigma_{b11}=\sqrt{\sigma_{a11x11}^2 + \sigma_{a12x21}^2 + \sigma_{a13x31}^2}$ is not valid? – Pete Mar 09 '16 at 02:09
  • Don't you want to propagate the uncertainties in A and b into x? I don't think you've done that by your analytical method. – Mark L. Stone Mar 09 '16 at 02:25
  • Yes you're correct, I want to propagate the uncertainties of A and b into x. In the first half of the script however, I pretended I know A and x, to calculate b so I have a system where everything is known so I can verify my implementation of your method is working. In the second half of the script, I ignored x, and solved for x. If my implementation was correct, the uncertainty of x should match. They don't. In the first half I created x to have uncertainty of $\pm 0.02$. The linear least squares gave me a distribution of x with uncertainties of 0.3 etc – Pete Mar 09 '16 at 02:37
  • I'll comment the script better and come back. – Pete Mar 09 '16 at 02:42
  • [This is a better commented version](http://pastebin.com/gu0UQgSU). – Pete Mar 09 '16 at 17:12