How to get a confidence interval on parameters that were fitted using multiple functions and datasets at once? I use matlab a lot to fit functions to datasets. I then need to determine some kind of error estimate on the fitted parameters. 
Say I have a dataset of $N$ points which is given by the points $(x_i,y_i),\ i=1,\ldots,N$. I need to fit a function $f(x,P)$ where $P$ is the parameter matrix
$$P = \left[\begin{array}{cc} a_1\\ a_2\\\vdots\\a_k\end{array}\right]$$
where $k$ is the number of parameters to be fitted. I use the root-mean-squared method to find estimates for parameters $a_1,\ldots,a_k$. I build an error function:
$$\varepsilon(P,x,y)= \frac{1}{N}\left(\sum\limits_{i=1}^N (f(x,P)-y)^2\right)^{\frac{1}{2}}$$
and then use the matlab function fminsearch to get the minimum of this function by optimizing $P$. (I know there are better matlab functions, but bare with me).
To get the estimates I calculate the Jacobian matrix $J_f$:
$$J_f = \left[\begin{array}{c}\frac{\partial f}{\partial a_1}(P',x) \\ \frac{\partial f}{\partial a_2}(P',x)\\ \vdots\end{array}\right]$$
where $P'$ denotes the fitted parameter matrix (best values to minimize the error function $\varepsilon$). I also calculate the variance on the residuals:
$$ \sigma_r = \frac{1}{N-k}\sum\limits_{i=1}^N r_i^2$$
with
$$ r_i = y_i - f(P',x_i), $$
the residuals. I can get the covariance matrix now:
$$ \operatorname{cov}(P) = \sigma_r(J_f^TJ_f)^{-1} $$
where the upper script $T$ stands for transposing and upper script $-1$ stands for the inverse matrix. TO get the error estimates I then need to get the point-wise square root of the covariance matrix:
$$s_P = \left[ \begin{array}{c} s_{a_1} \\ s_{a_2} \\ \vdots \end{array}\right] = \sqrt{\operatorname{cov}(P)}$$
The Question
Given the above method is at least a bit correct (I assume it is). Is there a similar approach for the case where not one function $f$ is fitted, but when there are multiple functions $f_1$ and $f_2$?
I'll explain what I mean by this. Let's say I have two data sets $(x_i,y_i)$ and $(x'_i,y'_i)$ which comply respectively to functions $f(x,P)$ and $f'(x',P)$ where $P$ is the parameter matrix. While $f$ and $f'$ are different functions, $P$ should be the same for both fits. Eg $f$ might be:
$$f(x,P)=a_1(1-exp(-x/a_2)exp(x/a_3)$$
while $f'$ is:
$$f'(x',P)=a_4(1+exp(x'/a_2))exp(x'/a_4)$$
I would therefore fit simultaneously with the error function:
$$\varepsilon(P,x,y) =  \frac{1}{N}\left(\sum\limits_{i=1}^N (f(x,P)-y)^2\right)^{\frac{1}{2}} + \frac{1}{N'}\left(\sum\limits_{i=1}^N (f'(x',P)-y')^2\right)^{\frac{1}{2}}$$
and using fminsearch again. This gives me the best estimate for $P$ using all data available (I actually have 4 datasets and 2 functions). How would I estimate the error/certainty bounds on the parameters in this case? (Preferrably using Jacobian again).
My english isn't top-notch, nor is my knowledge of statistics and therefore I may be mixing nomenclature a bit. Forgive me for this.
 A: The easiest thing for you to do would be to string the data sets together (rbind() in R jargon, append in Stata jargon), so that you have variables $x_1,y_1$ coming from the first data set, $x_2,y_2$ coming from the second data set, and an identifier $I$ that takes the value of 0 for the first data set and the value of 1 for the second data set. Construct your response variable $r=y_1 (1-I) + y_2 I$ that would correspond to $y_1$ or $y_2$ depending what data set a particular record came from. Likewise, construct the expanded fit function $f=f(x_1,P) (1-I) + f'(x_2,P) I$ . Then treat the whole thing as if it were a single data set.
Example: the first data set is
    x     y
    2.3   5.6
    1.7   4.5
    0.8   2.2

The second data set is
     x     y
     11.0  25.4
     8.3   21.2
     7.5   19.3

The combined data set is
    x1     y1    x2   y2    I   r 
    2.3   5.6     0    0    0   5.6
    1.7   4.5     0    0    0   4.5
    0.8   2.2     0    0    0   2.2
    0     0    11.0  25.4   1  25.4
    0     0     8.3  21.2   1  21.2
    0     0     7.5  19.3   1  19.3

I would expect your standard errors and confidence intervals to be wrong with this method though. You assume implicitly that the variance of the error terms is constant (a homoskedastic nonlinear regression model), and that is a strong assumption to make. You would be better off with a so-called sandwich formula, for which you need the second derivatives, as well. You can probably find all the theory in Ron Gallant's book or in Bates and White's book, but I doubt you'd be able to handle these references. Bruce Hansen's lecture notes (Ch. 6 in particular) might be helpful, too.
