How to constrain cumulative Gaussian parameters so that the function will intersect one given point? I am analyzing data from one study where participants had to choose (between two stimuli) the one with higher intensity. One way to look at the data is to fit the proportion of correct choices as a function of the absolute difference between the 2 intensities (let's call it delta $\Delta$). 
This gives me a function that predict the probability p of correct choice for any value of  $\Delta$. I use a cumulative Gaussian function, scaled so that the predicted probability is between $0.5$ (chance, meaning pure guessing) and $1$. Below there is a plot of a sample of data.
The problem is that often the best fitting function predict a probability higher than $0.5$ even at $\Delta=0$, which does not make much sense. The function should approach $0.5$ as $\Delta$ approach $0$, because for very small values of $\Delta$ participants are necessarily at chance (that is $lim_{\Delta\to0}p=0.5$). 
I would like to constrain the parameters so that the function will pass in $(0, 0.5)$, but I don't know how to do it. Can anyone help with this?
Any help is appreciated, thanks!
Please note that the desired function needs to predict a probability between 0.5 and 1 (probability of correct choice cannot be less than chance here), and that it needs to have both location and scale as free parameters. This is required to account also for cases where (differently from the plot below, which is presented only to demonstrate the problem stated above) the proportion of correct choices stays at chance until quite large values of $\Delta$ and then increase rapidly (with slope similar as the plot below or even more steep). 
It doesn't have to be necessarily a cumulative Gaussian, also another sigmoidal function could work (if using another function makes it easier to put the constraints).


I use R, and this is the code I use for fitting. This gives the cumulative Gaussian function, scaled so that the lower asymptote is at 0.5:
pnorm2AFC <- function(x,...){
    0.5 + 0.5*pnorm(x,...)
}

and I use this to compute the negative log-likelihood. d is a dataframe with the number of correct and error responses (nnyes and nno) for each values of delta. p here indicates the parameters (mean and standard deviation).
lnorm_2AFC <- function (p, d) {
    -sum(d$nyes * log(pnorm2AFC((d$delta - p[1])/p[2])) 
        + d$nno * log(1-pnorm2AFC((d$delta - p[1])/p[2])))
}

Then I find the parameters using optim()
par <- optim(par = c(0.2, 0.2), lnorm_2AFC, d= data)

 A: Why wouldn't you just fit a probit GLM without intercept?
Then $E(Y) = \Phi(0+b\Delta)$ and 
$E(Y|\Delta=0) = \Phi(0)=0.5$. 
Isn't that a solution to the posted problem?

Here's a plot for a model I fitted by transforming $p$ using an inverse normal cdf ($y=Φ^{-1}(p)$), then fitting a linear regression with no constant term, then transforming back:

It, too satisfies the conditions mentioned in the question
A: In the meantime I found one solution to my question so I will post it here. 
One way to constrain the function to pass in $(0, .5)$ is to use a cumulative Weibull function instead of a cumulative Gaussian. The Weibull is defined only for positive values, so that at $0$ the predicted probability is necessarily $0$. When scaled as described in the question above, the predicted probability at $0$ will be necessarily $0.5$. 
A: I've been bothered by this problem for several years, so happy to find out I'm not the only one who cares about this. I just came up with this solution after read your discussion with Mark L. Stone:

*

*simply add the constraint of mu >= 3*sigma when finding the maximum likelihood. this is because ±3 standard deviation covers 99.7% of the area under the normal distribution. passing exactly the (0,0.5) point might not be possible (Glen_b does have a point). so, I feel this is close enough to what we want.

below is how I fit my psychometric function in matlab, using the maximum likelihood method
    %x(1)=mu, x(2)=sigma, x(3)=lamda
    %constraint: 3*sigma-mu<=0
    %data=[intensity;ncorrect;Ntotal]
    %likelihoodfunction = prod(binopdf(ncorrect,Ntotal,cdf('Normal',intensity,mu,sigma).*(1-2*lamda)/2+1/2))
    options = optimoptions(@fmincon,'Algorithm','sqp','MaxIterations',1500,'Display','iter-detailed');
    negloglikelihood=@(x) -sum(log(binopdf(data(:,2),data(:,3),cdf('Normal',data(:,1),x(1),x(2)).*(1-2*x(3))/2+1/2)));
    A=[-1,3,0];
    b=0;
    x = fmincon(negloglikelihood,start,A,b,[],[],[],[],[],options);

below is an example figure of the fit

