Weird Function.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "from pymc3.distributions.dist_math import normal_lccdf\n",
    "import numpy as np\n",
    "import theano.tensor as T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://discourse.pymc.io/t/fit-interval-as-a-model-parameter/5453"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def weird_function(x, gate1, gate2, A, B):\n",
    "    return (\n",
    "        ((gate1 < x) & (x < gate2)) * (A * (x - gate1) + B)\n",
    "        + (x < gate1) * B * (1 - np.sin((gate1 - x) * 3) / 10)\n",
    "        + (x > gate2) * (B + (gate2 - gate1) * A + A / 5 * np.sin((x - gate2) * 10) - (x - gate2) ** 2 / 4)\n",
    "    )\n",
    "\n",
    "A = 2\n",
    "B = 3\n",
    "gate1 = 3\n",
    "gate2 = 7\n",
    "xdata = np.linspace(0, 10, 50)\n",
    "ydata = weird_function(xdata, gate1, gate2, A, B)\n",
    "yerror = np.random.normal(scale=.2, size=len(ydata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x117d29cc0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuIAAAHwCAYAAADjFQoyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3df3RcZ33n8c9XhGIbK9SDQLSxBQOSSqB2lxhpCmohIEgLpQEExOwpak+PgW63XSktUJK2tLSnPTY9LUSiu3S3qIS6PcUhGBa2HH4N2QZUdqw66TpdAho1Ck6gURECIjcRYPTsHzOTOsqMNPfO3PvcH+/XOTmDNXOl7yBZ/swz3+f7mHNOAAAAAOLV47sAAAAAII8I4gAAAIAHBHEAAADAA4I4AAAA4AFBHAAAAPCAIA4AAAB4QBAHAAAAPCCIAwAAAB4QxAEAAAAPCOIAAACABwRxAAAAwAOCOAAAAODBJb4LiIKZLUu6VNLdnksBAABAtj1F0v3OuWLQCzMZxCVdunv37sLll19e8F0IAAAAsuvOO+/Ugw8+GOrarAbxuy+//PLCmTNnfNcBAACADDt8+LBuu+22u8NcS484AAAA4AFBHAAAAPCAIA4AAAB4QBAHAAAAPCCIAwAAAB4QxAEAAAAPCOIAAACABwRxAAAAwAOCOAAAAOABQRwAAADwIKtH3AMAgIAWV9Y1v7Sq8xsXtHfXJRob7NNwf6/vsoDMIogDAJBz80urmilXdXp57RH3jRYLmh4f0thgn4fKgGyjNQUAgBw7uXBOk3OVpiFckk4vr2lyrqKbFu6JuTIg+wjiAADk1PzSqq4/dYc23faP23TSdafOan5pNZ7CgJwgiAMAkFMz5eqOIbxh00mz5Wq0BQE5QxAHACCHFlfWW7ajtFJZXtPiynpEFQH5QxAHACCHwraZ0J4CdA9TUwAAyKHzGxdiva4TjFVEVhHEAQDIob27wkWAsNeFwVhFZB2tKQAA5FDYABtX8GWsIvKAIA4AQA4N9/dqtFgIdE2pWIilJYSxisgLgjgAADk1PT6kHmvvsT0mTY0PRVtQHWMVkRcEcQAAcmpssE/HJg7uGMZ7TDo+cSiWthTGKiJP2KwJAECOHRkZ0P59ezRbrqrSJACXigVNbbMpstsTTToZq8gkFaQNQRwAgJwbG+zT2GBfoFAd1USTNI1VBDpFEAcAAJJqGzjbWVU+uXBu282UjYkmxycO6ZqRA4FqSMNYRaBb6BEHAABti3qiSdLHKgLdRBAHAABti3qiSZLHKgLdRhAHAABtiWuiSVLHKgLdRhAHAABt6WSiSRBJHKsIRIGdDQAAoC1xTjTpdKwikAYEcQAA0Ja4J5qEGasIpAlBHAAAtMXXRJN2xyoCaUOPOAAAaAsTTYDuIogDAIC2MdEE6B6COAAAaBsTTYDuoUccAICMimqTY9InmrC5E2lBEAcAIGPml1Y1U642PXxntFjQdBdCchInmsTxvIFuMufaPKc2RczszBVXXHHFmTNnfJcCAEBXtBt4Ty6c0/Wn7tj2GPpG28g1IwcirDheeX3e8O/w4cO67bbbbnPOHQ56LSviAAAkWJBV3vml1R3DqCRtOum6U2d12b7dmVghzuvzRvqxWRMAgIQ6uXBOk3OVpiFckk4vr2lyrqKbFu6RJM2UqzuG0YZNJ82Wq90q1au8Pm+kHyviAAAkUNBVXifXMrC3Ulle0+LKeqw93d3uKV9cWU/F8waaIYgDAJBAQVd5/9st/xzq68wvrcYSSKPaSDm/tBq6HoI4fKM1BQCAhAmzyvuVtQdCfa3zGxdCXRdE0BabIMLWH8fzBnZCEAcAIGHCrvKGsXdXtG+OB22xCfrcw9Yf9fMG2kEQBwAgYeJcrY16ekjUGynD1s/UFCQBQRwAgIQJu1r75MKeQI8vFQuR9kl3spGyXcP9vRotFgJ9jaifN9AugjgAAAkTdrX2P7/gaeqx9h7bY9LU+FCor9OuTjZSBjE9PpSo5w20iyAOAEDChF3lPTIyoGMTB3cMpY0TJqNuz4hrI+XYYF+injfQLoI4AAAJFHaV98jIgE4cLanUIsiXigWdOFqK5Zj3ODdSJul5A+1iyzAAAAnUWOXdaeJIs1XescE+jQ32df3wnKDi3kiZlOcNtIsgDgBAQh0ZGdD+fXs0W66q0mTTY6lY0NQ2B+EM9/d6DaCNFpsgGza7sZEyjudN2Ec3EMQBAEiwtK/yTo8PaXKu0tYIwzRspIzqhFDkE0EcAIAU8L26HVYnLTZJc3Lh3LbPo3FC6PGJQ/Sioy0EcQAAEKlOW2ySIOgJoZft253o54NkIIgDAIDIpb3FJswJoQRx7IQgDgAAYpPGFptOTghtPNe0vgBBtAjiAAAA2+jkhNCvr3+HzZ1oiQN9AAAAthH2hNAv/PM3NDlXabma3tjcedPCPZ2UhxQjiAMAAGwj7Amhn/7iStubO8OuuiPdCOIAAADbCNs60ubezoc2dyJ/6BEHAADYRpgTQoPaurkzDDaEpg9BHAAAYAdBTgg1tb8afrH5pdVQwZnTPtOL1hQAAIAdNE4I7bHtH9dj0lXP6A/1NcJsCj25cI4NoSlGEAcAAGjDkZEBnThaUqlYaHp/qVjQiaMl/fjTHh/q8wfdFBr0tE82hCYPrSkAACDXgvRWt3NC6BN6HxOqjqDtI5z2mX4EcQAAkEud9FZvd0JomM2dpWIhUH94N077hH+0pgAAgNyJurd6enxox37yhh6TpsaHAn3+Tk77RHIQxAEAQK7E0VsdZHPn8YlDgVtGwp72GfY6RIMgDgAAciVMb3UY7W7uvGbkQODPHfa0z7DXIRp8NwAAQG7E3VvdzubOMMJuumSzZrIQxAEAQG500lvdSXDebnNn2M8X9YZQRI/WFAAAkBtZ6q2OekMookcQBwAAuZGl3uqoN4Qiel0J4mb2ajN7t5l9zszuNzNnZn+1wzXPNbOPm9mamT1oZmfN7Foze1Q3agIAANgqa73VUW4IRfS69fLutyX9mKTzku6V9PTtHmxmL5f0IUkbkk5KWpP0s5LeJWlM0mu6VBcAAMBDsthbHdWGUESvW0H811QL4EuSni/pllYPNLNLJf25pO9LutI59w/1j79N0mclvdrMXuuc+0CXagMAAHjI9PiQJucqbY0wTFNvdbc3hCJ6XWlNcc7d4pyrOufamcr5aklPkPSBRgivf44N1VbWJemXu1EXAADAVvRWIyl87Dx4Yf32E03uu1XSA5Kea2aPcc59J76yAACID20Efh0ZGdD+fXs0W66q0qRNpVQsaGp8iBCOSPkI4j9Sv13ceodz7oKZLUt6pqSnSrpzu09kZmda3LVtjzoAAL7ML61qplxt2qM8WixomvAXG3qr4ZuPIP64+u23W9zf+PgPxlALAACxOblwTtefuqNlb/Lp5TVNzlV0fOIQUy5iRG81fEneUMwAnHOHm328vlJ+RczlAADQ0vzS6rYhvGHTSdedOqvL9u1mZRzIOB8H+jRWvB/X4v7Gx78VQy0AAMRiplxta0qHVAvjs+VqtAUB8M7HiviXJT1b0rCkh/V4m9klkoqSLki6K/7SAADovsWV9UBzqyWpsrymxZV1WiaQavTfb89HEP+spJ+T9NOS/mbLfc+TtEfSrUxMAQBkxfzSaujrCC1IIzYlt8dHa8rNklYlvdbMnt34oJntkvQH9T++x0NdAABE4vzGhVivA6KyuLKu980v693lqt43v6zFlfVHPObkwjlNzlVavgvU2JR808I9UZebeF1ZETezV0h6Rf2PT6rfPsfMbqz/71Xn3JslyTl3v5m9QbVA/r/N7AOqHXF/tWqjDW9W7dh7AAAyYe+ucP/chr0O6LZ2V7jZlBxMt/6G/wdJv7DlY0+t/ydJX5H05sYdzrmPmNnzJf2WpFdJ2iVpSdKvS5pt84ROAABSIWzQyHNAQXIEGbt58233Bt6UnOef864Ecefc2yW9PeA185Je2o2vDwBAkg3392q0WAi0YbNULNAfDu+CrHC/9UNnFXQlNe+bkn30iAMAkDvT40PqsfYe22PS1PhQtAUh99rp9w4ydjNsO0PYzcxZQPMZAAAxGBvs07GJgzuuLvaYdHziUK7frke02u33DjN2M4w8b0omiAMAEJMjIwPav2+PZstVVZoEnFKxoCnGuiFCQfq9/+278QTkPG9Kzu8zBwDAg7HBvodWGznoBHEKOtFk4or9sdSV5xeeBHEAADwY7u8leCNWQfq9N520EENbSt43JRPEAQDoAla4kWRh+r2/svZAqK9lJrUziJpNyQRxAAA6wlHeSIOwk0meXNgTKJCXigVNXHEZm5LbxPhCAABC4ihvpEXYySQjTykEHrt5ZGRAJ46WVCoWmj6uVCzoxNGSrhk5EKqmLGFFHACAEDjKG2kSdjLJMy+7VCPF4GM32ZTcHoI4AAAhBN34lvejvOFX2J+9RnAOO3aTTcnbI4gDABBQmI1veT/KG34N9/dqtFgI9HN78UQTVrijQRAHACCgsBvf5pdWCS3wZnp8SJNzlbbeyWk10YQV7u4iiAMAEFDYjW95PsobwXV79XlssE/HJoL3eyM6BHEAAAIKu/Etz0d5o31RjsQ8MjIQut8b3cdvBAAAAupk4xuwnZML57ZdsW6MxDw+cSj0+D/6vZODIA4AQECdbnwDmol7JCb93v5xoA8AACFMjw8FPugE2E6YkZhIN4I4AAAhNDa+7RTG2fiGdnQyEjNrFlfW9b75Zb27XNX75pcz+RwbaE0BACAkNr6hWxiJGe0m1aQiiAMA0AE2vqEb8j4SM45NqklEEAcAoAvY+IZO5HkkZtybVJMk/d89AAC6jNVtxC3PIzHDbFLNwvOWCOIAADwkjz2qSIa8jsTsZJNq2p+7xNQUAAAk1XpUJ+cqLUNBo0f1poV7Yq4MeZHHkZidbFLNAoI4ACD3gvaoZiUEIFnyOBIz75tUCeIAgNzjIBUkxZGRAZ04WlKpWGh6f6lY0ImjpcxMDsnzJlWJHnEAQM7lvUcVyZOnkZh53qQqEcQBADnHQSpIqjyMxMzrJtUGWlMAALmW9x5VwLc8blJtIIgDAHIt7z2qgG953KTawG8RAECu5b1HFUiCIyMD2r9vj2bLVVWatKmUigVNZXCOP0EcAJBree9RBZIiT5tUGwjiAIDM2+kf9unxIU3OVdoaYZi1HlUgafKwSbWBIA4AyKx2j6xv9KjudKhPFntUAfjDZk0AQCYFPbI+bwepAPCPFXEAQOYEPbL+sn27H1oZz1uPKgB/COIAgMwJc2T9xe0meepRBeAPrSkAgEzp5Mh6AIgTQRwAkCmdHFkPAHEiiAMAMoUj6wGkBUEcAJApHFkPIC0I4gCATOHIegBpQRAHAGRK48j6IDiyHoAPBHEAQOZMjw+px9p7LEfWA/CFIA4AyJzGkfU7hXGOrAfgEztTAACZdGRkQPv37dFsuapKk7nipWJBU+NDhHAA3hDEAQCZxZH1AJKMIA4AyDyOrAeQRPSIAwAAAB4QxAEAAAAPCOIAAACABwRxAAAAwAOCOAAAAOABQRwAAADwgCAOAAAAeMAccQBAqnA4D4CsIIgDAFJhfmlVM+WqTjc5rn60WNA0x9UDSBlaUwAAiXdy4Zwm5ypNQ7gknV5e0+RcRTct3BNzZQAQHkEcAJBo80uruv7UHdp02z9u00nXnTqr+aXVeAoDgA4RxAEAiTZTru4Ywhs2nTRbrkZbEAB0CUEcAJBYiyvrLdtRWqksr2lxZT2iigCgewjiAIDECttmQnsKgDQgiAMAEuv8xoVYrwOAOBHEAQCJtXdXuCm7Ya8DgDgRxAEAiRV2LjjzxAGkAUEcAJBYw/29Gi0WAl1TKhY4aRNAKhDEAQCJNj0+pB5r77E9Jk2ND0VbEAB0CU10AIBEGxvs07GJgzse6tNj0vGJQ7SlADmzuLKu+aVVnd+4oL27LtHYYF9q3hUjiAMAEu/IyID279uj2XJVlSZzxUvFgqbGhwjhQI7ML61qplxtetbAaLGg6RT8TiCIAwBSYWywT2ODfale/QLQHScXzm37Ltnp5TVNzlV0fOKQrhk5EG9xARDEAQCpMtzfS/AGcmx+aXXHVjVJ2nTSdafO6rJ9uxO7Ms5mTQAAAKTGTLm6Ywhv2HTSbLkabUEd8BrEzexnzOxTZnavmT1oZneZ2QfN7Dk+6wIAAEDyLK6sN+0J305leU2LK+sRVdQZb0HczN4h6X9JukLSJyTNSLpN0sslzZvZ63zVBgAAgOSZX1qN9bqoeekRN7MnSXqzpBVJh5xz/3rRfS+Q9FlJvy/pr3zUBwAAgOQ5v3Eh1uui5mtF/Mn1r125OIRLknPuFknrkp7gozAAAAAk095d4daQw14XNV9VVSV9V9KomfU55x56v8DMniepV9JHPNUGAIgR4wgBtCvs9JOkTk3xEsSdc2tm9lZJ75T0RTP7iKRvSHqapKslfVrSL+30eczsTIu7nt6tWgEA0cjCYRwA4jXc36vRYiHQhs1SsZDYF/feNms6526QNKHai4E3SLpO0msk3SPpxq0tKwCA7Di5cE6Tc5WW/5g2DuO4aeGemCsDkHTT40PqsfYe22PS1PhQtAV1wOfUlN+QdLOkG1VbCX+spMOS7pL012b2Rzt9Dufc4Wb/SfpShKUDADoQ9DCOpE47AODH2GCfjk0c3DGM95h0fOJQot9Z8xLEzexKSe+Q9FHn3K875+5yzj3gnLtN0islfVXSm8zsqT7qAwBEJ0uHcQDw48jIgE4cLalULDS9v1Qs6MTRUqKPt5f8bdZ8Wf32lq13OOceMLPTqgXyZ6m2Qg4AyIBODuNIao8nAD/GBvs0NtiX6g3fvoL4Y+q3rUYUNj7+3RhqAQDEpJPDONLyDyuAeA3396b294OvHvHP1W/faGaXXXyHmb1E0pikDUl/H3dhAIDoZO0wDgDohK8V8ZslfUbSiyTdaWYflnSfpMtVa1sxSdc5577hqT4AQASydhgHAHTC1xzxTTN7qaRfkfRa1frB90hak/RxSbPOuU/5qA0AEJ2sHcYBAJ3wtsTgnPuepBvq/wEAciBrh3EAQCe8zREHAORTlg7jAIBOEMQBALHK0mEcANAJdr8AAGJ3ZGRA+/ft0Wy5qkqTNpVSsaCp8SFCOIBMI4gDALomyMEaWTiMAwA6QRAHAHRsfmlVM+Vq002Yo8WCprdZ3U7zYRwA0Al6xAEAHTm5cE6Tc5WWk1BOL69pcq6imxbuibkyAEg2gjgAILT5pVVdf+oObbrtH7fppOtOnQ19xD0AZBFBHAAQ2ky5umMIb9h00my5Gm1BAJAiBHEAQCiLK+uBDuaRpMrymhZX1iOqCADShSAOAAglbJsJ7SkAUEMQBwCEcn7jQqzXAUDWEMQBAKHs3RVuAm7Y6wAgawjiAIBQwp56yWmZAFBDEAcAhDLc36vRYiHQNaVigcN7AKCOIA4ACG16fEg91t5je0yaGh+KtiAASBGCOAAgtLHBPh2bOLhjGO8x6fjEIdpSAOAi7JgBALS0uLKu+aVVnd+4oL27LtHYYN8jWkuOjAxo/749mi1XVWkyV7xULGhqfIgQDgBbEMQBAI8wv7SqmXK16YE9o8WCprcE67HBPo0N9rUV3AEANQRxAMDDnFw4p+tP3dHy6PrTy2uanKvo+MQhXTNy4GH3Dff3ErwBoE30iAMAHjK/tLptCG/YdNJ1p85ySiYAdIAgDgB4yEy5umMIb9h00my5Gm1BAJBhBHEAgKTaxsxmPeHbqSyvaXFlPaKKACDbCOIAAEkK3WZCewoAhEMQBwBIks5vXIj1OgDIO4I4AECStHdXuEFaYa8DgLwjiAMAJCn0gTsc1AMA4RDEAQCSajPAR4uFQNeUigXmhgNASARxAMBDpseH1GPtPbbHpKnxoWgLAoAMI4gDAB4yNtinYxMHdwzjPSYdnzhEWwoAdIAdNgCAhzkyMqD9+/ZotlxVpclc8VKxoKnxIUI4AHSIIA4AeISxwT6NDfZpcWVd80urOr9xQXt3XaKxwT56wgGgSwjiAICWhvt7Cd4AEBF6xAEAAAAPCOIAAACABwRxAAAAwAOCOAAAAOABQRwAAADwgCAOAAAAeEAQBwAAADxgjjgA5AgH9ABAchDEASAH5pdWNVOu6nSTI+tHiwVNc2Q9AMSO1hQAyLiTC+c0OVdpGsIl6fTymibnKrpp4Z6YKwOAfCOIA0CGzS+t6vpTd2jTbf+4TSddd+qs5pdW4ykMAEAQB4AsmylXdwzhDZtOmi1Xoy0IAPAQgjgAZNTiynrLdpRWKstrWlxZj6giAMDFCOIAkFFh20xoTwGAeBDEASCjzm9ciPU6AEAwBHEAyKi9u8JNqA17HQAgGII4AGRU2LngzBMHgHgQxAEgo4b7ezVaLAS6plQscNImAMSE9x8BIKXaOa5+enxIk3OVtkYY9pg0NT4UUbUAgK0I4gCQMkGOqx8b7NOxiYM7HurTY9LxiUO0pQBAjGhNAYAUCXNc/ZGRAZ04WlKpRZtKqVjQiaMlXTNyIJKaAQDNsSIOACkR9Lj6y/btftjK+NhgX1vtLACAeBDEASAlwhxXv7XVZLi/l+ANAAlBawoApADH1QNA9hDEASAFOK4eALKHIA4AKcBx9QCQPQRxAEgBjqsHgOwhiANACnBcPQBkD0EcAFKA4+oBIHsI4gCQEtPjQ+qx9h7LcfUAkHwEcQBIicZx9TuFcY6rB4B0YBcPAKTIkZEB7d+3R7PlqipN5oqXigVNjQ8RwgEgBQjiAJAyHFcPANlAEAeAlOK4egBIN3rEAQAAAA8I4gAAAIAHBHEAAADAA4I4AAAA4AFBHAAAAPCAIA4AAAB44D2Im9m4mX3YzO4zs++Y2dfM7JNm9lLftQEAAABR8TpH3Mz+SNJbJN0r6aOSViU9QdJhSVdK+ri34gAAAIAIeQviZvYG1UL4+yW90Tn33S33P9pLYQAAAEAMvLSmmNljJP2hpHNqEsIlyTn3vdgLAwAAAGLia0X8xaq1oNwgadPMfkbSj0rakHTaOfcFT3UBAAAAsfAVxEfqtxuSblcthD/EzG6V9Grn3Ne3+yRmdqbFXU/vuEIAiNniyrrml1Z1fuOC9u66RGODfRru7/VdFgAgIr6C+BPrt2+R9EVJPynpHyUVJf2xpKskfVC1DZsAkGnzS6uaKVd1enntEfeNFguaHh/S2GCfh8oAAFHyFcQbvekXJF3tnLu7/uc7zOyVkr4s6flm9pzt2lScc4ebfby+Un5FF+sFgEicXDin60/doU3X/P7Ty2uanKvo+MQhXTNyIN7iAACR8jVH/Fv129svCuGSJOfcA5I+Wf/jaJxFAUCc5pdWtw3hDZtOuu7UWc0vrcZTGAAgFr6C+Jfrt99qcf8367e7Y6gFALyYKVd3DOENm06aLVejLQgAECtfQbwsyUl6hpk1q6GxeXM5vpIAID6LK+tNe8K3U1le0+LKekQVAQDi5qVH3Dn3FTP7mKSrJU1LelfjPjO7StJPqbZa/gkf9QFAp3aagBK2zWR+aZVJKgCQET6PuP8VSc+S9M76HPHbVZua8gpJ35f0eufctz3WBwCBtTsB5fzGhVCfP+x1AIDk8RbEnXP3mtlhSb+j2sr48yTdL+ljko455077qg0AwggyAWXvrnC/fsNeBwBIHl894pIk59zXnXP/xTn3ZOfcDzjn+pxzrySEA0iboBNQ9vzAo0J9HeaJA0B2eA3iAJAVQSegnLrtqxotFgJ9jVKxQH84AGQIQRwAOhR2AsqrrrhMPdbe43tMmhofClEdACCpCOIA0KGwE1Ae+O73dWzi4I5hvMek4xOHaEsBgIxh1w8AdKiTCSi/OFbU/n17NFuuqtJkVb1ULGiqPmkFAJAtBHEA6FCnE1DGBvs0Nti34+xxAEC2EMQBoENhV6u3Xjfc30vwBoAcoUccADo03N/LBBQAQGAEcQDogunxISagAAACIYgDQBeMDfYxAQUAEAg94gDQJUdGBpiAAgBoG0EcALqICSgAgHYRxAEgAkxAAQDshB5xAAAAwAOCOAAAAOABQRwAAADwgCAOAAAAeEAQBwAAADwgiAMAAAAeEMQBAAAADwjiAAAAgAcc6AMAbeCkTABAtxHEAWAb80urmilXdXp57RH3jRYLmh4f0thgn4fKAABpR2sKALRwcuGcJucqTUO4JJ1eXtPkXEU3LdwTc2UAgCwgiANAE/NLq7r+1B3adNs/btNJ1506q/ml1XgKAwBkBkEcAJqYKVd3DOENm06aLVejLQgAkDkEcQDYYnFlvWU7SiuV5TUtrqxHVBEAIIsI4gCwRdg2E9pTAABBEMQBYIvzGxdivQ4AkE8EcQDYYu+ucJNdw14HAMgngjgAbBF2LjjzxAEAQRDEAWCL4f5ejRYLga4pFQuctAkACIQgDgBNTI8Pqcfae2yPSVPjQ9EWBADIHII4ADQxNtinYxMHdwzjPSYdnzhEWwoAIDB2FgFAC0dGBrR/3x7NlquqNJkrXioWNDU+RAgHAIRCEAeAbYwN9mlssE+LK+uaX1rV+Y0L2rvrEo0N9tETDgDoCEEcANow3N9L8AYAdBU94gAAAIAHBHEAAADAA4I4AAAA4AFBHAAAAPCAIA4AAAB4QBAHAAAAPGB8IYDUY8Y3ACCNCOIAUmt+aVUz5apONzn1crRY0DSnXgIAEowgDiCVTi6c0/Wn7tCma37/6eU1Tc5VdHzikK4ZOfCw+1hBBwAkAUEcQOrML61uG8IbNp103amzumzfbo0N9rGCDgBIFDZrAkidmXJ1xxDesOmk2XJVJxfOaXKu0jSES/++gn7Twj1drBQAgNZYEQeQKosr6y3DdCuV5TWdvntNLuAKOgAAUWJFHECqzC+thrpupxDe0FhBBwAgagRxAKlyfuNC5F+jsrymxZX1yL8OACDfCOIAUmXvrng66sKuvAMA0C6COIBUiat3O46VdwBAvhHEAaTKcH+vRouFyL9OXCvvAID8IogDSJ3p8SH1WHuPbfNhj8DUFABA1AjiAFJnbLBPxyYO7hjGe0x6x6sOBV5BLxULnLQJAIgcQRxAKh0ZGdCJoyWVWoTsUrGgE+NMKHsAABSWSURBVEdLumbkQKAV9B6TpsaHulgpAADN0QQJILXGBvs0NtinxZV1zS+t6vzGBe3ddYnGBvsetqLdWEG//tQd257I2WPS8YlDtKUAAGJBEAeQesP9vTu2khwZGdD+fXs0W66q0uRkzlKxoKnxIUI4ACA2BHEAudHuCjoAAHEgiAPInXZW0AEAiBqbNQEAAAAPCOIAAACABwRxAAAAwAOCOAAAAOABQRwAAADwgKkpABKH8YIAgDwgiANIjPmlVc2Uqzrd5MCd0WJB0xy4AwDIEFpTACTCyYVzmpyrNA3hknR6eU2TcxXdtHBPzJUBABANgjgA7+aXVnX9qTu06bZ/3KaTrjt1VvNLq/EUBgBAhAjiALybKVd3DOENm06aLVejLQgAgBgQxAF4tbiy3rIdpZXK8poWV9YjqggAgHgQxAF4FbbNhPYUAEDaEcQBeHV+40Ks1wEAkBSJCeJm9jozc/X/Xu+7HgDx2Lsr3BTVsNcBAJAUiQjiZnZA0p9KOu+7FgDxCjsXnHniAIC08x7EzcwkvU/SNyT9medyAMRsuL9Xo8VCoGtKxQInbQIAUs97EJc0JemFkn5R0r95rgWAB9PjQ+qx9h7bY9LU+FC0BQEAEAOvQdzMLpd0XNKMc+5Wn7UA8GdssE/HJg7uGMZ7TDo+cYi2FABAJnjb7WRml0g6IemcpN8M+TnOtLjr6WHrAuDHkZEB7d+3R7PlqipN5oqXigVNjQ8RwgEAmeFz7MDvSHqWpJ9wzj3osQ4ACTE22KexwT4trqxrfmlV5zcuaO+uSzQ22EdPOAAgc7wEcTMrqbYK/ifOuS+E/TzOucMtPv8ZSVeE/bwA/Bru7yV4AwAyL/Ye8XpLyl9KWpT0tri/PgAAAJAEPjZr7pU0LOlySRsXHeLjJP1u/TF/Xv/YDR7qAwAAACLnozXlO5LmWtx3hWp945+X9GVJodtWAAAAgCSLPYjXN2Y2PcLezN6uWhB/v3PuvXHWBQAAAMQpCQf6AAAAALnjc3whgBxgFCEAAM0lKog7594u6e2eywDQBfNLq5opV3W6yeE8o8WCpjmcBwCQc7SmAOi6kwvnNDlXaRrCJen08pom5yq6aeGemCsDACA5COIAump+aVXXn7pDm277x2066bpTZzW/tBpPYQAAJAxBHEBXzZSrO4bwhk0nzZar0RYEAEBCEcQBdM3iynrLdpRWKstrWlxZj6giAACSiyAOoGvCtpnQngIAyCOCOICuOb9xIdbrAABIM4I4gK7ZuyvcRNSw1wEAkGYEcQBdE3YuOPPEAQB5RBAH0DXD/b0aLRYCXVMqFjhpEwCQS7wfjMTjiPR0mR4f0uRcpa0Rhj0mTY0PRV8UAAAJRBBHYnFEejqNDfbp2MTBHQ/16THp+MQhvocAgNyiNQWJxBHp6XZkZEAnjpZUatGmUioWdOJoSdeMHIi5MgAAkoMVcSRO0CPSL9u3m1XVBBob7NPYYB+tRQAAtEAQR+KEOSKdIJ5cw/29BG8AAJqgNQWJwhHpAAAgLwjiSBSOSAcAAHlBawpit13PMEekAwCAvCCIIzbtjCP0dUQ6GwoBAEDcCOKIxcmFc9tOQmmMI7z2RcOhPn/YzZpxzirPStjPyvMAAMA3gjgiF2Qc4Q2fWdTTn9SrL93X/ubLsEekt/vi4PjEoY7mXWflYKKsPA8AAJKCzZqIXNBxhFLt1MV2hD0iPeis8rCbQbNyMFFWngcAAElCEEekwowj/NJ965p+0dCOYbyTI9LDzCrfanFlXe+bX9a7y1W9b375ESMU4wr7UcvK8wAAIGloTUGkwoayS3c9WieOljRbrqrSJMiXigVNhWyF6GRW+XB/b9stGlk5mCgrzwMAgKQhiHcJG9ia62QcYVRHpHcyq/z2c99se9NpJ2E/KTp90QIAAFojiHeIDWzb68Y4wm4fkR72xcH/++r9OnX7vW21aLzr04uhvsb80mqiAmwnL1qS9DwAAEgigngH4pq6kWZhX4RE+eIl7IuDhbvX2m7RaPNhj9DpwUTdfveAA5YAAIgOQTykoBvYLtu3O5cr48P9vRotFgK1N4QdR3ix7QJp2O/DV9Ye6KimdoR9kRDVOzO+DlgCACAP+NcyJDawtW96fEiTc5W2/v8KO46wod1AGvTFwZMLe2IJ4mF+RqJ8ZyaJ72gAAJAVjC8MoZMNbHk0NtinYxMHIx1HKAWbdT09vvN4xIvrGikWQtUURJh3AqIeLdh4RyOIbryjAQBAHhDEQ+hkA1teHRkZ0ImjJZVahLpSsaATR0uhe+mDBlJJgV4cPPOHLw1VV5tZP/Q7Ad2Yh76ToC9aOnlHAwCAPKE1JQQ2sIUT1ThCKVwgPflLz9H+fXvamlUe9t2Ma188pJnPbF9b2HcC4hot2HhHY6cXOp2+owEAQN4QxENgA9vDBQ3W3R5H2EkgbffFQdhNp9Pjw3r2kwuRHEwU52jBIyMDbb9oAQAA7clmMowYG9hqkjJDvRuBtJ0XB2E3nUb1TkDc78xE+Y4GAAB5RBAPwddIviRJ0gz1uAJppy0a3X4nwNc7M91+HgAA5BVBPKQ4R/I1JGUlMmkz1OMMpElq0eCdGQAA0o0gHlKcG9iS0gLSkLQZ6nEH0qS0aPDODAAA6UYQ70Acq6NJagGR4pvUEYSvQJqEFg0f78wAAIDuIIh3KMrV0aS1gDRqCntdlKE1r4GU0YIAAKQXQbxLolgdTVoLiJTcGep5DqRJ6lsHAADtI4gnVBJbQKRkz1DPcyBNSt86AABoH0E8oZLaApL0SR15D6RJ6FsHAADtIYgnlK8WkKhOmIw7HBJIAQBA0hHEEyruFpAgIxLzujESAACgm3p8F4Dm4mwBOblwTpNzlZar3I0RiTct3PPQ1zg2cVA9tv3nzeLGSAAAgG4hiCdUowUkiDAtIEFHJDZ614+MDOjE0ZJKLWosFQs6cbQUy2xzAACANKI1JcHiaAHpZERi3jdGAgAAdIIgnmBRz8bu1ohENkYCAAAERxBPuE5nY2+3Wp3UEYkAAAB5QBD3qN2WjjAtIO1MQUnqKZmID21FAAD4QxD3IMiowIu12wJycuHctu0sjSkoLzv0w4Frl+I5JRPRCvszCAAAuoepKTELOiowqCBTUD72f78W6msQ0NIt6p9BAADQHoJ4jMKOCgwiyBQUJ6k34Oq2j1My0T1x/AwCAID2EMRjFGZUYBBhpqCsb1yQ7XAwTwOnZKZf1D+DAACgfQTxmHQyKrBdYVcvX3bohzglMwfi+BkEAADtI4jHpJNRge0KO81k+Im9nJKZA3H8DAIAgPYx/iImcYwKDDvNpDG2jlMyk6eb3w/GVQIAkCwE8Zh0EpLbFbZt5OLrOCUzGaIYLxjHzyAAAGgfrSkx6UZI3slwf69GW7SXtMIUlOSJarxgHD+DAACgfQTxmMQVkqfHh3bceNnAFJTkiXK8IC/UAABIFoJ4jOIIyWODfTo2cZApKCkV9XhBXqgBAJAcBPEYxRWSj4wMMAUlheIYL8gLNQAAkoNdWDE7MjKg/fv2aLZcVaVJ6CoVC5oKsRFvK6agpE8n4wWDfE/j+hkEAADbI4h7EGdIZgpKesQ5XpAXagAA+EcQ94iQjIv5GC/IzyAAAP7QIw4kBOMFAQDIF4I4kBCMFwQAIF8I4kCCMF4QAID8IIgDCcJ4QQAA8oPNmkDCMF4QAIB8IIgDCcR4QQAAso8gDiQY4wUBAMguesQBAAAAD7wEcTN7vJm93sw+bGZLZvagmX3bzD5vZkfNjBcIAAAAyDRfrSmvkfQeSf8i6RZJ5yT1S5qQ9F5JLzGz1zjnnKf6AAAAgEj5CuKLkq6W9LfOuc3GB83sNyWdlvQq1UL5h/yUBwAAAETLSwuIc+6zzrmPXRzC6x+/T9Kf1f94ZeyFAQAAADFJYi/29+q3F7xWAQAAAEQoUeMLzewSST9f/+Mn2nj8mRZ3Pb1rRQEAAAARSNqK+HFJPyrp4865T/ouBgAAAIhKYlbEzWxK0pskfUnSZDvXOOcOt/hcZyRd0b3qAAAAgO5KxIq4mf2qpBlJX5T0AufcmueSAAAAgEiZ71HdZnatpHdJ+idJ4865f+3C5/zG7t27C5dffnnH9QEAAACt3HnnnXrwwQfXnHOPD3qt1yBuZm9VrS/8HyW92Dm32qXPuyzpUkl3d+PzBdDYJPqlmL8u4sX3OR/4PucD3+fs43ucDz6/z0+RdL9zrhj0Qm9B3MzeJun3JZ2RdFUW2lEaU1xa9a4jG/g+5wPf53zg+5x9fI/zIa3fZy+bNc3sF1QL4d+X9DlJU2a29WF3O+dujLk0AAAAIBa+pqY0lu4fJenaFo/5O0k3xlINAAAAEDNfR9y/3TlnO/x3pY/aAAAAgDgkYnwhAAAAkDcEcQAAAMAD73PEAQAAgDxiRRwAAADwgCAOAAAAeEAQBwAAADwgiAMAAAAeEMQBAAAADwjiAAAAgAcEcQAAAMADgngXmNl+M/sLM/uamX3HzO42sxvMbJ/v2tA5M3u8mb3ezD5sZktm9qCZfdvMPm9mR82Mv0cZZmavMzNX/+/1vutB95jZeP3v9X31391fM7NPmtlLfdeGzpnZz5jZp8zs3vrv7bvM7INm9hzftSEYM3u1mb3bzD5nZvfXfx//1Q7XPNfMPm5ma/Xv/1kzu9bMHhVX3e24xHcBaWdmT5P095KeKOl/SvqSpFFJ05J+2szGnHPf8FgiOvcaSe+R9C+SbpF0TlK/pAlJ75X0EjN7jeN0rMwxswOS/lTSeUl7PZeDLjKzP5L0Fkn3SvqopFVJT5B0WNKVkj7urTh0zMzeIek3JH1D0kdU+/4OSnq5pFeZ2c8757YNckiU35b0Y6r9Lr5X0tO3e7CZvVzShyRtSDopaU3Sz0p6l6Qx1f5dTwRO1uyQmX1S0lWSppxz777o4++U9GuS/rtz7j/5qg+dM7MXSnqspL91zm1e9PEnSTot6YCkVzvnPuSpRETAzEzSpyUVJZ2S9GZJb3DOvddrYeiYmb1B0v+Q9H5Jb3TOfXfL/Y92zn3PS3HoWP1381clfV3SIefcv1503wskfVbSsnPuqZ5KRED179u9kpYkPV+1RbG/ds69rsljL60/7nGSxpxz/1D/+C7VvvfPkfQfnXMfiKn8bfGWegfqq+FXSbpb0n/dcvfvSvo3SZNm9tiYS0MXOec+65z72MUhvP7x+yT9Wf2PV8ZeGKI2JemFkn5Rtb/LyAAze4ykP1Ttna1HhHBJIoSn3pNVyzeVi0O4JDnnbpG0rtq7H0gJ59wtzrlqm+88v1q17+8HGiG8/jk2VFtZl6RfjqDMUAjinXlB/fZTTULauqR5SXsk/XjchSE2jX+wL3itAl1lZpdLOi5pxjl3q+960FUvVu0f6VOSNut9xG81s2l6hzOjKum7kkbNrO/iO8zseZJ6JX3GR2GIxQvrt59oct+tkh6Q9Nz6i3Lv6BHvzI/Ubxdb3F9VbcV8WFI5looQGzO7RNLP1//Y7C88Uqj+fT2h2orpb3ouB903Ur/dkHS7pB+9+E4zu1W1VrOvx10YusM5t2Zmb5X0TklfNLOPqNYr/jRJV6vWcvZLHktEtFpmM+fcBTNblvRMSU+VdGechTVDEO/M4+q3325xf+PjPxhDLYjfcdX+Ef+4c+6TvotB1/yOpGdJ+gnn3IO+i0HXPbF++xZJX5T0k5L+UbW9AH+s2uLJB0W7Wao5524ws7sl/YWkN1x015KkG7e2rCBTUpXNaE0BQjCzKUlvUm1KzqTnctAlZlZSbRX8T5xzX/BdDyLR+HfvgqSrnXOfd86dd87dIemVqm0Iez5tKulmZr8h6WZJN6q2Ev5Y1Sbi3CXpr+tTcwDvCOKdabyqelyL+xsf/1YMtSAmZvarkmZUW017gXNuzXNJ6IJ6S8pfqvZ25ts8l4PoNH4f3+6cu/viO5xzD0hqvLs1GmdR6B4zu1LSOyR91Dn36865u5xzDzjnblPtxdZXJb3JzJiakk2pymYE8c58uX473OL+ofptqx5ypIyZXSvp3ZL+SbUQfp/nktA9e1X7u3y5pI2LDvFxqk1BkqQ/r3/sBm9VolON39ut/hH+Zv12dwy1IBovq9/esvWO+out06rln2fFWRRi0zKb1Rdciqq9I3ZXnEW1Qo94Zxp/ya8ys54tM6Z7VRsa/4Ck/+OjOHRXffPPcdX6SV/snFv1XBK66zuS5lrcd4Vq/2h/XrVf8rStpFdZkpP0jK2/t+samzeX4y0LXdSYhtFqRGHj448YXYlM+Kykn5P005L+Zst9z1Ntmt2tzrnvxF1YM6yId8A598+SPiXpKZJ+Zcvdv6daT9oJ5xwziFPOzN6mWgg/I2mcEJ49zrkHnXOvb/afaicvStL76x876bNWhOec+4qkj0kaUO0E5IeY2VWSfkq11XImIaXX5+q3bzSzyy6+w8xeotoi2YZqp2Ije25W7STV15rZsxsfrB/o8wf1P77HR2HNcLJmh5occX+npJJqM8YXJT2XI+7Tzcx+QbUNP99XrS2l2U7su51zN8ZYFmJkZm9XrT2FkzUzwMz2q/Z7+4BqK+S3q/Z29StUWy1/LSflppeZ9ajW6/8i1Q7v+bCk+1RrO3uZJJN0rXNuxluRCMTMXqHa309JepJqL5jv0r+/6Fp1zr15y+NvVu0F1wdUO+L+atVGG94s6Zo2DweKHEG8C8zsgKTfV+1tkMdL+hfV/uL/nnPum9tdi+S7KIRt5++cc1dGXw18IIhnj5k9QbVRlVdL+iFJ96v2j/ox59xpn7Whc2b2aNXeqX6tpGeo1o6wplp/+Kxz7lMey0NAbfw7/BXn3FO2XDMm6bdUO9J+l2qjK/9Cte//96OpNDiCOAAAAOABPeIAAACABwRxAAAAwAOCOAAAAOABQRwAAADwgCAOAAAAeEAQBwAAADwgiAMAAAAeEMQBAAAADwjiAAAAgAcEcQAAAMADgjgAAADgAUEcAAAA8IAgDgAAAHhAEAcAAAA8IIgDAAAAHhDEAQAAAA8I4gAAAIAH/x/+mOn8JK2/UwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 248,
       "width": 369
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(xdata, ydata + yerror)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BoundedLine(pm.Continuous):\n",
    "    def __init__(self, A, B, sigma, lower, upper, xdata, *args, **kwargs):\n",
    "        super(BoundedLine, self).__init__(*args, **kwargs)\n",
    "        self.A = A = T.as_tensor_variable(A)\n",
    "        self.B = B = T.as_tensor_variable(B)\n",
    "        self.sigma = sigma = T.as_tensor_variable(sigma)\n",
    "        self.lower = lower = T.as_tensor_variable(lower)\n",
    "        self.upper = upper = T.as_tensor_variable(upper)\n",
    "        self.xdata = xdata = T.as_tensor_variable(xdata)\n",
    "\n",
    "    def logp(self, x):\n",
    "        A = self.A\n",
    "        B = self.B\n",
    "        sigma = self.sigma\n",
    "        lower = self.lower\n",
    "        upper = self.upper\n",
    "        xdata = self.xdata\n",
    "        n = pm.Normal.dist(mu=A * (xdata-lower) + B, sd=sigma)\n",
    "        n2 = pm.Normal.dist(mu=5, sd=30)\n",
    "     \n",
    "        return T.sum(T.switch(T.and_(T.ge(xdata, lower), \n",
    "                                     T.le(xdata, upper)),\n",
    "                              n.logp(x),\n",
    "                              n2.logp(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [upper, lower, sigma, B, A]\n",
      "Sampling 4 chains, 0 divergences: 100%|██████████| 6000/6000 [04:43<00:00,  9.33draws/s]\n",
      "The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n",
      "The acceptance probability does not match the target. It is 0.8913020215112039, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n",
      "The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n",
      "The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n",
      "The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.\n",
      "The estimated number of effective samples is smaller than 200 for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    A = pm.Normal('A', mu=2, sigma=5)\n",
    "    B = pm.Normal('B', mu=3, sigma=5)\n",
    "\n",
    "    sigma = pm.HalfNormal('sigma', .2)\n",
    "    \n",
    "    lower = pm.Normal('lower', mu=3, sd=1, testval=3)\n",
    "    upper = pm.Normal('upper', mu=7, sd=1, testval=7)\n",
    "    \n",
    "    y = BoundedLine('obs', A=A, B=B, \n",
    "                    sigma=sigma, \n",
    "                    lower=lower, upper=upper, xdata=xdata,\n",
    "                    observed=ydata + yerror)\n",
    "\n",
    "    trace = pm.sample(draws=1000, tune=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {