Different fitmodels for different intervals A number of points were measured and plotted on an XY graph. I want to have the best fit line for this graph.
The best fit model for the first part of the graph is described by some function $f_{1}\left(x\right)$ from $ 0 \leq x \leq a$.
The best fit model for the second part of the graph is described by function $f_{2}\left(x\right)$ from $ a < x$.
Another problem is that $a$ is unknown.
How can I obtain the best fit line for this graph with different fit models?
 A: You can solve this problem with the minimisation of errors.
In the following I will assume that your measurement errors are normally distributed. That means that given the function $f_1$, the probability of measuring a value $y$ at the argument $x$ is
$$p(y)  = \mathcal{N} \exp\left[-\frac{1}{2\sigma^2}(y-f_1(x))^2\right]$$
where $\mathcal{N}$ is a normalisation constant. The error is given by $\sigma$.
I will now summarise how the minimisation of errors work, and the explain how this relates to your question.
Imagine you want to fit a line with points $x_1,x_2,...,x_n$ and $y_1,y_2,...,y_n$ (which we will call our data), and we assume that all these points are independent. Furthermore, we denote with $\Lambda$ all the parameters that describe your function $f_1$. As an example, if $f_1$ is a linear function $f_1(x) = a x +b $, the parameters are $\Lambda = \{a,b\}$.
Your best fit is then given by (using Bayes theorem)
$$p(\Lambda|\mathrm{data})  = \frac{p(\Lambda)}{p(\mathrm{data})}p(\mathrm{data}|\Lambda)$$
For simplicity, we will assume a flat prior on these parameters $\Lambda$, and drop the prefactor $\frac{p(\Lambda)}{p(\mathrm{data})}$. In principle you can find this by normalising your function.
Now, we use the assumption that the data are independent, we can thus factorize the likelihood as
$$p(\Lambda|\mathrm{data})  = \mathcal{N'} \prod_{i=1}^np(x_i,y_i|\Lambda)$$
Now, we can use the likelihood we have written above, to rewrite our expression as
$$p(\Lambda|\mathrm{data})  = \mathcal{N''} \prod_{i=1}^n \exp\left[-\frac{1}{2\sigma^2}(y_i-f_1(x_i))^2\right] $$
$$p(\Lambda|\mathrm{data})  = \mathcal{N''} \exp\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-f_1(x_i))^2\right] $$
Now, if you want to obtain your best fit parameters, (for example the parameters that maximise your log likelihood), you can derive the function we have above with respect to the parameter and obtain a value.
For example, if $f_1$ constant, say $f_1 = c$. Then, we derive the above logarithm of the expression with respect to $c$ and obtain
$$ \partial_c(\log p(\Lambda|\mathrm{data})) = \frac{1}{\sigma^2}\sum_{i=1}^n (y_i-c)$$
To obtain the maximum, we set this expression 0, this gives
$$ c = \frac{1}{n}\sum_{i=1}^n y_i$$
So the best fit parameter would be the average of all values.
Now, in your case you have to functions that are characterised by some unknown number of parameters. For each of the parameters, you can obtain the parameters that maximize the log likelihood.   In the same way I have shown above in an example.
The trick is to realise that the function $f_1(x)$ has to be replaced by
$$f_1(x) \rightarrow f_1(x) \Theta(-x+a) + f_2(x) \Theta(x-a) 
$$
where $\Theta$ is the Theta Heaviside function.
Now, $a$ is also part of the parameters you are trying to optimise such that the (log) likelihood is large.
I expect that this is difficult to do analytically (especially given that the Heaviside function is not differentiable), but you can rely on numerical methods to maximise your likelihood.
If your parameter space is small enough, you can evaluate the above log likelihood on a grid of the parameter space.
To be explicit, your log likelihood would look like
$$ \log p(\Lambda|\mathrm{data}) = \frac{1}{2\sigma^2}\sum_{i=1}^n \left[y_i-\left(f_1(x) \Theta(-x+a) + f_2(x) \Theta(x-a) \right)\right]^2$$
where $\Lambda$ denotes all your parameters, including the parameter $a$ that defines where $f_1$ or $f_2$ is applied.
For example, if the two functions are constant, we have in total 3 free parameters. If you choose a grid spacing of $n$ per parameter, you will have to evaluate the log likelihood for $n^3$ points. The point that gives you the largest log likelihood is then your best fit point.
For larger parameter space, this evaluation of a grid quickly becomes unfeasible. You have to rely then on other methods like Monte Carlo Markov Chains (MCMCs) for example. Depending on the complexity of the problem, you can write your own likelihood and then run an MCMC analysis with bilby. That is just one software of many that allows for effective likelihood characterisation.
A: If you don't know the breakpoint $a$, then you have two fits to do: one on the piecewise fit and a second on the actual breakpoint. I think it would be better to break them into separate tasks with finding the fit for each breakpoint.
A naive algorithm would be something like,
def fit_at_breakpoint(break: int, x: array, y: array):
    x1, y1 = x[0:break], y[0:break]
    x2, y2 = x[break:size(x)], y[break:size(y)]
    fit1 = fitter(x1, y1)
    fit2 = fitter(x2, y2)
    # compute L^2 norm:
    l2 = sum([(yi - fit1(xi))**2 for xi, yi in zip(x1, y1)])
    l2 += sum([(yi - fit2(xi))**2 for xi, yi in zip(x2, y2)])
    return sqrt(l2)


def fit_breakpoints(x: array, y: array):
    min_break = (0, 100)
    for i = 1 to size(x):
        # with appropriate protections for sizes of fits
        l2_fit = fit_at_breakpoint(i, x, y)
        if l2_fit < min_break[1] then
            min_break = (i, l2_fit)
    return min_break

So basically you iterate over each index and compute the $L^2$ norm (or other relevant norm) for the fits when broken at that point. Whichever index minimizes the $L^2$ norm is your $a$ and the two fits are your best approximations to $f_1(x)$ and $f_2(x)$.
You might be able to improve the convergence by using your eyes and finding a rough breakpoint estimate and starting your fit_breakpoints function at that index, rather than 1.
Extending the two functions above to handle multiple breaks could be done by adding an additional loop over the inner and outer breaks, though this would be much more difficult a task. It probably would be better to use some multi-dimensional optimization method at that point, than the above naive approach.
