# GLM link function for bimodal probit fitting?

I am trying to model a set of data I have physical reason to believe can be represented by a bimodal normal cumulative distribution function (Technically it is a bimodal log-normal CDF, but I think I can handle transforming the predictor by log(x) ). I can successfully find the best fit to the data using MLE, e.g. $$\hat{y}(x;\mu_1,\sigma_1,\mu_2,\sigma_2) = \frac{1}{2} + \frac{w}{2}erf \left(\frac{x-\mu_1}{\sqrt{2}\sigma_1} \right) + \frac{(1-w)}{2}erf \left(\frac{x-\mu_2}{\sqrt{2}\sigma_2} \right)$$ However, I also want to establish a prediction interval on the data for theoretical future measurements of y at some x. The only way I have come up with to do this is to set up a GLM with what I would interpret as a sort of bimodal probit model, e.g. $$\hat{y} = \Phi ( b_0 + b_1x ) * b_2 + \Phi ( b_3 + b_4x ) * (1-b_2)$$ However, I can't invert the CDF to get the right side to be some sort of linear function like I can with a mono-modal CDF, e.g. $$G(\hat{y}) = \Phi^{-1}(\hat{y}) = b_0 + b_1x$$ Am I thinking about this incorrectly? Is there a way to do a bimodal CDF link function? Perhaps there is a better way of determining how well a dataset fits a given bimodal CDF and establish a prediction interval that does not use a GLM?

I think your chosen function may be inherently nonlinear, though it's partially linearizable (each of $(b_0,b_1)$, $(b_3,b_4)$ and $b_2$ may be made linear in a transformed model, given the others).