0
$\begingroup$

I have some sets of experimental data for some functions $f_i$ from $I=[0,1]$ onto itself, which should satisfy the following physical constraints:

  1. $f_i(0)=1$
  2. $f_i(x) \in I= [0,1] \; \forall x \in I $
  3. $f_i(1)=1$

Of course experimental data may on occasion slightly violate the constraints, because of measurement errors. I need to perform fits to different sets of data $S_i=\{ (x_1,y_1),\ldots (x_{k_i}, y_{k_i})\}$: sample size $k_i$ is usually from 10 to 20.

Given the constraints, I thought of fitting a low degree rational function for each dataset, forcing it to pass for (0,1) and (1,1):

$$\hat{f} (x) = \frac{1+p_1x+p_2x^2}{1+q_1x+(p_1+p_2-q_1)x^2}$$

Thus I have only 3 fitting coefficients. For each data set, I fit a model of this form, using MATLAB’s nlinfit function, . Sometimes I get good results, but other times I get poles in $I$. Is there a way to guarantee that the rational function doesn't have poles in $I$ ?

There may be also an alternative solution: each value $y_i$ is actually the ratio of two positive values $z_i$ and $w_i$, such that, except for experimental errors, $ z_i \leq w_i \ \forall i$, and if $x_i=0$ or $x_i=1$, then $z_i=w_i$. That's the reason behind the constraints 1, 2, 3. So another possible fitting strategy may be to fit a polynomial $\hat{p}(x)$ to the $(x_i,z_i)$, another polynomial $\hat{q}(x)$ to the $(x_i,w_i)$, and compute

$$\hat{f}(x) = \frac{\hat{p}(x)}{\hat{q}(x)}$$

However, I need some way to make sure that $\hat{p}(x)$ and $\hat{q}(x)$ are positive in I, or at least that $\hat{q}(x)>0 \ \forall x \in I$. Some constrained optimization problem, I guess, but I don't know how to setup it... is there a way to do that, preferably in MATLAB? Thanks a lot!

EDIT: I found out how to impose that $\hat{q}(x) >0 \ \forall x \in I$. By setting $x=t_1$, $x^2=t_2$, $ 1+q_1x+(p_1+p_2-q_1)x^2 $ becomes $1+q_1t_1+(p_1+p_2-q_1)t_2=l(t_1,t_2) $ which is linear in $t_1$ and $t_2$. This function is positive in $I^2=[0,1]^2$ if it positive in the corners of $I^2$. In reality, since $t_2=x^2=<x=t_1$, we only need to check the corners of the bottom left triangle, i.e., [0,0], [1,0], [1,1], leading to conditions which corresponds to the conditions

  1. $1 >0$
  2. $1+q_1>0$
  3. $1+p_1+p_2-q_1+q_1>0$

These are two linear constraints, so I can use fmincon with linear constraints to find the fit coefficients. However, using an Optimization tool, instead than a tool from the Statistics and Machine Learning toolbox, I lose any chance to get prediction bounds and confidence intervals for my coefficients, together with the fitted coefficients. Any idea how to remedy this?

$\endgroup$

2 Answers 2

Reset to default
0
$\begingroup$

Your problem is that the denominator can show roots. But this would not be the case if $$\Delta=-4 ({p_1}+{p_2})+{q_1}^2+4 {q_1}<0$$ I suppose that the easiest way to do the work would be to compute the sum of the squares of the residuals and try to minimize it using FMINCON with this nonlinear constraint (this seems to be part of the Optimization Toolbox).

I must underline that I am not a Matlab user.

$\endgroup$
4
  • $\begingroup$ good point, but that's a bit too strong. You're imposing the constraint that that the denominator has no roots $r_1$, $r_2$ anywhere in the real line. I only need the roots to be outside $I$, and of course the weaker the constraint, the better the fit. I think there should be a way to impose my weaker constraint. I'll wait for further suggestions, otherwise I'll accept your answer. $\endgroup$ Commented May 18, 2015 at 9:40
  • $\begingroup$ You are correct ! What you should need is that the roots of denominator be outside the range $(0,1)$. Let me think about it. May be a penalty function ? $\endgroup$ Commented May 18, 2015 at 9:44
  • $\begingroup$ Why not to consider that $q_1$ is fixed at an arbitrary value ? Solve the problem and check for the denominator. Identify an area (by trial and arror) in which the solution is acceptable. This would be a manual iteration process. Otherwise, you have a complex constraint : either $\Delta<0$ either, if $\Delta>0$, $r_1<0$ and $r_2>1$ must be both satisfied. Optimizer can handle that (even if not very pleasant). Please, let me know. By the way, could you send me a data set (.TXT format). My e-mail address is in my profile. Cheers :-) $\endgroup$ Commented May 18, 2015 at 10:59
  • $\begingroup$ The fit procedure must be automated, because it will be repeated a lot of times for many different datasets. But hey, great news! I found out how to impose that $r_1$ and $r_2$ are outside $I$, using just linear constraints. See edit. Regarding the data, sorry but I can't share them.. $\endgroup$ Commented May 18, 2015 at 15:13
0
$\begingroup$

If you accept a fitting with a criterion different from the classical least mean square, the method is very simple : $$ y-1=p_1(x-x^2y) +p_2(x^2-x^2y) +q_1(x^2y-xy)$$ From data $(x_1, y_1), (x_2, y_2), ... , (x_k, y_k), ... , (x_n, y_n)$, compute :

$F_k=y_k-1$

$u_k=x_k-x_k^2y_k$

$v_k=x_k^2-x_k^2y_k$

$w_k=x_k^2y_k-x_ky_k$

Then, the linear regression : $$\hat F=p_1u +p_2v +q_1w$$ leads to the approximates of $p_1$ , $p_2$ and $q_1$.

Anyways, if the criterion of fitting is imperatively the least mean square (wich is not always the best criterion in real life) the above method will give very good first estimates, in order to start a usual itterative fitting method of non-linear regression.

$\endgroup$

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.