<p><strong>EDIT: Important note</strong>
This has been bothering me for a while now. The answers given by myself and Abraham here are correct in the sense that they add variability to x. HOWEVER: Note that you <strong>cannot</strong> simply add uncertainty in this way to cancel out the errors you have in your x-values, so that you regress against "true x". The methods in this answer can show you how adding errors to x affects your regression if you have the true x. If you have a mismeasured x, these answers will not help you. Having errors in the x-values is a very tricky problem to solve, as it leads to "attenuation" and an "errors-in-variables effect". The short version is: having unbiased, random errors in x leads to <em>bias</em> in your regression estimates. If you have this problem, check out Carroll, R.J., Ruppert, D., Crainiceanu, C.M. and Stefanski, L.A., 2006. <em>Measurement error in nonlinear models: a modern perspective</em>. Chapman and Hall/CRC., or for a Bayesian approach, Gustafson, P., 2003. <em>Measurement error and misclassification in statistics and epidemiology: impacts and Bayesian adjustments</em>. CRC Press. I ended up solving my specific problem using Carroll et al.'s SIMEX method along with PyMC3. The details are in Carstens, H., Xia, X. and Yadavalli, S., 2017. <em>Low-cost energy meter calibration method for measurement and verification.</em> Applied energy, 188, pp.563-575. It is also available on ArXiv</p>

<hr>

<p>I converted Abraham Flaxman's answer above into PyMC3, in case someone needs it. Some very minor changes, but can be confusing nevertheless.</p>

<p>The first is that the deterministic decorator <code>@Deterministic</code> is replaced by a distribution-like call function <code>var=pymc3.Deterministic()</code>. Second, when generating a vector of normally distributed random variables,</p>

<pre><code>rvs = pymc2.rnormal(mu=mu, tau=tau)
</code></pre>

<p>is replaced by </p>

<pre><code>rvs = pymc3.Normal('var_name', mu=mu, tau=tau,shape=size(var)).random()
</code></pre>

<p>The complete code is as follows:</p>

<pre><code>import numpy as np
from pymc3 import *
import matplotlib.pyplot as plt

# set random seed for reproducibility
np.random.seed(12345)

x = np.arange(5,400,10)*1e3

# Parameters for gaussian
amp_true = 0.2
size_true = 1.8
ps_true = 0.1

#Gaussian function
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true )

# add noise to the data points
noise = np.random.normal(size=len(x)) * .02 
f = f_true + noise 
f_error = np.ones_like(f_true)*0.05*f.max()

with Model() as model3:
    amp = Uniform('amp', 0.05, 0.4, testval= 0.15)
    size = Uniform('size', 0.5, 2.5, testval= 1.0)
    ps = Normal('ps', 0.13, 40, testval=0.15)

    gauss=Deterministic('gauss',amp*np.exp(-1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.)))+ps)

    y =Normal('y', mu=gauss, tau=1.0/f_error**2, observed=f)

    start=find_MAP()
    step=NUTS()
    trace=sample(2000,start=start)

# extract and plot results
y_min = np.percentile(trace.gauss,2.5,axis=0)
y_max = np.percentile(trace.gauss,97.5,axis=0)
y_fit = np.percentile(trace.gauss,50,axis=0)
plt.plot(x,f_true,'b', marker='None', ls='-', lw=1, label='True')
plt.errorbar(x,f,yerr=f_error, color='r', marker='.', ls='None', label='Observed')
plt.plot(x,y_fit,'k', marker='+', ls='None', ms=5, mew=1, label='Fit')
plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5)
plt.legend()
</code></pre>

<p>Which results in</p>

<p><a href="https://i.stack.imgur.com/mvhbd.png" rel="noreferrer">y_error</a></p>

<p>For errors in x (note the 'x' suffix to variables):</p>

<pre><code># define the model/function to be fitted in PyMC3:
with Model() as modelx:

    x_obsx = pm3.Normal('x_obsx',mu=x, tau=(1e4)**-2, shape=40)

    ampx = Uniform('ampx', 0.05, 0.4, testval=0.15)
    sizex = Uniform('sizex', 0.5, 2.5, testval=1.0)
    psx = Normal('psx', 0.13, 40, testval=0.15)

    x_pred = Normal('x_pred', mu=x_obsx, tau=(1e4)**-2*np.ones_like(x_obsx),testval=5*np.ones_like(x_obsx),shape=40) # this allows error in x_obs

    gauss=Deterministic('gauss',ampx*np.exp(-1*(np.pi**2*sizex*x_pred/(3600.*180.))**2/(4.*np.log(2.)))+psx)

    y = Normal('y', mu=gauss, tau=1.0/f_error**2, observed=f)

    start=find_MAP()
    step=NUTS()
    tracex=sample(20000,start=start)
</code></pre>

<p>Which results in:</p>

<p><a href="https://i.stack.imgur.com/ZOh71.png" rel="noreferrer">x_error_graph</a></p>

<p>the last observation is that when doing </p>

<pre><code>traceplot(tracex[100:])
plt.tight_layout();
</code></pre>

<p>(result not shown), we can see that <code>sizex</code> seems to be suffering from 'attenuation' or 'regression dilution' due to the error in the measurement of <code>x</code>.</p>
Fit a non-linear function to data/observations with pyMCMC/pyMC
<p>I am trying to fit some data with a Gaussian (and more complex) function(s). I have created a small example below. </p>

<p>My first question is, <em>am I doing it right?</em> </p>

<p>My second question is, <em>how do I add an error in the x-direction, i.e. in the x-position of the observations/data?</em></p>

<p>It is very hard to find nice guides on how to do this kind of regression in pyMC. Perhaps because its easier to use some least squares, or similar approach, I however have many parameters in the end and need to see how well we can constrain them and compare different models, pyMC seemed like the good choice for that.</p>

<pre><code>import pymc
import numpy as np
import matplotlib.pyplot as plt; plt.ion()

x = np.arange(5,400,10)*1e3

# Parameters for gaussian
amp_true = 0.2
size_true = 1.8
ps_true = 0.1

# Gaussian function
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true )

# add noise to the data points
noise = np.random.normal(size=len(x)) * .02 
f = f_true + noise 
f_error = np.ones_like(f_true)*0.05*f.max()

# define the model/function to be fitted.
def model(x, f): 
    amp = pymc.Uniform('amp', 0.05, 0.4, value= 0.15)
    size = pymc.Uniform('size', 0.5, 2.5, value= 1.0)
    ps = pymc.Normal('ps', 0.13, 40, value=0.15)

    @pymc.deterministic(plot=False)
    def gauss(x=x, amp=amp, size=size, ps=ps):
        e = -1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.))
        return amp*np.exp(e)+ps
    y = pymc.Normal('y', mu=gauss, tau=1.0/f_error**2, value=f, observed=True)
    return locals()

MDL = pymc.MCMC(model(x,f))
MDL.sample(1e4)

# extract and plot results
y_min = MDL.stats()['gauss']['quantiles'][2.5]
y_max = MDL.stats()['gauss']['quantiles'][97.5]
y_fit = MDL.stats()['gauss']['mean']
plt.plot(x,f_true,'b', marker='None', ls='-', lw=1, label='True')
plt.errorbar(x,f,yerr=f_error, color='r', marker='.', ls='None', label='Observed')
plt.plot(x,y_fit,'k', marker='+', ls='None', ms=5, mew=2, label='Fit')
plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5)
plt.legend()
</code></pre>

<p>I realize that I might have to run more iterations, use burn in and thinning in the end. The figure plotting the data and the fit is seen here below. </p>

<p><img src="https://i.stack.imgur.com/iKjXf.png" alt="Resulting figure from the code."></p>

<p>The pymc.Matplot.plot(MDL) figures looks like this, showing nicely peaked distributions. This is good, right?</p>

<p><img src="https://i.stack.imgur.com/U71o2.png" alt="enter image description here"></p>