{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Some setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import sympy"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"x,a,b,x1,B=SR.var('x,a,b, x1,B')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"f(x)=x^(a-1)*(1-x)^(b-1)/beta(a,b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The peak of $f(x)$ is at"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a - 1)/(a + b - 2)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"m(a,b)=(a-1)/(a+b-2);\n",
"m(a,b).show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see this from a picture"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAHUCAYAAACznbW8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xmc1uP+x/H3PS0zlZlbqZhqShKa6qRlKJRsEQo/jj2HClGRLOckx8Fx5FgTY186ROQgOUhO27RN0yqStNGiaNF9t07L3L8/PqeItnvmvu/r+73v1/PxuB/XiZnud6ZT77m+1xKIRCIRAQAAwLfSXAcAAABA2VDoAAAAfI5CBwAA4HMUOgAAAJ+j0AEAAPgchQ4AAMDnKHQAAAA+R6EDAADwOQodAACAz1HoAAAAfI5CBwAA4HO+KHSRSEThcFhcOwsAAPB7vih0GzZsUDAY1IYNG1xHAQAA8BxfFDoAAADsG4UOAADA5yh0AAAAPkehAwAA8DkKHQAAgM9R6AAAAHyuvOsAQFktWSJNmCB98YW0YoW0erW0c6eUkSHVqiU1aCDl5UknnigFg67TAgAQe54udPn5+crPz9fOnTtdR4HH/PCD9PLL0jvvSF9/bf/sqKOkevWkhoes1NnfvaCP69yor+Zma/hw6eefpfLlpbPOki69VPrjH6UqVdz+GgAAiJVAxAfXL4TDYQWDQYVCIWVlZbmOA4eWLpXuvVcaMsRm4C65ROrcWTr9dOnQQ//3QeGwNGWK1KaNlJWlSERasEAaNUoaNkyaOFGqWlXq1cteNWo4/SUBAFBmFDr4wpYt0gMPSE8+aY9N+/WTrruudI9Qv//efp6XXpLS0qS//lW69VYpPT32uQEASAQ2RcDzpk2TWrSwEvaXv0iLFkl9+uynzK1cKd13n417Ua+eNHCgzfZ16ybdfbfUtKk0aVLcfgkAAMQVhQ6eFYlIzzwjnXSSdMgh0syZ1tMOOeQAn7hypXT//fssdLscdpgVuy++kKpXl9q1k/r3l7Zti9kvAQCAhKDQwZO2b5e6d5d697bX5MlSbm583qtxY6mgwB7pPvKI1L69tGpVfN4LAIB4oNDBc7Zskf7v/6TXX5f+9S/piSekChXi+57ly9vs3MSJ0nffSa1aSdOnx/c9AQCIFQodPGXzZqljR2nMGOmjj6Rrrkns+594ohW5OnWktm2ljz9O7PsDAFAaFDp4RnGxzcxNny59/rl0zjml/ImqVpWuusrGUqhVSxo3zorlhRfaWXcAAHgZx5bAE3bskC6/XPrPf6RPPrFz5VzbsUPq2tXOvHvxRVvTBwCAF3n6pgikjr59peHDpfffj0GZ27pVWr7cnptmZJT6pylfXho8WMrKkq6/3nbdXn99GbMBABAHFDo498IL0tNPS889Z7c+lNnXX0stW0ozZtgBdmWQlmbZAgHpxhulzEybSQQAwEsodHCqoMCu3+rZU+rRw3WavQsEpKeekkIhqUsXm7E791zXqQAA+AWbIuDM6tU223XKKXbAr5elpUmvvmpF7tJL7TBiAAC8gkIHJ0pK7EiSHTukN9+09WpeV7689NZb0jHHSJ06HfAiCgAAEsbThS4/P1+5ubnKy8tzHQUx9vjj0siR0htv2DEhflGlip2Pt3OndMEFdggyAACucWwJEu6rr2zPwq232lVbfjRjht39ev750ttv2zo7AABc8fQMHZLPjh3SdddJDRrY3al+1bKlXU02bJj31/8BAJIfhQ4J9eij0syZdr5bGY6I27/586U2bWyMo4svlu64Q7rrLmnSpLi+FQAA+0WhQ8J88410333SnXdKJ5wQxzfatEkqLLQxzh56SGrdWrrsMumnn+L+dgAA7BWFDgkRiUi9e9vlDffd5zpN7FSoYHe9bt9uZ9SVlLhOBABIRRQ6JMR770n//a80aFAcH7U6UquWracbNcpulQAAINEodIi7TZvsrtZOnaTzznOdJj7OPlu65Rbpz3+2XbwAACQShQ5x949/2PqyhO0GPfJIO+DuyCMT9Ibm4Yelo4+WrrpKKi5O6FsDAFIchQ5xtWSJHSL85z9LRx2VoDetVk26+mobE6hSJbv14ptvpLvvTuhbAwBSHIUOcXXPPdJhh9nRHgmzerWUn29jgjVrZjOSTz4pTZ6c8LcHAKQoCh3iZuZMu/v0vvvsyqyEWbZM6tXLRgduu03Ky5O6dZO2bnUSAQCQYih0iJs//1k67jipa1fXSRKrXDnplVekRYtstg4AgHij0CEuxoyxY0oeflgqX951msRr0kTq399+/V984ToNACDZUegQc5GIPWZt1Urq3Nl1Gnf69bMZym7d7A5bAADixdOFLj8/X7m5ucrLy3MdBVEYO1aaMEH629+kQMBBgMxMqUMHGx2qWNEevc6aZZskAACIl0AkEom4DnEg4XBYwWBQoVBIWVlZruNgPyIRqX17O0x42jRHhc5j+vSRXn7ZjjOpU8d1GgBAMvL0DB38Z9w4qaDA4eycJO3cKYXDNnrA/fdLhxwi3X676yQAgGRFoUNM3X+/1KKFdP75DkN88YUUDHpmN0IwKD36qDRsmDR6tOs0AIBkRKFDzEycKI0fL917L49af+vqq6VTTpF695a2bXOdBgCQbCh0iJlHH5UaNZI6dXKdxHsCAemZZ6T586Wnn3adBgCQbCh0iIlvvpFGjJDuuENK43fVXjVrJt18sx3p8sMPrtMAAJIJf/UiJp54QjriCOmqq1wn8ba//12qVEm6807XSQAAyYRChzL78Ufp9delW26R0tNdp5HUtKn00082esyhh0oDBtgdt0VFrtMAAJIF59ChzO65Rxo4UFq2TKpa1XUa79u5U2re3Mrd+PFsIAEAlB0zdCiTrVul55+36608U+YWLbI7xxYtcp1kr8qVkx57zG7TGD7cdRoAQDKIqtANGDBAeXl5yszMVM2aNXXhhRdq/vz5+/2cwYMHKxAI/O61devWMgWHN7zzjrR2rdSrl+skvxIKSR99ZKNHdeggnX22dNddHGMCACi7qArd+PHj1bNnTxUWFurzzz/Xjh071KFDB23atGm/n5eVlaWVK1fu8crIyChTcLgXidgRHGefLTVs6DqN/zz6qLR4sfTCC66TAAD8rnw0Hzxy5Mg9fvzaa6+pZs2amjFjhtq1a7fPzwsEAjriiCNKlxCeVVQkzZhhk2GIXtOm0nXX2e0aXbrYmjoAAEqjTGvoQv97pFWtWrX9ftzGjRtVr1491alTR+eff75mzZq1348vLi5WOBze4wXvyc+X6teXOnZ0ncS/HnhA2rLFdr4CAFBapS50kUhEffv21SmnnKImTZrs8+OOO+44DR48WCNGjNDQoUOVkZGhk08+WQsWLNjn5wwYMEDBYHD3Kycnp7QxESc//WTr5266yRb5e0rt2tLjj9vocbVq2Zl0Tz1lu4QBACiNUh9b0rNnT3388ceaOHGi6tSpc9CfV1JSohYtWqhdu3YaNGjQXj+muLhYxcXFu38cDoeVk5PDsSUeMmCAzS4tXy4ddpjrNP62YYN01FHSRRdJL77oOg0AwI9KNUPXu3dvjRgxQmPHjo2qzElSWlqa8vLy9jtDl56erqysrD1e8I6dO+2okiuu8GiZ+/ln6d13bfSBzEypXz/p1Vel/fzfAgCAfYqq0EUiEfXq1Uvvv/++xowZo/r160f9hpFIRLNnz1Z2dnbUnwtv+O9/paVLpRtucJ1kH5YskS691EafuOkmuzrtvvtcJwEA+FFUha5nz54aMmSI3nrrLWVmZmrVqlVatWqVtmzZsvtjrrnmGvXr12/3j++//3599tlnWrx4sWbPnq1u3bpp9uzZ6tGjR+x+FUioV16RcnOlE090nSR5VKok3XuvNHSoNGeO6zQAAL+JqtA999xzCoVCat++vbKzs3e/3nnnnd0fs3TpUq1cuXL3j9evX68bbrhBjRo1UocOHbRixQoVFBTohBNOiN2vAgmzZo3dbtCtG1dWxdp119laur/+1XUSAIDfcJcrojJwoN1usGKFVKOG6zT7MHOm1LKlHZLXooXrNFF56y3pqqukKVOk1q1dpwEA+AV3ueKgRSL2uLVzZw+XOcmeXzZvbqPPXH651KSJ1L+/6yQAAD+h0OGgTZ8uffWVPW71tEaNbJauUSPXSaKWlmbHwYwZI02Y4DoNAMAveOSKg9ajh/Txx9J333nwMOEkUlJiE4w1atiOYgAADoQZOhyUzZttB+a11/qgzM2aJaWn2+hDaWm243X0aGnSJNdpAAB+QKHDQRkxQgqHrdB5XiQibdtmo09ddJGtpXvgAddJAAB+QKHDQXnzTalNG6lBA9dJUkNamh1fMmqUVFjoOg0AwOsodDig1aulkSPtOA0kziWX2AHOzNIBAA7E04UuPz9fubm5ysvLcx0lpQ0bZuOll7rNkWrS0qR77pE+/VQqKnKdBgDgZexyxQG1aSNVry599JHrJAdpyxZp8WK7dsGHZ9H92s6dtpauQQPpP/9xnQYA4FWenqGDe4sW2Rquq692nSQKlSpJjRv7vsxJtqP4nnvsuJjp012nAQB4FYUO+/Xmm9Ihh0idOrlOEoXvv5e6d7cxCVx2mdSwofTgg66TAAC8ikKHfYpEpCFDpIsvlipXdp0mCmvX2h1la9e6ThIT5ctLf/mL9OGH0rx5rtMAALyIQod9mj5dWrDAZ49bk9RVV0m1a0uPPOI6CQDAiyh02Kd33pEOP1w67TTXSZCeLvXtazOmS5e6TgMA8BoKHfYqEpH+/W/p//7PB1d9pYjrr5cyM6UnnnCdBADgNRQ67NW0aban4I9/dJ2kFA4/3BadHX646yQxlZkp9e4tvfSStGaN6zQAAC+h0GGv3n1XqllTatfOdZJSqF1bGjDAxiTTu7fNnj7zjOskAAAvodDhdyIRK3S+fdy6YYM0bpyNSaZ6dXv0OmiQtHGj6zQAAK+g0OF3pk/38eNWybbmnnaajUno9tutq770kuskAACvoNDhd959V6pRw6ePW1NA3bp2jMnjj0vbtrlOAwDwAgod9vDrx63ly7tOg3256y5pxQo7xgQAAAod9jBjhvTddz5+3JoicnOlzp2lRx+VSkpcpwEAuObpQpefn6/c3Fzl5eW5jpIy3n3XFt6feqrrJGVQoYLtcK1QwXWSuLrjDumbb6SRI10nAQC4FohEIhHXIQ4kHA4rGAwqFAopKyvLdZykFYlIRx8tnXmm9MILrtPgQCIR6cQT7Xy60aNdpwEAuOTpGTok1ldfSYsX2/o5eF8gYDtex4yRZs92nQYA4BKFDrt9+KHN9rRv7zpJGX35pVSnjo1J7uKLbdcr14EBQGqj0GG34cOljh3tInhf277dtoBu3+46SdyVLy/deqs0dKj9kgEAqYlCB0nS8uW2w/WCC1wnQbS6d5cqV5aeftp1EgCAKxQ6SJJGjLDZnnPPdZ0E0crKsuvAXniB68AAIFVR6CDJ1s+deqp06KGuk6A0brnFrgN77TXXSQAALlDooFBIGjs2iR63Nmxov6CGDV0nSZi6de0w6IEDpZ07XacBACQahQ4aOdL2DyRNodu1VTcz03WShLr9djt25sMPXScBACQahQ768EOpeXOb5UkKK1ZI/fql3LbPVq2kdu2kxx93nQQAkGgUuhS3fbv0ySdJNDsnST/+KD38sI0p5vbbpcmTpcJC10kAAIlEoUtx48fbGrqkKnQp7PzzbekgBw0DQGqh0KW4Dz+U6tWTmjVznQSxkJZmBw2//760bJnrNACARKHQpbBIRPrPf2xWJxBwnQaxcs01UpUq0nPPuU4CAEgUTxe6/Px85ebmKi8vz3WUpPTNN9J330nnnec6SYwddpjUrZuNKSgzU+raVXrxRWnLFtdpAACJEIhEIhHXIQ4kHA4rGAwqFAopKyvLdZyk8fjj0j33SOvWSZUquU6DWFq0yNbSvfSSdVsAQHLz9Awd4uuTT6TTT0/CMrdlizR3bkpPTzVoYI/SBw2yR+sAgORGoUtR4bA0YYLUsaPrJHEwb57UpImNKeyWW6Q5c6SCAtdJAADxRqFLUaNH2xl0557rOgni5YwzpEaNbJYOAJDcKHQp6tNPpWOPlY46ynUSxEsgYLN0w4fb5hcAQPKi0KWgSMTWzzE7l/y6dLFdr88+6zoJACCeKHQp6Msv7ZrTpC10gYBUsSKH68nOo+ve3Xa7btrkOg0AIF4odCnok0/sL/q2bV0niZPmzaXiYhuhXr1sE8ybb7pOAgCIFwpdCvrkE+nMM6X0dNdJkAhHHil17swRJgCQzCh0Kebnn6XJk5P0uJJd5s2TWrRI+WNLfu2WW+xovjFjXCcBAMQDhS7FfP65tHNnkhe6LVukWbNS+mDh32rf3o7me+op10kAAPFAoUsxn30m5eZKdeu6ToJECgSk3r2l//yHI0wAIBlR6FJIJCKNGiWdfbbrJHDhqqukrCzp+eddJwEAxBqFLoV88420fLl01lmuk8CFKlWka6+VXn5Z2rrVdRoAQCx5utDl5+crNzdXeXl5rqMkhVGj7Hi2du1cJ4mz+vWlYcNsxB5uvllau9b+8wAAkkcgEvH+QQbhcFjBYFChUEhZWVmu4/jW+efbPoHRo10ngUsdOkihkDR1quskAIBY8fQMHWKnuFgaO9b+Mk96P/4oPfGEjfidXr2koiJp2jTXSQAAsUKhSxFTpkibN6fI+rkVK6Tbb7cRv3PeeVK9elJ+vuskAIBYodCliFGjpOrVpeOPd50ErpUrJ/XoIb39trRmjes0AIBYoNCliM8/t9m5NL7ikNStmx1j8+qrrpMAAGKBv95TwJo10owZKbJ+DgelRg3pssuk556zm0MAAP5GoUsBo0fbbExKrJ+TpGBQ6tTJRuxTz552a8Snn7pOAgAoK44tSQHdu9umiLlzXSeBl0QiUl6ezdZR6gDA35ihS3K7rvtKqcet27dLq1fbiH0KBGyWbuRIaeFC12kAAGVBoUtyCxZIy5ZJZ57pOkkCffmlVLOmjdivyy+XqlWztXQAAP+Ke6EbMGCA8vLylJmZqZo1a+rCCy/U/Pnz4/22+J8xY+yYiqS/7gulUqmSdN110uDBdosIAMCf4l7oxo8fr549e6qwsFCff/65duzYoQ4dOmjTpk3xfmvIboc44QQpM9N1EnjVjTdK69ZJ//636yQAgNIqH+83GDly5B4/fu2111SzZk3NmDFD7fYxbVRcXKzi4uLdPw6Hw3HNmKxKSqzQXX+96yTwsoYN7ZH8c89JXbq4TgMAKI2Er6ELhUKSpGrVqu3zYwYMGKBgMLj7lZOTk6h4SWXuXNsbcPrprpPA63r0sJ3QX3zhOgkAoDQSemxJJBLRBRdcoJ9//lkTJkzY58ftbYYuJyeHY0uiNGiQdOed0vr1tlYqZezcKW3aJFWpYgsIcUDbt9v9rhdcwAYJAPCjhM7Q9erVS3PmzNHQoUP3+3Hp6enKysra44XojRkjnXRSipU5yUpcVhZlLgoVKtij+SFDpA0bXKcBAEQrYYWud+/eGjFihMaOHas6deok6m1T1s6d0rhx0mmnuU7iwIIF0tln24iD1r27tHmz9OabrpMAAKIV90IXiUTUq1cvvf/++xozZozq168f77eEpFmzpFAoRdfPbdhgpykz1RSVnBy7Me255+xAagCAf8S90PXs2VNDhgzRW2+9pczMTK1atUqrVq3SFg69iquxY6XKle3IEuBg9eghzZkjFRa6TgIAiEbcC91zzz2nUCik9u3bKzs7e/frnXfeifdbp7QxY6S2baWKFV0ngZ906CDVr8/GCADwm4Q8ct3b69prr433W6esbdukCRNSdP0cyiQtzQ4aHjZMWrvWdRoAwMHiLtckNG2andqRkuvnJFsM9swzNiJqXbvaGrrBg10nAQAcLApdEho7VgoGpebNXSdxpEYNqWdPGxG1GjWkSy6Rnn/ebhsBAHgfhS4JjRkjnXqqVD7uF7t51Lp1dqDaunWuk/hWjx7SwoX2ewkA4H0UuiSzZYs0eXKKr5/77ju7lPS771wn8a1TTpEaN2ZzBAD4BYUuyUyZIhUXp/D6OcREIGCzdB9+KK1Y4ToNAOBAKHRJZuxYqXp1qUkT10ngd126SOnp0iuvuE4CADgQCl2SGTfO1s+l8ZVFGQWD0pVXSi++KO3Y4ToNAGB/PP3Xfn5+vnJzc5WXl+c6ii9s2SIVFVmhS2lVqkitW9uIMrnpJnvk+vHHrpMAAPYnEIl4/9bGcDisYDCoUCikrKws13E8a/x4qX17afZsqVkz12mQLE48UapaVRo50nUSAMC+eHqGDtEpKJAOPZT1c4itHj2kzz6TFi92nQQAsC8UuiRSUGDHTZQr5zqJYzNn2jbNmTNdJ0kKl11m3yi88ILrJACAfaHQJYnt2+38uXbtXCdBsqlcWfrTn6RXX7UjcQAA3kOhSxIzZkibN7MhAvFx443SmjXS+++7TgIA2BsKXZIoKLBNnSl7fyviqlEj+2aBx64A4E0UuiRRUCCddJJUoYLrJEhWPXrYTup581wnAQD8FoUuCezcKU2cyPq53XJzpQULbETMXHSR3ULCLB0AeA+FLgl8+aUUClHodsvIkI4+2kbETHq61LWr9K9/2SHWAADvoNAlgYICqWJF6YQTXCfxiCVLpKuvthExdcMN0vr10rBhrpMAAH6NQpcECgrsNH8mpP7n55+lN9+0ETHVoIHUoYP0/POukwAAfo1C53ORiBU6HrciUW68USostCvmAADeQKHzufnzpdWrKXRInE6dpOxsNkcAgJdQ6HyuoMCu+mrTxnUSpIoKFaTu3aUhQ6QNG1ynAQBIHi90+fn5ys3NVV5enusonlVQILVoIWVmuk7iIdnZ0t/+ZiPiont3u5lk6FDXSQAAkhSIRCIR1yEOJBwOKxgMKhQKKSsry3Ucz4hEpLp17fL0xx5znQappnNnaflyu3YuEHCdBgBSm6dn6LB/339vf6Gyfu43wmHps89sRNzceKM0a5Y0bZrrJAAACp2PFRTYeMopbnN4zsKF0jnn2Ii4OeccmyHmCBMAcI9C52MTJ0pNmkjVqrlOglRUrpwdNPz223bYMADAHQqdj02cKJ18susUSGVdu0rbt0tvvOE6CQCkNgqdT61bJ82bR6GDW9nZ0gUX2GNX72+vAoDkRaHzqcmTbaTQ7UV6ut1RlZ7uOklK6NFD+vprmzEGALhBofOpSZOkI46Q6td3ncSDGje2DRGNG7tOkhJOP106+mhujgAAlyh0PjVpks3Ocf4XXEtLsyNM3n1XWrPGdRoASE0UOh/ats3O/uJx6z7MmSPVqGEjEuLaa20cPNhlCgBIXRQ6H5o5U9q6lfPn9mnHDpsq2rHDdZKUUb26dMkl9ti1pMR1GgBIPRQ6H5o4UapcWTr+eNdJgF/06GFLF8eOdZ0EAFIPhc6HJk2STjhBqlDBdRLgF6ecIuXmcnMEALhAofOZSOSXDRGAlwQCNks3fLi0cqXrNACQWjxd6PLz85Wbm6u8vDzXUTxj4UJp9WoK3X4dc4wd1HfMMa6TpJwuXWzm+NVXXScBgNQSiES8f757OBxWMBhUKBRSVlaW6zhODR5s1y2tWycdeqjrNMDvde0qjR4tLV5s970CAOLP0zN0+L1Jk+y8XMrcfixfLvXtayMSrkcPaelS6bPPXCcBgNRBofMZ1s8dhJ9+kp580kYkXF6e1Lw5myMAIJEodD6ydq00bx7nz8HbAgG7OeLjj22mDgAQfxQ6H5k82UZm6OB1V15pZyW+/LLrJACQGih0PjJpkpSdLR15pOskwP5lZkpXX22Fbvt212kAIPlR6Hxk1/q5QMB1Eo+rXl26+WYb4UyPHnYe3X/+4zoJACQ/ji3xieJiKRiUHn5Y6tPHdRrg4LRubb9v2fEKAPHFDJ1PzJ5tpe6kk1wn8YHNm6WZM22EUz16SKNGSYsWuU4CAMmNQucTU6ZIGRnS8ce7TuID33wjtWxpI5y69FI7M/Gll1wnAYDkRqHziSlTrKNUrOg6CXDwKleW/vQnuwqsuNh1GgBIXhQ6nygstPVIgN/ceKPdP/zBB66TAEDyotD5wA8/2AGtbdq4TgJEr1EjqV07bo4AgHjydKHLz89Xbm6u8vLyXEdxqrDQRmboDlJamh2Elubp394ppUcPafx4ljUCQLxwbIkP3Hmn9M47XKME/youlurUscOGn3zSdRoASD5MYfjAlCnMzsHf0tOlrl2lwYOlLVtcpwGA5EOh87ht26QZM1g/F5Wvv5YaN7YRnnH99dL69dKwYa6TAEDyodB53BdfSFu3UuiisnWrlbmtW10nwa8cfbR01llsjgCAeKDQeVxhoZ0917y56yRA2fXoYb+nv/jCdRIASC4UOo+bMkVq0cLWIAF+16mTlJ0tvfCC6yQAkFwodB7HgcJIJhUqSN26SW+8IW3Y4DoNACQPCp2H/fijtGQJ6+eidtRR0ocf2gjPuf56afNmaehQ10kAIHlQ6DxsyhQbKXRROvRQqXNnG+E5detK555rmyO8fwomAPgDhc7DCgulWrXsQFZEYdUqacAAG+FJPXpIs2ZJ06e7TgIAySHqQldQUKBOnTqpVq1aCgQCGj58+H4/fty4cQoEAr97fcMdQAc0ZYrNzgUCrpP4zA8/SHffbSM86ZxzbKaOI0wAIDaiLnSbNm1Ss2bN9Mwzz0T1efPnz9fKlSt3vxo2bBjtW6eUHTukadPYEIHkVK6craUbOtQOGwYAlE35aD+hY8eO6tixY9RvVLNmTR3KmqaDNmeOXZHE+jkkq27dpPvuk4YMkXr1cp0GAPwtYWvomjdvruzsbJ1xxhkaO3bsfj+2uLhY4XB4j1eqmTJFKl/ezqADklF2tnThhWyOAIBYiHuhy87O1osvvqj33ntP77//vo499lidccYZKigo2OfnDBgwQMFgcPcrJycn3jE9p7DQboeoVMl1Eh869FDpkkvY5eoDPXpIc+dKkya5TgIA/haIREr/vXEgENAHH3ygCy+8MKrP69SpkwKBgEaMGLHXf19cXKzi4uLdPw6Hw8rJyVEoFFJWVlZp4/rK0UdL550nPfWU6yRA/JSUSMd5S1nVAAAgAElEQVQcY2tFhwxxnQYA/MvJsSWtW7fWggUL9vnv09PTlZWVtccrlaxeLS1axIaIUtu2TVq+3EZ4WlqadOON0rvvSmvWuE4DAP7lpNDNmjVL2dnZLt7aFwoLbWRDRCl99ZWUk2MjPO/aa23817+cxgAAX4t6l+vGjRu1cOHC3T9esmSJZs+erWrVqqlu3brq16+fVqxYoddff12SNHDgQB155JFq3Lixtm3bpiFDhui9997Te++9F7tfRZIpLJQOP1yqV891EiD+atSwJY/PPy/ddpvN2gEAohN1oZs+fbpOO+203T/u27evJOlPf/qTBg8erJUrV2rp0qW7//22bdt0xx13aMWKFapUqZIaN26sjz/+WOeee24M4icnDhRGqunRQ3rrLWn0aOmss1ynAQD/KdOmiEQJh8MKBoMpsSlixw7bnHnvvdJdd7lO41MzZ0otW0ozZnDui09EItIf/iA1aCAd4PIZAMBe8HDDY+bOlTZtYkMEUksgYIcLf/SR9KsJfgDAQaLQecyUKXYtUqtWrpP42PHHS1u32gjfuOoq6ZBDuN8VAEqDQucxU6fao6fKlV0n8bG0NCk9ndX1PnPIIbbj9aWXrI8DAA4ef+N5TFGRdOKJrlP43LffSu3b2whfuflmO4/u3XddJwEAf6HQeUg4LM2bJ51wguskPrdxozR+vI3wlWOPlc48U8rPd50EAPyFQuch06fbbj9m6JDKevWypQczZrhOAgD+QaHzkKIiKTPTZimAVHX++VLduszSAUA0KHQeUlRku1vLlXOdBHCnXDk7aHjoUGntWtdpAMAfPF3o8vPzlZubq7y8PNdREmLqVNbPxUTdurZVsm5d10lQSt27SyUl0quvuk4CAP7ATREesWKFVKeO9P770kUXuU4DuNelizRpkrRgAbPWAHAgnp6hSyVFRTYyQxcDa9ZIL79sI3yrVy9pyRJp5EjXSQDA+yh0HlFUJNWqJdWu7TpJEli6VLr+eu6Q8rkTTrAredkcAQAHRqHziKlTOa4E+LVAQOrZ02boFi1ynQYAvI1C5wE7d9oZdDxuBfZ0+eVS1arSc8+5TgIA3kah84D586UNGyh0wG9VqiR17Sq98oq0ebPrNADgXRQ6DygqssdLrVq5TpIkDjlEOvVUG+F7N90khUJ2Lh0AYO84tsQDbrpJKiiQ5s51nQTwpvPPl374wa4DCwRcpwEA72GGzgOKinjcGlMlJVJxsY1ICj17SrNmSYWFrpMAgDdR6BzbskWaM4dCF1OzZ0sZGTYiKZx9tnT00dLTT7tOAgDeRKFzbPZsaccOCh2wP2lpdtDwu+/arSoAgD1R6BybOlVKT5f+8AfXSQBvu+46m3h9/nnXSQDAeyh0jhUVSS1aSBUquE4CeFtWlpW6F16Qtm51nQYAvIVC5xgbIoCD17u3tHo1R5gAwG95utDl5+crNzdXeXl5rqPExdq1dqURhS7GmjSRli2zEUmlYUPp3HOlQYMk7x+4BACJwzl0Dn36qf3ltHCh1KCB6zSAP4waZbtex4+X2rVznQYAvMHTM3TJrqhIqlZNOuoo10mSzOLF0h//aCOSzllnSY0aSU895ToJAHgHhc6hXevnOPk+xtavl/79bxuRdAIBW0s3fLj0/feu0wCAN1DoHIlE2BABlNY110iZmVJ+vuskAOANFDpHliyR1qyRTjzRdRLAf6pUkbp3l156Sdq0yXUaAHCPQudIUZGNSbqBF4i7Xr2kcFgaMsR1EgBwj0LnSFGRVL++VKOG6yRJqFYt6aGHbETSOvJI6YILOMIEACQKnTNFRTxujZsjjpD69bMRSe2WW6Svv5ZGj3adBADcotA5sH27NGMGGyLiZv16acQIdrmmgFNPtXuQOcIEQKqj0Dnw1Vd2FyWFLk4WL7ZncZxDl/QCAZul+/hjO6AbAFIVhc6BoiKpXDmpeXPXSQD/u/JKO6D76addJwEAdyh0DhQV2WOiypVdJwH8r1Il6cYbpVdflUIh12kAwA0KnQNTp/K4FYilnj2l4mI7lw4AUpGnC11+fr5yc3OVl0SHtW3YYLvyKHRxlJEh5ebaiJRQq5Y9eh00yDYdAUCqCUQi3j/BKRwOKxgMKhQKKSsry3WcMhk3TjrtNOnLL6UmTVynAZLHF19Ixx8vvfWWdMUVrtMAQGJ5eoYuGU2dKh1yiNSokeskQHJp1kw680zpiSc4aBhA6qHQJVhRkdSqle1yRZzMni1lZdmIlNK3rzR9ujRhguskAJBYFLoEKypi/VzclZTYYsWSEtdJkGDnnGPLJ594wnUSAEgsCl0C/fCDtHw5hQ6Il0BAuu02uyhkwQLXaQAgcSh0CVRUZCN3uALxc/XVUvXq0sCBrpMAQOJQ6BKoqEjKzpZq13adBEheGRl2Lt1rr0nr1rlOAwCJQaFLoF3r5wIB10mS3HHHSTNm2IiUdPPNtoTy+eddJwGAxKDQJUhJiTRtGuvnEqJyZalFC+5WS2E1akjXXGP3uxYXu04DAPFHoUuQ+fOlcJj1cwmxdKk9c1u61HUSOHTbbdKqVdLbb7tOAgDxR6FLkF0bIlq1cpsjJaxZIz37rI1IWY0aSeeey0HDAFIDhS5BiopsSVcw6DoJkDpuv12aM0f6/HPXSQAgvih0CVJUxONWINFOO01q2VJ65BHXSQAgvih0CbB1q10czoYIILECAemuu6TRo23jMwAkK08Xuvz8fOXm5iovL891lDKZPVvavp1ClzA1a9qK+Jo1XSeBB1x8sXTUUczSAUhugUjE+8uFw+GwgsGgQqGQsrKyXMeJ2qBBNksQDksVK7pOA6Se556TevWSvv1WatDAdRoAiD1Pz9Ali6IiqXlzylzCbNwoTZliIyDp2mulww6THn/cdRIAiA8KXQJMncrj1oT69lvppJNsBCRVqiTdcotdB/bTT67TAEDsUejibN06aeFCCh3g2s03S+XK2e0RAJBsKHRxNm2ajRQ6wK1q1aTrr5fy83kaDyD5UOjirKhIqlpVOvpo10kA9O0rbdggvfyy6yQAEFsUujibOlXKy7PzsJAg5ctL1avbCPxKTo505ZV2Hdj27a7TAEDsUOjiKBLhhggn/vAHafVqG4HfuPNOadkyaehQ10kAIHYodHH0/ffWK1g/B3hHkybS+edLDz8slZS4TgMAsUGhi6OiIhspdAk2d64tWpw713USeFT//tK8edIHH7hOAgCxkZBCV1BQoE6dOqlWrVoKBAIaPnx4It7WuaIi6cgjuYEq4YqLpUWLbAT2onVr6bTTpIcesqURAOB3CSl0mzZtUrNmzfTMM88k4u08o6iI2TnAq/r3l2bOlD77zHUSACi7hGwD7Nixozp27JiIt/KMHTukGTOkBx5wnQTA3px+um1Y+sc/pHPOcZ0GAMrGk2voiouLFQ6H93j5zdy50ubNzNABXhUI2CzdxIlSQYHrNABQNp4sdAMGDFAwGNz9ysnJcR0pakVFds1Qixauk6Sgo4+WRo7kNGcc0HnnSU2b2iwdAPiZJwtdv379FAqFdr+WLVvmOlLUiorseIQqVVwnSUFZWdLZZ9sI7EdamnT33dKoUb9c0wcAfuTJQpeenq6srKw9Xn7DhgiHVq6U7rvPRuAA/vhHm8x96CHXSQCg9DxZ6Pxu0ybpq68odM6sXCndfz+FDgelXDnpL3+Rhg/n6EIA/pWQQrdx40bNnj1bs2fPliQtWbJEs2fP1tKlSxPx9gk3Y4adQE+hA/yhSxe753XAANdJAKB0ElLopk+frubNm6t58+aSpL59+6p58+a69957E/H2CVdUJFWuLOXmuk4C4GBUrCjddZfd7/rtt67TAED0ElLo2rdvr0gk8rvX4MGDE/H2CVdUJLVsKZVPyCl/AGKhe3fpiCOkBx90nQQAoscaujgoKrIDS+FI1arSVVfZCBykjAypXz/pzTel+fNdpwGA6AQiEe/fZBgOhxUMBhUKhTy/4/XHH+27/GHDbPccAP/YutV2vJ52mvTGG67TAMDBY4YuxnadZcWGCIe2bpUWLrQRiMKuWbq33mKWDoC/UOhibOpUqWZNqW5d10lS2NdfSw0b2ghEqVs3KTtb+vvfXScBgINHoYuxXQcKBwKukwAojYwMuz1i6FDpm29cpwGAg0Ohi6FIhA0RQDLo1k2qVYtZOgD+QaGLoYULpfXrWT8H+F16OrN0APyFQhdDRUU2tmrlNgeAsuvaVapTh1k6AP7g6UKXn5+v3Nxc5eXluY5yUIqKbC1+tWquk6S4Fi3s+XeLFq6TwMd+PUs3b57rNACwf5xDF0Nt2kgNGkhDhrhOAiAWiovtm7RTTrGjTADAqzw9Q+cn27ZJs2axfs4T5s+3ds1BYiij9HSpf3/p7beZpQPgbRS6GJkzx76bZ4erB2zaJBUW2giU0XXXSTk50n33uU4CAPtGoYuRoiKpQgWpWTPXSQDEUsWK0r332nV+M2e6TgMAe0ehi5GiIitzGRmukwCItT/9STr2WHv8CgBeRKGLkV03RABIPuXL2/ElI0dKBQWu0wDA71HoYiAUssNHKXQeceSR0htv2AjEyMUX20k4/frZqTgA4CUUuhiYPt3+gKfQeUS1atLVV3MgIGIqLU166CFp8mTp449dpwGAPVHoYmDqVCkYtDU28IDVq6X8fBuBGOrQQTr1VFtLV1LiOg0A/IJCFwOFhXZcSRr/Nb1h2TKpVy8bgRgKBKQBA+yYorffdp0GAH5BBSmjSOSXQgcg+bVpI3XqZEeZbN/uOg0AGApdGS1ZYk/2Wrd2nQRAovzjH9LixdIrr7hOAgCGQldGhYU2MkMHpI6mTaUrr5QeeEDavNl1GgDweKHLz89Xbm6u8vLyXEfZp6lT7fLuww5znQS7ZWba6vXMTNdJkMQeeEBas0YaONB1EgCQApGI909UCofDCgaDCoVCysrKch1nDyeeaLtbX3/ddRIAidanj/Tqq9LChVLNmq7TAEhlnp6h87qtW6VZs3jc6jk7d0rhsI1AHP31r7a7/YEHXCcBkOoodGUwa5btcmNDhMd88YUdDPjFF66TIMkddpidSff889L8+a7TAEhlFLoyKCyUMjKkP/zBdRIArvTuLdWpI/3lL66TAEhlFLoyKCyUWrWSKlRwnQSAKxkZdiXY8OHShAmu0wBIVRS6Mpg6lcetAKTLL7dv7m6/nSvBALhBoSullSul779nQwQA2xjx2GPStGnSsGGu0wBIRRS6Upo61UZm6DyoaVPpp59sBBLk1FOlzp2lfv2k4mLXaQCkGgpdKRUWSrVr22JoeEyFClKNGixuRML985/SsmXSoEGukwBINRS6UiosZHbOsxYtsqmSRYtcJ0GKOe446eabpb//XVq1ynUaAKmEQlcKO3ZI06dT6DwrFJI++shGIMHuv1+qWNEevQJAolDoSmHuXGnTJjZEAPi9qlWlBx+UBg+WiopcpwGQKih0pVBYKJUrJ7Vs6ToJAC+6/nqpWTPplls4xgRAYlDoSqGw0P6wrlzZdRIAXlSunG2MmDpVeuMN12kApAJPF7r8/Hzl5uYqLy/PdZQ9sCHC42rXlh5/3EbAkXbtpMsusyvBNmxwnQZAsgtEIpGI6xAHEg6HFQwGFQqFlJWV5TTL+vW2Rub116UuXZxGAeBxy5ZJxx5r973+85+u0wBIZp6eofOiXYuc2RDhYT//LL37ro2AQzk5NkP35JPSggWu0wBIZhS6KBUW2gxdw4auk2CfliyRLr3URsCxO++UatWSbrvNdRIAyYxCF6Vd6+cCAddJAPhBpUp2z+vHH9sLAOKBQheFkhIrdDxuBRCNiy+WzjpL6tVL2rzZdRoAyYhCF4X5821Z1sknu04CwE8CASk/X1q50g4dBoBYo9BFYfJkKS2NGTrPq1RJat7cRsAjGjaU7r5bevRRu20GAGKJY0ui0LWrNGuWvQAgWsXFUtOmUna2NG4ca3EBxA4zdFGYPFk66STXKQD4VXq69OyzUkGBnWUJALFCoTtIa9bYGjrWz/nArFn2NydTqfCgM8+UrrxSuuMOae1a12kAJAsK3UGaMsVGZuh8IBKRtm2zEfCgxx+Xtm+3Q4cBIBYodAdp8mRb91KvnuskAPzuiCOkhx6SXn5ZmjTJdRoAyYBCd5AmTbLHrSxiBhALN95oO+a7d5e2bnWdBoDfUegOwrZt0rRpPG4FEDvlykmvvCItWsTZdADKztOFLj8/X7m5ucrLy3OaY9Ys+w6aDRE+0aiR9NVXNgIe1rixdM890sMPS7Nnu04DwM84h+4gPPmkHQgaCkkVKyb87QEksW3bpFatpPLlpaIiGwEgWp6eofOKyZOlE06gzPnG99/bwqTvv3edBDigihXt0esXX0iPPeY6DQC/otAdQCRiGyJYP+cja9fa35Ac8gWfyMuT+vaV7rvPzrsEgGhR6A7g++/tQm0KHYB4uv9+KSdH6tZNKilxnQaA31DoDmDXGVFt2rjNASC5Va78y7l0Tz3lOg0Av6HQHcDkydKxx0rVq7tOAiDZnXqq1KeP1K+f9PXXrtMA8BMK3QFMnsxxJb5z+OF2p9Lhh7tOAkTtoYek+vWla66x68EA4GBQ6PZjwwZpzhwet/pO7drSgAE2Aj5TqZL0+ut2Lt0//uE6DQC/oNDtx+TJtji5bVvXSRCVDRukceNsBHwoL88OHH7wQWn6dNdpAPgBhW4/CgqkmjWlY45xnQRRWbBAOu00GwGf6t9fOv54e/S6ZYvrNAC8rlSF7tlnn1X9+vWVkZGhli1basKECfv82MGDBysQCPzutdUHt1FPmCC1aycFAq6TAEg1FSrYo9fFi+2mGgDYn6gL3TvvvKM+ffqof//+mjVrltq2bauOHTtq6dKl+/ycrKwsrVy5co9XRkZGmYLH29at0tSpPG4F4E5urvTPf0oDB0qffuo6DQAvi7rQPfHEE+rWrZu6d++uRo0aaeDAgcrJydFzzz23z88JBAI64ogj9nh5XVGR3bHYrp3rJABS2S23SOeeK/3pT3bIOQDsTVSFbtu2bZoxY4Y6dOiwxz/v0KGDJk+evM/P27hxo+rVq6c6dero/PPP16xZs/b7PsXFxQqHw3u8Em3CBCkYlJo2Tfhbo6wqVLAdrhUquE4ClFkgIL32mlSunK2n4xYJAHsTVaFbs2aNdu7cqcN/c77X4YcfrlWrVu31c4477jgNHjxYI0aM0NChQ5WRkaGTTz5ZC/azYH3AgAEKBoO7Xzk5OdHEjImCAjt/rly5hL81yqppU2n5cto4kkbNmtIbb0ijR0uPPeY6DQAvKtWmiMBvdglEIpHf/bNdWrdurauvvlrNmjVT27ZtNWzYMB1zzDF6+umn9/nz9+vXT6FQaPdr2bJlpYlZajt22PU7PG4F4BVnninddZftfi0qcp0GgNdEVeiqV6+ucuXK/W427qeffvrdrN0+3zAtTXl5efudoUtPT1dWVtYer0SaNUvatIlC51tffinVqWMjkET+/nepRQvpiiskBytRAHhYVIWuYsWKatmypT7//PM9/vnnn3+uk0466aB+jkgkotmzZys7Ozuat06oCRPstPaWLV0nQals3y6tWMG9SUg6FSpIb70lrV4t9eghRSKuEwHwivLRfkLfvn3VpUsXtWrVSm3atNGLL76opUuXqkePHpKka665RrVr19aAAQMkSffff79at26thg0bKhwOa9CgQZo9e7by8/Nj+yuJoYICqXVrqWJF10kAYE8NGkgvvCBdeaV0yinSzTe7TgTAC6IudJdddpnWrl2rBx54QCtXrlSTJk30ySefqF69epKkpUuXKi3tl4m/9evX64YbbtCqVasUDAbVvHlzFRQU6IQTTojdryKGSkpshq53b9dJAGDvrrjCribs08cewbZu7ToRANcCkYj3J+3D4bCCwaBCoVDc19PNnSs1aSL997/SGWfE9a0QLzNn2vPyGTPsbzsgCW3bJrVvLy1bZr/la9RwnQiAS9zl+hsFBVL58nzH62sNG0pjx9oIJKmKFaVhw6TiYpux27nTdSIALlHofqOgQGrVSqpSxXUSlFpmpk1dZGa6TgLEVZ060ttv2/cv997rOg0Alyh0vxKJSOPGcVyJ761YIfXrZyOQ5E4/XXroIXuNGOE6DQBXKHS/8s030qpV9gckfOzHH6WHH7YRSAF33SVdeKFdDfbtt67TAHCBQvcrY8bY+rlTTnGdBAAOXiAgDR4sZWdLnTtL69e7TgQg0Sh0vzJ2rG2GYP0cAL8JBu2R608/SZddZlcYAkgdFLr/KSmxQsfjVgB+1bCh9O9/S6NHS7ff7joNgETydKHLz89Xbm6u8vLy4v5ec+ZI69ZR6JLCYYdJ3brZCKSY00+Xnn5aGjRIevFF12kAJAoHC//PE09I/fvb2pP09Li8BQAkTM+eVuhGjZJOO811GgDx5ukZukQaO1Y6+WTKXFLYssWu/NiyxXUSwJmBA6VTT5X+7/+kr792nQZAvFHoZIuHx4/ncWvSmDfP7m+bN891EsCZChVsPV2dOlLHjtIPP7hOBCCeKHSyKz83bKDQAUguhx4qffqpXQt27rlSOOw6EYB4odDJzp/LzLQrvwAgmdSpY6VuyRLpkkukbdtcJwIQDxQ62fq5tm3tUGEASDZNm0rDh9vSkuuvt2sOASSXlC90xcXSxInsAksqgYBUsaKNACTZn3H/+pf0+uvSn/9MqQOSTcrPSU2caJshO3RwnQQx07y5NXUAe7j8crviuE8fW193992uEwGIlZQvdJ99ZvcfNm3qOgkAxN+tt9rmiP79be1w796uEwGIhZR/5PrZZzY7x9O5JDJvntSiBceWAPtwzz3SHXdIt9wiDR7sOg2AWEjpGbqVK+3Krz//2XUSxNSWLdKsWRwsDOxDICA98ojN1HXrJlWpIv3xj65TASiLlC50o0bZH2xnneU6CQAkViAgPfustGmTdMUV9s8odYB/pXSh++wzezJXo4brJACQeOXK2SPXSMRKXUmJdNllrlMBKI2ULXQlJdLnn0s33OA6CQC4U768HWVSrpx05ZV2q8SVV7pOBSBani50+fn5ys/P186dO2P+c8+cKa1ZI519dsx/arhWv740bJiNAA6oXDnptdds7NLFSl2XLq5TAYhGIBLx/vGS4XBYwWBQoVBIWVlZMfk5//EP6eGHpXXr7BJrAEh1JSX21OLVV6WXX5a6dnWdCMDB8vQMXTyNGiWdfjplLin9+KP05pvSVVdJhx/uOg3gG2lp0osv2p+L3bpJa9dKd97pOhWAg5GS59CFQtLkydI557hOgrhYsUK6/XYbAUQlLc12v/bvL911l728/xwHQErO0I0cKe3YIZ13nuskAOA9gYD04IN2AkCfPrbe+MUXbQMFAG9Kyf97fvSR1KyZVLeu6yQA4F233ioddph03XVW6oYOtUOIAXhPyj1y3bFD+uQTqVMn10kAwPuuvloaMUIaM0Zq146VDIBXpVyhmzxZ+vlnCl1SCwbtCxwMuk4CJIWOHaVJk6SffpJOOMFu1gPgLSlX6D76SDriCKlVK9dJEDcNGtiUQoMGrpMASaNZM6moSKpVSzrlFPu/GADvSMlCd955tpMLSWr7dmn1ahsBxEx2tjR+vM3YXXih9Pjj7IAFvCKlas2CBdL8+TxuTXpffinVrGkjgJiqXNkuYvnLX6Q77pCuuUbavNl1KgApVeg++khKT5fOPNN1EgDwr7Q06aGHpLfekt5/X2rTRlq0yHUqILWlVKF77z2pQwe23QNALFxxhVRYaDN0LVvaN80A3EiZQrdihe1wveQS10kAIHk0bSpNny61by917mw3S2zb5joVkHpSptB98IHdT8j6OQCIrWDQHr0+8oj05JPSySdLCxe6TgWkFk8Xuvz8fOXm5iovL6/MP9e//21r56pWjUEweFuzZnZhb7NmrpMAKSMtTbrzTmnKFGn9eql5c+mNN1ynAlJHIBLx/qbzcDisYDCoUCikrKysqD//xx9tu/3LL0tdu8YhIABgtw0bpN69pX/9S7rySunpp6Vq1VynApKbp2foYuWDD+y7xwsucJ0ECbFggXT22TYCSLjMTGnwYOnNN+2qxdxcafhw16mA5JYShW7oUOmMM+ySaaSADRukUaNsBODMlVdKc+dKJ54oXXSR7Ypds8Z1KiA5JX2h+/57qaDALpgGACRWrVo2O/fmm/Z9Vm6u9O673DABxFrSF7q33rKTzS+6yHUSAEhNgcAvs3Vt20qXXiqdey47YYFYSupCF4nYLqsLL5QOOcR1GgBIbUccYQe8f/ihNG+e1KSJ9Le/SVu2uE4G+F9SF7pZs+wPDR63ppicHOmZZ2wE4DmdO0tff213wT78sD2GfecdHsMCZZHUhW7IELuj/ayzXCdBQtWoIfXsaSMAT6pcWXrwQenLL+22icsvtzthJ050nQzwp6QtdNu22SLcK66Qypd3nQYJtW6dtfl161wnAXAAxxwjjRghjRkjbd9ua+wuvphTh4BoJW2hGz5c+ukn6frrXSdBwn33ndSli40AfOG006Rp02zd8/TpUqNG0rXXUuyAg5W0he755+07vcaNXScBAByMtDRb8zx/vt0JO2qUdNxx0jXXSN9+6zod4G1JWei++UYaO1bq0cN1EgBAtDIy7OqwxYulp56SRo+2GbsuXaQ5c1ynA7wpKQvdiy9K1avbOgwAgD9lZEi9ekmLFkmDBknjx0vNmtlGt08/lUpKXCcEvCPpCt3GjXaH4LXXSunprtPAiSpVpNatbQTgexkZtnF90SK7ynH9ejuYuHFj6YUX7M99INV5utDl5+crNzdXeXl5B/05r74qhcP2XR1S1LHHSlOm2AggaVSoYMebFBVJEybYY9ibbpKys6Ubb5RmzHCdEHAnEIl4/yjHcDisYDCoUCikrKysfX7c9u3S0UfbZoghQxIYEADgxAtJrZMAAApbSURBVPff2zfyr7wirVghNW9upxtcfrlUtarrdEDieHqGLlrvvistXSrdeafrJHBq5ky7PHLmTNdJAMRZvXrS/ffbKUUffWQXxPTqZdeMXXSR/b3A1WJIBUlT6HbssFPHzz7bFs0CAFJH+fLS+efbPbHLl0v//KeNl15q5e7qq+0eWdbbIVklTaF7/XW7t/XBB10nAQC4lJ0t9eljBxV/841022123Mkll9gJCJ0722Pa1atdJwViJynW0G3ZYtfHtGkjDRvmICC8ZeZMqWVLWyHdooXrNAA8YuFCm8H74ANp8mT7Zy1aSB062Oukk6SKFd1mBEorKWboBg2SVq5kdg4AsG9HHy3dfrs0caL9nfHqq7YZ/qWX7OqxatXsse0jj9hG+W3bXCcGDp7vZ+gWLZKaNrUt608+6SggvGXrVls8U6eOHWAFAPtRUiLNnm1Xjf33v1bmNm+2Pz5OPFE65RSbvWvVSqpZ03VaYO98XegiETsxfOFC6auvpEMOcRgSAJAUtm+3gjdhgs3mTZz4y3q7unWt2LVqZUekNG5s3zsGAm4zA74udE89ZQtfR4603a2AJGnJEumvf5X+/nepfn3XaQD4XCRif6xMn/7La8YMO8RekrKyrNg1afLLeOyxUq1aUlpSLGyCH/i20I0da4tYb71VeuwxxwHhLWyKABBnJSV29t3cufb66isb582TiovtY9LT7XvKBg2ko4765dWggc30ZWY6/SUgyfiy0E2ZYjNyrVtLH39s18EAu1HoADiyY4et7V6wQFq82F6LFv3yv7du/eVjMzNtFq92bRt3/e8jjrDjVQ477JdXpUrufk3wh/KuA0SjpMSud+nZU8rLk95/nzIHAPCO8uXtceverpIuKZF+/NEK3rJl0g8/2HVlP/xgs32TJ9uPd83w/VqlSr8vedWqWSnc36tyZdvcUamSjRkZlhHJx9mXNRKJaMOGDXv9d8XFxSr+1e/ooUPt3pbWrcOaP99O/H78cfs/x641DMBuu46C37iR3yAAPKVKFekPf7DX3kQiUigkrVtnr59/tnHt2l/+2bp10qpVdmjyxo3Shg02HuwVZ+XK/VLufv2qWNH+XfnypXsFAvZKS/vlf+/tn0kH/ph9/TwulfX9zzvPNtAcrMzMTAWieFNnj1x3PUYFAADAnvZ1mcK+OCt00czQLV++Um3anKCvv/5atWvXTki+vLw8TZs2LSHvlezvFw6HlZOTo2XLlkX1m7Mskvm/p4v342vo7/fj6+fv9+Pr5//3K83XMNoZOmePXAOBQNS/MTMzMxP2m7lcuXIJe69UeD9JysrK4uvn0/fbha+hf99P4uvn5/eT+Pr5+f12iefXkBNy9qFnz568n48l+3/PZP/6Scn/3zTZv4bJ/t+Trx/v5zW+OLZk+fLlu6cq60SzohCesL+r2+APfA39ja+fv/H1879EfA19MUOXnp6+xwh/SU9P19/+9je+fj7G19Df+Pr5G18//0vE19AXM3R8dwIAALBvvpihAwAAwL5R6AAAAHyOQgcAAOBzvlhDt+sQ4mgP2QMAAEgFvih0AAAA2DceuSImnn32WdWvX18ZGRlq2bKlJkyYsM+Pfemll9S2bVtVrVpVVatW1ZlnnqmioqIEpsVvRfP1+7W3335bgUBAF154YZwT4kCi/RquX79ePXv2VHZ2tjIyMtSoUSN98sknCUqL34r26zdw4EAde+yxqlSpknJycnTbbbdp69atCUqLg1VQUKBOnTqpVq1aCgQCGj58eNzei0KHMnvnnXfUp08f9e/fX7NmzVLbtm3VsWNHLV26dK8fP27cOF1xxRUaO3aspkyZ8v/t3V9IU30cBvCHpktQb0QdQ1RwNwamoJKUk1ByN1HdRVDiRSLqAiW6mIgt/AdleKUJ2i66CAuEIEJCI+0i72QDQZnamIF/sKCLoRfz6Pe9qb3v9PXPmZ0zDz4fOBf+ZNtzeDj63TjnDDk5OXA4HFhZWdE5OQHq+/tjeXkZjx49QkVFhU5J6SBqOwyHw6iurkYwGMTo6Cj8fj+Gh4d1+65siqa2v9evX8PlcsHtdmN+fh4ejwdv375Fa2urzsnpKJubmygqKkJ/f7/2LyZEJ3Tp0iVpaGiIWsvPzxeXy3WsxyuKIqmpqfLq1Sst4tERYulPURQpLy+Xly9fSm1trdy6dUvrmHQItR0ODg5KXl6ehMNhPeLREdT253Q6paqqKmrt4cOHYrfbNctIJwdA3r17p9nz8xM6OpFwOIyZmRk4HI6odYfDgenp6WM9x9bWFra3t5GWlqZFRDpErP11dHQgIyMD9+/f1zoiHSGWDt+/f4/Lly/D6XTCYrGgoKAAPT092NnZ0SMy/Ucs/dntdszMzEROVQkEAhgbG8P169c1z0unV0K8A5Cx/fz5Ezs7O7BYLFHrFosF6+vrx3oOl8uFrKwsXLt2TYuIdIhY+vv69Ss8Hg98Pp8eEekIsXQYCATw+fNn3L17F2NjY1hcXITT6YSiKHj8+LEesem3WPq7c+cOfvz4AbvdDhGBoihobGyEy+XSIzKdUhzo6K/YezsZETnWLWaePXuGkZERTE1NISkpSat4dITj9hcKhXDv3j0MDw8jPT1dr3h0DGqOwd3dXWRmZmJoaAgmkwklJSVYXV1Fb28vB7o4UdPf1NQUuru78eLFC5SVlWFpaQnNzc2wWq1ob2/XIy6dQhzo6ETS09NhMpn2vZPc2NjY945zr+fPn6OnpwefPn1CYWGhljHpAGr7+/btG4LBIG7cuBFZ293dBQAkJCTA7/fDZrNpG5qixHIMWq1WJCYmwmQyRdYuXLiA9fV1hMNhmM1mTTPTv2Lpr729HTU1NairqwMAXLx4EZubm6ivr0dbWxvOnePZVGcRW6cTMZvNKCkpwcTERNT6xMQErly5cuDjent70dnZiY8fP6K0tFTrmHQAtf3l5+djdnYWPp8vst28eROVlZXw+XzIzs7WKzr9FssxWF5ejqWlpcgwDgALCwuwWq0c5nQWS39bW1v7hjaTyQQRgfDWsmeXZpdb0Jnx5s0bSUxMFI/HI3Nzc9LS0iLJyckSDAZFRKSmpibqaq2nT5+K2WyW0dFRWVtbi2yhUCheu3Cmqe1vL17lGn9qO/z+/bukpKTIgwcPxO/3y4cPHyQzM1O6urritQtnmtr+3G63pKamysjIiAQCARkfHxebzSa3b9+O1y7QAUKhkHi9XvF6vQJA+vr6xOv1yvLy8l9/LQ509FcMDAxIbm6umM1mKS4uli9fvkR+d/XqVamtrY38nJubKwD2bW63W//gJCLq+tuLA93poLbD6elpKSsrk/Pnz0teXp50d3eLoig6p6Y/1PS3vb0tT548EZvNJklJSZKdnS1NTU3y69evOCSnw0xOTv7v/7vD/qbGil/9RURERGRwPIeOiIiIyOA40BEREREZHAc6IiIiIoPjQEdERERkcBzoiIiIiAyOAx0RERGRwXGgIyIiIjI4DnREREREBseBjoiIiMjgONARERERGRwHOiIiIiKD+wf+WU/cI9sIPQAAAABJRU5ErkJggg==\n",
"text/plain": [
"Graphics object consisting of 2 graphics primitives"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_ex(x) = f(x).subs(a==3).subs(b==7);\n",
"m_ex = m(3,7);\n",
"plot(f_ex(x),(x,0,1))+line([(m_ex,0),(m_ex,f_ex(m_ex))], color='red',linestyle='--')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We define $f^{-1}(x)$ to be the inverse of left and right side of $f(x)$ for $1>x>m(a,b)$ and $m(a,b)>x>0$ respectively."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"finv=function('finv', latex_name=\"f^{-1}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $r(x)=f(f^{-1}(x))$, which reflects a point $x$ to throught the curve of $f(x)$. We want to show that $r(x)$ is convex if $b>a>1$. See my [previous post](./the-convexity-property-related-to-beta-densities.html) for pictures."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"r(x)=finv(f(x))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"r(x).show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is $f(x)$'s direvative."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"((a - 1)*x^(a - 2)*(-x + 1)^(b - 1) - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2))/beta(a, b)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df(x)=diff(f(x),x).combine();\n",
"df(x).show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is $f^{-1}(x)$'s direvative."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"-beta(a, b)/((b - 1)*(-finv(x) + 1)^(b - 2)*finv(x)^(a - 1) - (a - 1)*(-finv(x) + 1)^(b - 1)*finv(x)^(a - 2))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"dfinv(x)=1/df(finv(x));\n",
"dfinv(x).show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## First direvative of $r(x)$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/beta(a, b) - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/beta(a, b))*D[0](finv)(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd1=diff(r(x),x);\n",
"rd1.show()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"D[0](finv)(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"to_be_replace=rd1.operands()[1];\n",
"to_be_replace.show()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"-((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/beta(a, b) - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/beta(a, b))*beta(a, b)/((b - 1)*(-finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) + 1)^(b - 2)*finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))^(a - 1) - (a - 1)*(-finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) + 1)^(b - 1)*finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))^(a - 2))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd1_0=rd1.subs(to_be_replace==dfinv(f(x)));\n",
"rd1_0.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make it a easier to manipulate by changing some symbobls. Let $B=\\beta(a,b)$ and $x_1$ be the reflection point of $x$."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/B - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/B)*B/((a - 1)*x1^(a - 2)*(-x1 + 1)^(b - 1) - (b - 1)*x1^(a - 1)*(-x1 + 1)^(b - 2))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd1_1=rd1_0.subs(finv(f(x))==x1).subs(beta(a,b)==B);\n",
"rd1_1.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Simplify using sympy."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"rd1_2=rd1_1._sympy_().simplify()._sage_();"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"((a - 1)*x^(a - 2)*(-x + 1)^(b - 1) - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2))/((a - 1)*x1^(a - 2)*(-x1 + 1)^(b - 1) - (b - 1)*x1^(a - 1)*(-x1 + 1)^(b - 2))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd1_2.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we have three equalities."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"x^(a - 1)*(-x + 1)^(b - 1) == x1^(a - 1)*(-x1 + 1)^(b - 1)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"eq1=beta(a,b)*f(x)==beta(a,b)*f(x1);\n",
"eq1.show()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"x1^(a - 1) == x^(a - 1)*(-x + 1)^(b - 1)/(-x1 + 1)^(b - 1)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"eq2=eq1.op[1]/eq1.op[1,1]==eq1.op[0]/eq1.op[1,1];\n",
"eq2.show()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(-x1 + 1)^(b - 1) == x^(a - 1)*(-x + 1)^(b - 1)/x1^(a - 1)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"eq3=eq1.op[1]/eq1.op[1,0]==eq1.op[0]/eq1.op[1,0];\n",
"eq3.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we simplify again use these equations."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a*(x - 1) + b*x - 2*x + 1)*(x1 - 1)*x1/((a*(x1 - 1) + b*x1 - 2*x1 + 1)*(x - 1)*x)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd1_3=rd1_2._sympy_().subs(eq2.op[0],eq2.op[1]).subs(eq3.op[0],eq3.op[1]).simplify().together()._sage_()\n",
"rd1_3.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The final form for $r'(x)$ is"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a*(x - 1) + b*x - 2*x + 1)*(finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) - 1)*finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))/((a*(finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) - 1) + b*finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) - 2*finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) + 1)*(x - 1)*x)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd1_4=rd1_3.subs(x1==finv(f(x)));\n",
"rd1_4.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Second direvative"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"a*(x - 1) + b*x - 2*x + 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"1/(a*(finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) - 1) + b*finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) - 2*finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) + 1)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/beta(a, b) - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/beta(a, b)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"1/(x - 1)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"1/x"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"finv(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b)) - 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"D[0](finv)(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd2_0=rd1_4.diff(x);\n",
"[o.show() for o in rd2_0.op[0].op];"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"D[0](finv)(x^(a - 1)*(-x + 1)^(b - 1)/beta(a, b))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"to_be_replace_2=rd2_0.op[0].op[-1];\n",
"to_be_replace_2.show();"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a*(x - 1) + b*x - 2*x + 1)*((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/B - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/B)*B*(x1 - 1)/(((a - 1)*x1^(a - 2)*(-x1 + 1)^(b - 1) - (b - 1)*x1^(a - 1)*(-x1 + 1)^(b - 2))*(a*(x1 - 1) + b*x1 - 2*x1 + 1)*(x - 1)*x) + (a*(x - 1) + b*x - 2*x + 1)*((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/B - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/B)*B*x1/(((a - 1)*x1^(a - 2)*(-x1 + 1)^(b - 1) - (b - 1)*x1^(a - 1)*(-x1 + 1)^(b - 2))*(a*(x1 - 1) + b*x1 - 2*x1 + 1)*(x - 1)*x) - (((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/B - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/B)*B*a/((a - 1)*x1^(a - 2)*(-x1 + 1)^(b - 1) - (b - 1)*x1^(a - 1)*(-x1 + 1)^(b - 2)) + ((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/B - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/B)*B*b/((a - 1)*x1^(a - 2)*(-x1 + 1)^(b - 1) - (b - 1)*x1^(a - 1)*(-x1 + 1)^(b - 2)) - 2*((a - 1)*x^(a - 2)*(-x + 1)^(b - 1)/B - (b - 1)*x^(a - 1)*(-x + 1)^(b - 2)/B)*B/((a - 1)*x1^(a - 2)*(-x1 + 1)^(b - 1) - (b - 1)*x1^(a - 1)*(-x1 + 1)^(b - 2)))*(a*(x - 1) + b*x - 2*x + 1)*(x1 - 1)*x1/((a*(x1 - 1) + b*x1 - 2*x1 + 1)^2*(x - 1)*x) + (a + b - 2)*(x1 - 1)*x1/((a*(x1 - 1) + b*x1 - 2*x1 + 1)*(x - 1)*x) - (a*(x - 1) + b*x - 2*x + 1)*(x1 - 1)*x1/((a*(x1 - 1) + b*x1 - 2*x1 + 1)*(x - 1)*x^2) - (a*(x - 1) + b*x - 2*x + 1)*(x1 - 1)*x1/((a*(x1 - 1) + b*x1 - 2*x1 + 1)*(x - 1)^2*x)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd2_1=rd2_0.subs(to_be_replace_2==dfinv(f(x))).subs(finv(f(x))==x1).subs(beta(a,b)==B);\n",
"rd2_1.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a lot, but again we can use `eq2,eq3` to simplify."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"rd2_2=rd2_1._sympy_().subs(eq2.op[0],eq2.op[1]).subs(eq3.op[0],eq3.op[1]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some more simplification with `sympy`, and study the denominator and numerator separately."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"rd2_numer, rd2_denom=sympy.fraction(rd2_2.together())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The denominator."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a*(x1 - 1) + b*x1 - 2*x1 + 1)^2*(x - 1)^2*x^2*((b - 1)/(x1 - 1) + (a - 1)/x1)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd2_denom1=rd2_denom.simplify()._sage_();\n",
"rd2_denom1.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want the 2nd direvative to be positive when $b>a>1$, so we can ingore the first 3 terms, which are positive."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(b - 1)/(x1 - 1) + (a - 1)/x1"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd2_denom2=rd2_denom1.op[-1];\n",
"rd2_denom2.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The numerator."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"a^2*b*x^2 + a*b^2*x^2 - a^2*b*x1^2 - a*b^2*x1^2 - 2*a^2*b*x - a^2*x^2 - 4*a*b*x^2 - b^2*x^2 + 2*a^2*b*x1 + a^2*x1^2 + 4*a*b*x1^2 + b^2*x1^2 + 2*a^2*x + 4*a*b*x + 3*a*x^2 + 3*b*x^2 - 2*a^2*x1 - 4*a*b*x1 - 3*a*x1^2 - 3*b*x1^2 - 4*a*x - 2*b*x - 2*x^2 + 4*a*x1 + 2*b*x1 + 2*x1^2 + 2*x - 2*x1"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rd2_numer1=rd2_numer.simplify()._sage_();\n",
"rd2_numer1.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Put the numerator and denominator back and separte them again, which actually simplifies again.\n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"numer3,denom3=sympy.fraction((rd2_numer1/rd2_denom2)._sympy_().together())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The denominator."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a + b - 2)*x1 - a + 1"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"denom3.expand().collect(x1)._sage_().show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is $>0$ if $x_1 > m(a,b)= (a-1)/(a+b-2)$, which is our assumption. So the denominator does not matter."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Factor the numerator"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(a*x + b*x + a*x1 + b*x1 - 2*a - 2*x - 2*x1 + 2)*(a - 1)*(b - 1)*(x - x1)*(x1 - 1)*x1"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"numer4=numer3.factor()._sage_();\n",
"numer4.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we restrict that $x