Rational function regression without poles in a interval, or polynomial regression with positivity constraint

Question

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?

Answer 1

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.

Answer 2

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.