Suppose you've done a (robust) Bayesian multiple linear regression, and now you want the posterior distribution on the predicted value of \(y\) for some probe value of \( \langle x_1,x_2,x_3, ... \rangle \). That is, not the posterior distribution on the
mean of the predicted value, but the posterior distribution on the predicted value itself. I showed how to do this for simple linear regression in
a previous post; in this post I show how to do it for multiple linear regression. (A lot of commenters and emailers have asked me to do this.)
The basic idea is simple: At each step in the MCMC chain, use the parameter values to randomly generate a simulated datum \(y\) at the probed value of \(x\). Then examine the resulting distribution of simulated \(y\) values; that
is the posterior distribution of the predicted \(y\) values.
To implement the idea, the first programming choice is whether to simulate the \(y\) value with JAGS (or Stan or whatever) while it is generating the MCMC chain, or to simulate the \(y\) value after the MCMC chain has previously been generated. There are pros and cons of each option. Generating the value by JAGS has the benefit of keeping the code that generates the \(y\) value close to the code that expresses the model, so there is less chance of mistakenly simulating data by a different model than is being fit to the data. On the other hand, this method requires us to pre-specify all the \(x\) values we want to probe. If you want to choose the probed \(x\) values after JAGS has already generated the MCMC chain, then you'll need to re-express the model outside of JAGS, in R, and run the risk of mistakenly expressing it differently (e.g., using precision instead of standard deviation, or thinking that
y=rt(...) in R will use the same syntax as
y~dt(...) in JAGS). I will show an implementation in which JAGS simulated the \(y\) values while generating the MCMC chain.
To illustrate, I'll use the example of scholastic aptitude test (SAT) scores from Chapter 18 of DBDA2E, illustrated below:
Running robust multiple linear regression yields a posterior distribution as shown below:
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA0gAAAHgCAIAAAAQVdomAAAgAElEQVR4nO3duXHbTByHYXSiDpRJAWtQ5tyjiBWoBM045bgKpUrYg5KvAZbgIvAFuPbGAtgb7zMOTBLHCtw/9kccZNcDAACgCV3uBgAAACAMgh0AAEAjCHYAAACNINgBAAA0gmAHAADQCIIdAABAIwh2AAAAjSDYAQAANIJgBwAA0AiCHQAAQCMIdgAAAI0g2AEAADSCYAcAANAIgh0AAEAjCHYledwuXdd1l9vDf/rrPW6btiuzVYBsrLZRjh57v2ZbNWAw9khV1B4qjxeWFqw1YypmamlAsCvJpmBX5qhQZqsAmZzqMmU7igVlsccq76MNu1ZIsAuMYFcSgh2QwthNp0pTHiZuBcWCQhh65BSZ4hSHowQ2VQfBTkawK4kU7OZ+LR5cGPut+rFm6s7ipEIhzsu9T6+PL0rHLeTKtSxqbpXQhGnttlYBhVGTnPzY1ckHa9VhqFlpimE2gh3KYuyR4pNbh5I9o5irLcpc84qUYDdNNM9sKdhmEexKYgp2qsvtYSkJfYapB+vnna5349ko5QiG1yv2BjBgoVDW4NX3/eqVRpurw6d4gPxWj9htHUo2j2LOtpgWJ7dLHIjmOe0F2yyCXUmMwU4dFZThRTlYZ3z4UD5dic8pH7iEMjUuytYq5TFjFcqmjUVCl3V28j3VIR7smGdVHwO52a9wUzut31CydRQztEV4Sb5USVqP8GCYb2mfq2CbRbAriflU7PSq8oT80HQxuKnX95al6a0wL0qfTX6GYId6WI4ZuDr5tuoQn1grHSA3c7AzhDi/oWTrKOaedSEVoTzEXdW5nAXbLIJdSSIHO+n486FgJy7JcnlS+K0DRDKPZ2LpmTs5wQ7tWu2R24aS0MHOWHzy4fCJdvWdcbZmEexKcjzYma8d2BXszIviiB3qt3KWx+OInWd1GIKd6V6LKH8ksFXZwU4+j2s8Favfk7T1y2HbQLAryYFgZ7ktSCoBqW+rd7ELw45rUfOxDeVIIdfYoR7a5W1yN3Z2co/qMNcs19ihdLuCnWMoCRrs5Moz39RhOEPlKthmEexKsifYLc+t3hWrfGhxfEfrvrtija0CCmQ+QWPt/tbKM1aHrWa5KxZl2xXsXEPJ5lHM1RZz/ZguIzdeH2Qo2GYR7EqyKdiJ/dUUu0zfY6f1ZssM9lfmRiy1rCxWbxVQInl/v7OTu8+tumv2HAcPUI+dwU6YVSmJPaOYoy1yQQpTaDd1OALkGcYlgh02YSxC8+jkACpGsMMmjHloHp0cQMUIdtiEMQ/No5MDqBjBDgAAoBEEOwAAgEYQ7AAAABpBsAMAAGgEwQ4AAKARBDsAAIBGEOwAAAAaQbADAABoBMEOAACgEQQ7AACARhDsAAAAGkGwAwAAaATBDgAAoBEEOwAAgEYQ7AAAABpBsGvJ43bpuuvdPdGf//4N/3au5H7tZpfbY8ur4ou2BQDFWym0ucT2V5ngfjWt63G7CHW0VvRALl6j0kwsnEN1tKlApInVaWsctgh2Ldkc7LZWjlwsahm4XzW+Xn6FAJp0wW6sGWVd+lBDKaFQOYLdpgJZGZeqHLYIdi2JHezGclnWMDwxPna/urXAgWIlCnbL8CStaxhphMFFqzygGOmD3bYCMQ1Twsx1DlsEu+Qet8vQb6bd9tiH5r24+nFA/PShdTBp36/3QcO8+ytnbrm6gut99dXxQfkfddCKtIX2Ksz7HiDYzWORVEPS36VOXd3og2LkHpVmAYLdxgLRxiV52jqHLYJdckO3u16F47uX210+3Cv1KoXyWUJ85Xq9uOd9uT1iBLvL7bH26vzJZ2lVfdWCmmQttLfvw8FOPBC+ltnqHH1QjGKKJdg1djJHgQwNth2xq3TYItglN3V8+QSmFH/kB2qfkz9MyT1wmXqZd6yN78+nruuu98OnYvVD3MMz7leNlz1wkAHRpC20179THT1uT13X/foKND75BDtyHY5JPirNM9uKpQ8Z7FYKRB2bhHKrdNgi2CWnfECQjmr10oEvdUp5av0ImTi98OpcGx8fL+J5ouGl91+GjitVtdZ2aZLbVa1486vadQ/GQ3xAKGkLTRyBPj5euu73uzw+bSo0tSGukaTOi4BQlOSj0mxrsTx9/Oy55M5VIHp4s15xV8uwRbBLTu0Y6o77Lmch073Xl9vDuMMXZjB+0Oi6rrsYj3tvKBX5oLTxr7G+6t4SQECZC208wdTHPWL3uF0oIRxWXrH0WuEY60g5mqDVwmqBaNFNz3LuTVUkgl1ytQc7vTn2YWd9UCq+QlCr8saqwMGOQ3UIpbxi6Y8HO58CMQ1CDRyPINglt7GENhz0Fhcl9GlbefTbT8Ua/xbrdOKr1vphYEIcaQttdUAKfCqWy+oQUPJRaba1WHxPxXoWiDvYVTtsEeyS8y8hr8tU51eNl6lehNtgf96eO/0au803T2j3P01PuF/V9wd8RQNiSltowoHwx9tzp182dKjoDF9QTKpDOMlHpXnFtmLpjxTOhgLR7vmT/r5ahy2CXXIbSmjlxnL11fUby9WrvLeeJNJujzBVhOXVlb8FCCtroQ3HFeIEO9sJreJHGxSrvGLp9xfOWoHIf5xh2BL/njqHLYJdcptKqFf6nbbvXl41fRWk3Gdf/wa4n9y3OaZX114GwklbaE9yoYU58GBqtv1CJQoKe2UdlRxXCu0pnNUC0T8rSU0y3oBRWZUR7NrnvsZua7AD0G8pq2PBDmiKZ3WsTgkHgl37CHZAcAQ7YAeCXQIEu/ZROUBwBDtgB4JdAgS7BlE5QGwEO2AHhqcECHYNonKA2Ah2wA4MTwkQ7BpE5QCxEeyAHRieEiDYNYjKAWILFeyoOJwKw1MCBLsGUTlAbAQ7YIdQwxOF40CwaxDBDoiNYAfsQLBLgGDXiH31QLAD9iHYATsQ7BIg2DWCYAekRLADdiDYJUCwawTBDkiJYAfsQLBLgGDXCIIdkBLBDtiBYJcAwa4RBDsgJYIdsAPBLgGCXSMIdkBKBDtgB4JdAgS7RhDsgJQIdsAOBLsECHaNINgBKRHsgB0IdgkQ7BpBsANSItgBOxDsEiDYNYJgB6REsAN2INglQLBrBMEOSIlgB/jYXQ4Uzm4Eu0ZErRyKB1AQ7AAfBLv0CHa1ClIPBDvAU6SyosTQtgTBjqFKQbCrFcEOSIlgB+xAsEuPYFcrgh2QEsEO2IFglx7BrlYEOyAlgh2wA8EuPYJdrQh2QEoEO2AHgl16BLtaEeyAlAh2wA4Eu/QIdrUi2AEpEeyAHQh26RHsakWwA1JKEOyoOLSHYJcewa5WBDsgJYIdsAPBLj2CXa0IdkBKBDtgB4JdegS7WhHsgJQIdsAOBLv0CHa1Sh/sqBycGcEO2IFglx7BrlYEOyAlgh2wA8EuPYJdrQh2QEoEO8BTjHKgcPwR7GpFsANSItgBngh2eRHsakWwA1Ii2AGeCHZ5EexqRbADUkof7BiuUCmCXV4Eu1oR7ICUCHaAJ4JdXgS7WhHsgJQIdoCnvMGOwiHY1YpgB6REsAM8EezyItjVimAHxJa4rBif0AaCXV4Eu2qkrweCHU6OYAfsQLDLi2BXDYIdkBjBDtiBYJcXwa4aBDsgMYIdsAPBLi+CXTUIdkBiBDtgB4JdXgS7ahDsgMQIdsAOBLu8CHbVINgBiRHsgB0IdnkR7KpBsAMSI9gBOxDs8iLYVYNgByRGsAN2INjlRbCrBsEOSIxgB+xAsMuLYFcNgh2QGMEO2IFglxfBrhoEOyAxgh2wA8EuL4JdNQh2QGIEO2AHgl1eBLtqFBXsqBycAcEO8JG+HPZNeRIEu2oQ7IDECHaAD4JdUQh21SDYAbEVVVZUGWpBsCsKwa4aBDsgtqLKyv0QKEfGctg05UkQ7KqRsR5WZwTaUFRZuR8C5chYDpumPAmCXTUy1sPqjEAbiior90OgHBnLYdOUJ0Gwq0bGelidEWhDUWXlfgiUI2M5bJryJAh25SqnHlZnBNpQVFm5HwLlyFgOm6Y8CYJducqph9UZgTYUVVbuh0A5MpbDpilPgmBXrnLqYXVGoA1FlZX7IVCOjOWwacqTINiVq5x6WJ0RaENRZeV+CJQjYzlsmvIkCHblKqceVmcE2lBUWbkfAuXIWA6bpjwJgl25yqmH1RmBNhRVVu6HQDkylsOmKU+CYFeucuphdUagDUWVlfshkFE55bBpypMg2JWrnHpYnRFoQ1Fl5X4IZFROORyZslUEu3KVUw9bZwQqVXJZUWUoRznlcGTKVhHsylVOPWydEahUyWVFlaEc5ZTDkSlbRbArVzn1sHVGoFIll5VjRiCxcsrhyJStItiVq5x62DojUKmSy8oxI5BYOeVwZMpWEezKVU49bJ0RqEVFZeWYEUisnHI4MmWrCHZlKbMets4I1KKisnLMCCRWTjkcmbJVBLuylFkPW2cEalFRWTlmBBIrpxyOTNkqgl1ZyqyHrTMCtaiorBwzAomVUw5HpmwVwa4k92s3u9w+5C74/qtTPH38TK8+3p7HJ693qe8+bpfu+fPD2cs/Pl667vfdUA9fr9Na+r7/+HhRW9B12ow/t0s3rzHr1gTWiBX3/PmQS0B80VFxr3+lsnJU3PxwqLh3Qz16VZw848/t0nWX2yPjZsR5yCXzIXdyfZC63B5zRxVLRi2H588PZwh73C5Dz9de/boKhekoGWGun7fnsd5bRbArxeN2sXVHZSDRh5n7daqxv7+77vL2/U/s9HMVxQl2UqES7FALU8Vd3/+Te7JXxb3cHr1PxQUKdl0nfX4j2CER2yDVi4FJNge791/TZ6e/v7vuZRik+qnPDyXTxwp24yDVE+yQ1vg5aPkoMzzx62t8+Lg9aR90TD1b6+Xakb9jwU4djf78/d113VyoBDtUQq24sSdPFffn+/PiV3GvwujlrrhdwU452DBX3EWqOIIdohtL5nrv+6E3TiUzPvz+fNIOYE//Hzt2b+zk2pG/Y8FO/Hj2rx+jZCd8ACPYIY3HuHMWu+P7r27Zrd+vQnhyDDM/t8s8OH29amdm++DBbsqgUgMIdijVnNue5Ir7M557nQYP+eB3b6+4t+euu959Ki5MsJsGVKniCHaI7SF1tD/zcbh57JgOxfUrwe7n7XkcpPr+/mr5+NQHDHZTyUz1TrBDGtZgN4a5sVuPnzw6+Xzrz+1iOGL3/kv4LBU12A2niccwR7BD0daC3Vhx4wjhUXHzEbvVigsW7KQqI9ghCWuwG4+Ejd11KRkx5P28PRuO2N2vwjHyqMFuCHPj0XSCHRIZjnILRwjGw97jMKNfxy2eRTJcY/f9+aR9eHIGOyuPYCceXCTYoWjzePCqxDX5ogL9MnB7xb3cHr1PxcnBbkPF6csRDi4S7JDGOEgNXU0smSHYGUtmPoBtuMbu+/MifV5aDXYrJeMMdkPzhpeWYKfUZjMIdqXQOu7l7WO+yHS4jls8FSseJDPcFTscPJDqYTz0TbDDqc39Vuv2L7fbfMHDcBn4xafihouKfCqOYIeqaYPUyzBI3R7TCVbtKJ1wyal6V+z7r/EahqUcxksaCHZHEexKshyXe7k9ht26ePeQ1AVtr/b9cvBAuE5IvMTBGOyOnIol2KEaUtcVKu7texg8zFey9vaK+zPe27RecXKwO3IqVgt28hesADH8WY7SjSUzFIVwV7g+sliOyX1/Pg0zLjfJSidVe0OwO3IqlmCHfCxpaVOwewxXM4xdefrAtBwJDx/suMYO1TCOQFpa2hTslkvu3BUXLthp19gR7BCfJS0Zxw53sFsuuRvK5GMuQNMdsiGCnfkaO4IdorANM8P3m4y9/PvzIvRd906/v1/n5xMFO+6KRT2swe778yJU3JN2M5PzHtWrT8UFC3b6XbEEO8SnFs735zBIPcbTREvXXYlZf3/P/T9RsLPcFUuwQxRz73ztpMsRhisKpH231K3vr53wLXfLcsT79ZIEu/HI/Nw2gh2KZqu44QIg6UtMhBHCUXFvz/K3sMYOdmPFXY0VR9Ehmvvr9Mmn73uxZISHYq/+epUvm/sjzDh38hTBbhqkpjGOYIfI5l5luJ5a3Fkbbos1XQwkHDyYHoa9xs6MX55ALVwVJ3638HIt0WrFCbUT/ho7M/WXJwh2iM82SPXq5XdLyRguv5sO1wkPw15jZx2keoId0tD67rLnVuthOO49vWy6jMDw/ajed8XyW7E4BWvF6V9BN1wOsVZxht+KtVRcoGBn+q1Ygh3iU5OTfgBbGqSMF2QbvpHY+67Y8L8VS7BDFPpO39it/R/GmDLgKoC8GigrqgxpVFEOAadsBsEusyrqIeAqgPQaKyuqDGlUUQ7xGlMvgl1mVdRDwFUA6TVWVlQZ0qiiHOI1pl4Eu8yqqIeAqwDSa6ysts4I7FNFOcRrTL0IdnuJt6m6fs/nfu267nqfH//575/44yqddiVpgfWwY0bLD6KLW03YKF6b1DkvGqfW0aaJ//ynXNY9Fl1fW1lZZzT9dfr4dL86t6H1Z5v4wbLa+A5PKxP/+U++19XyXSSOh/5TJqg415TD11gKFbRhQ8s1pX+DRfr6IdjtYdoBmveX43ssjjHyLtg2zPQl1cPWGcdbk5RgJ1eOsl1WNqnhZbLdieh1tGli0xcxWL+LpNiyss5o+C4k6WsjBmMNEexa5z88rU5susnU/KuSGcshwJSm/YNPpzftlwxblGBXBe3dHJ7Qamd5h4WX3n91nZTkvuRvfSyxHjbNuIwycrAb/vC5iw+TTQ/vV7n7K1tUnlh9iLYZ62jLxA/158nVLyytoKzsD4dvLxe+Y2+swOUrHnrxKMKmD0RqXaJ8K/vSLRMv30g8dKR58OpLKofDUw4n0ISv3Ju+A9xdKub90uN2KeCoA8Fuu8ftou3sDAdkxzdYeUUtlSrqwXvGr9duqPyvVzXYjV/0/1A2kG3QkF7UhxfnvGiJuY62TPy46T925P6Ny8LKyvnQ9HuD8s9Mzx9EN53O7ol1jdi0s5QnVrvc9LthhpdqGMjML5n+qPdfa6Vi2y/dryXUDMFuO2uwE5+b32/5jb9fO8tFdSXXg/eMX69jmNOCnfbjm3/G75x0nMKetqdpg5fxuQixWepobWLLz3/Jwc7ye1/lldXmGcWfaer7+3XcbtuCHbGuDZvex/dfy9Hf3tivxCPfZZRDjCmH3wb0uHRerSkx72U6Ddv3BLtdtDOB47u4PCOkdj3Yvbx9S1+Q3eLNE+ZgZ/pVDFNgUzavKcRxyO4MrHW0NvFasHP93lfBZeU54/3VdI3dtmA3FCIfnSpm2peuT2z63bzlbj9r5stYDhGm/Hp1bjXHfsl0vWuGMiLY7aFdHnm53a7ymUP5erDpjX3cLl338vSsvvHKeaLS6mH7jFqw+/u78wp2phu0TOMRwa599jpandgZ7H7enl2/91VwWfnM+HO7LH+7Mlz5BzsO19Vv952xenobr7GZX+1LKofQU477B2VGeVPZ9ktDNBC3dp6RimC3l3ysdXn31PdRD3ad/NN1Q80YfziylHrYPuPuI3by9h1e5IjdCTnraHVie7B7vD1bf2jc8TDGlBFm/Hl7Vi94EF71D3Zbr8ZDHsohBssOcUNKt93MJ+3Y41zDkKDieu0+X+Xi1KGCjMdZxC3uu18yzpEEwS6M+e21fmPA8EFnOhUr9rmPjxf9Y1BfUj1sn/HoNXZ978jK86uMPM2y15FhF7lSdEInvGjf8VZVWTln/P58Wrlj0TvYkesq4RfsNmQL23ikjFbLDaQFlMOmKV3B7nF7sl8ZZdzgHtueYFcv4fjryhhj6jqnCHbTXbEfy5Q/t8t0iff3pyW6Dc9xV+zphA92f39vGpCKLCv7jH9/e1z/5BvsyHUVs34MNu0utRek8Ui4nVwZreoNdtYp//7uust8zEWfUdhg9v1SMccgCHY7KDljfK9dh8GXs0LDV2ppp2Ib+h675Y8yfY/dvJuQH/5ol2pLG5nvsTu9lbAh90DtVOyUe7Qp6yory4ym8wCmKT2DHZ+aqqbf9uI4Fesej6yv1nsq1vzS2qc++9Y2XmOnbfzkn5IIdnts+W5pJdid4pcnjMHO+MsT79LgZN+m/PLE2Rm/L2p5Qu6BSrC7v2p9q+u6Rr7Hbohrpr/u3TilUje2n0SivOrl3pcqhWP40QXh0K/l1b6kcjg85Zdl/zANT9Odf9Ztre2XrBs/FYLdTlLS0N5zofd8qcFO7Unmj9p9SfWwfUZTsFOHWGlYXd+m/FbsqR0IdsZvINB6YA1lZXpo+jWk4a/bGew4YNcAx75ULxzpWMPvd6WPmV7tSyqHo1NaK0gKdpbfX7Z+ZUPeoYpgF0UVw0NRq9BfAvzV2OeztA0QNVAOWRpTOIJdFEV1wSpWob8EuNXe5wtsG84me5dLuYqAjSkcwS6KorpgFavQXwLcau/zBbYNZ1BUl0u5iniNKQ3BLoqSu2CZq1idElDU3ucLbBvOoKgul3IVaRpTAoJdGP7vuuU6zVNwbxmGHLjt6zy5e31m/hsq51uLhHz6Q+5umxPB7lyCxPncnTanINuwnOJBbP5lRbCz8d9Q7i2MihwvnNzdNqeDQ3wJQxXBLrXcnTan3NseZ5S712eWe/OjSrm7bU65t30ALfwNdcndaXPKve1xRrl7fWa5Nz+qlLvb5pR72wfQwt9Ql9ydNqfc2x5nlLvXZ5Z786NKubttTrm3fQAt/A0AAADoCXYAAADNOFewy32IF3XL3X9TyL2N0ZTc3TmP3Fsd1TvaA4P041oc314BFdWYnvasKa09kRT4Z9IkHzSpHKX94bTHrb32lPX3xFbU+1dUY3ras6a09kRS4J9Jk3zQpHKU9ofTHrf22lPW3xNbUe9fUY3pac+a0toTSYF/Jk3yQZPKUdofTnvc2mtPWX9PbEW9f0U1pqc9a0prTyQF/pk0yQdNKkdpfzjtcWuvPWX9PbEV9f4V1Zie9qwprT2RFPhn0iQfNKkcpf3htMetvfaU9ffEVtT7V1RjetqzprT2RFLgn0mTfNCkcpT2h9Met/baU9bfE1tR719Rjelpz5rS2hNJgX8mTfJBk8pR2h9Oe9zaa09Zf09sRb1/RTWmpz1rSmtPJAX+mTTJB00qR2l/OO1xa689Zf09AAAA2I1gBwAA0AiCHQAAQCMIdgAAAI0g2AEAADSCYAcAANCI5oLd43bpuq7rusvtYXj2el95Mn5j5iel59M0ZvXPT7ihfNqTclup7fFuZIL3Lg2vPzZ1gwrasEVtH6E0TAWTo0n3q16rxb2J8Rh3IMJr6TaC+Y0wNcc+SdL2JN4+5toxT5Bi+6y2Z/v2aSzY3a/TGyFW2f06bpD5P7YnEzSmv1/11aVpjLho/20Sp21e7Um4rdT2+Dcy/nuXht87krhBBW3YoraPsM5lN5O1SY/bZR4DldYV9CZGdL9ax92kG8HyRqgNSvZeeLQnbScx1446Sbq+ut6eHdunrWCnBqhhG83/EScwPpmgMeKKZ3EbM6T9y/VqbM/aNgnfNu/2JNpWpvYc3GiV8X9HEsregEVx20cqjHn15Wwx9763Rfdrd7lY/sCMG8G0B7U9m609abePuXbWW5mvPXu2zwmCnf+TCRrTP24XbV2RG/O439XDhlk3lG97Um0rQ3sObrTaeL8jKduUuwFiWwrcPovpY3wxTVqO0RTTpMget8v1bht0820E8WCZ/HSeN8LcnpydxHgILNv2Mbdn1/ZpK9gJ22X4kG3d36XoTKbGTJ/spLP3aXp2McHOtz2Jt9XqNXbtBrtRUcElewN0RW0fsVHlpKjxWqC0e7bcxiBQVLBT3giJab+asT3FBd8c28fRHoJd3y9XGV5u93F7ZByPzY2Z1zWdySfYWVedcFsR7IoKLtkboCtq+8ytKKxJ48pT7tkyu18NwVqQcSPYDkjp+9XM7cmwfaTaUV/IsH2c7SHYLXKfijU2xr+FwVUX7CRp2+PfyDqHrukTh353bzEpIXsDdEVtn/FNFNZbQJMmdVfHKrF8lpE/X7Azl7PX6gpoT47to9aOVQHtIdhJls1RwDXv5lWkvcTY1j9ybaj19ojib6vVYFdCR4pq2zsSW/YGaEraPqbP9+VssfPcPHGfv8ljpv2Rxd08sWWCBO1JvX2sx8YMUmyftfZw84SwCcSDvnfT3cLGJ6M3RlzZY7lGM35j1D7hv01itW21PYm31Xqwy9WRUvH4Y1PK3gBVOdvHfnFQtiZl3LMVwj7oJt0IljdiwwSJ25N4+1hqx9CaPsn2WW3Pru3TWLATPkFJ22p+VtwqxiejN2Y5DJu4Mep+x3+bxGmbR3uSbiuPYJetI6Xh9cemlL0BsmK2j36YSPyUkWuLZduzFULL/cJeP+lGML4RYnss71S29qTcPrbaybV9fNqzY/s0F+wAAADOimAHAADQCIIdAABAIwh2AAAAjSDYAQAANIJgBwAA0AiCHQAAQCMIdm2TvyWnkV/0ARKjjoDjqKNECHbtGopoKR71MYB11BFwHHWUEMGuVcO3Z8tfVD2UUuvfAA+EQx0Bx1FHSRHsWnW/Gj4OGZ8EYEMdAcdRR0kR7Bpl+oBkfhKADXUEHEcdpUWwa5SpZh63C2UEbEAdAcdRR2kR7BolFNLw347LGYCtqCPgOOooLYJdo4xHubmkAdiEOgKOo47SItg1ynz5wvAspQT4oY6A46ijtAh2jbJdl6rfYS5+ZySHxgGRfx0t01NFgMy7joQTtYxH+xHsGuUupPkz0v0qTCU9AOBbR9KzlBAg86sj6QgexXQAwa5RtkKSSke7Lel+5cA4sPCqo77vxSPfjEWAzKuO1CML3Da7G8GuUT6Hvh+3i5Lj9GeAM/M8hTQVDkMRYOA9HomTMBztRrBrlKWQpPHIMApxyA4Q+NSROjnBDpBtrCPhVappD4Jdo5ZrUJecNj43V4qhbiglQOBTR+rk1A8g21hHveliB/gj2DVq+YQk3mYklwnBDnTGGz4AACAASURBVHDzqSN1cuoHkG2sI+6cOIhg1yifn+HjVCzgtvHnLAl2gMGWOuJY3XEEu0Z5Bzv95gnGJWBEsAOO864j8/cIYSOCXaO8ComvOwGcCHbAcX51tLHaYEWwa5TnNQp8QTHgsPFaH4IdYOBTR3y7STgEu0b5/74yPykG2Gz8nXKCHWDgUUfSj4ktSHp7EOwAAAAaQbADAABoBMEOAACgEQQ7AACARhDsAAAAGkGwAwAAaATBDgAAoBEEOwAAgEYQ7AAAABpBsAMAAGgEwQ4AAKARBDsAAIBGEOwAAAAaQbADAABoBMEOAACgEQQ7AACARhDsAAAAGkGwAwAAaATBDgAAoBEEOwAAgEYQ7AAAABpBsAMAAGgEwQ4AAKARBDsAAIBGEOwAAAAaQbADAABoBMGubo/bpVtc7wcXd7+GWQ6QwVwNUgcuolOLjdAb9LhdKDrsNPYnVaIONRXd5fYwNkt7/viq0hWzuzCL2LOYEewqJqe6ED2+4J4KrFkKQhxNiujUjmBXRPtQL0uwCxuq7KaqUzpwzFyXplpWC7PgyiXY1UupnBCFVHBPBdaY9/xFdGpHI4poH+pl6EDWw2gxGJNdtFx3uVwSlQvBDjmopaOXknRITy4x6VOeuoilp4pLSLKTAHYzX5jgGPe0fj0NHrf7NMX838vtMRfN5fawnhharyzh/+rBlutVb2zBwwdKYOwg4pOmXj30S2shzLMLPdTaAw3JTh+MzHXhaJt5Jdf7PMfD9LppGZaX5r/R8IFQK0x5Z9KJ+4MSK5NgV6+VI9OGE7XKwT39Jbmn6pOR7VCweZd/uwq9Vd39ms5cTS/qRSPv9m02VZYr2N21waLg0QNFWD1iZ+rVvbsQtl23p2WtOYX11sWNE1vaZlvF9W5Ym3Ows79k+Rsvt4cl2OW9mHETgl3NtD6rH2pTykevC+nTluFjnjLkldiJgb7vxZ4t9lbzpxV1BBsePwyHDeRbMpTjdFL5eFaW+xo7pdDuV6U5gMx+jd3YiUy9eqUQpkfqJxb34TQlTCld2lQX5rZZVyDNM8+xPti5alL7G92FqW/VEsdEgl31zAeNfT7nS7lQHXJsxylK7MVA38s7/GUkkErBcOmP4dOM8eiHnAyV0UAuC2dlrd08IQ1b5DqsMQc7w0cTwwl+WyHo3do9okid1noQQK8L260XrsVrD3dewLp2cHz1JqeCD6YT7JohXwDk6nPGzEawQ+3k/f2UiW6GYGcb4kzX72wYVLwqa3XAENZArsOq1Xxh6tUrhaB3vJW7IYR16LnOXhfGttkWrrFcbOG5bQh2KJD527C69WBnPvlkCXaMKaiGOYJJ44jfEbudwc6vsjy+7kTOpNQgXI4Eu1BH7PRD5PqsxmN0PuOM/SpX06Xh3ttmX7DT7wch2CEg7RS/dDmBer/70ivNlzjoFWI+A1ViJwb6vrdfw23ox+5r7PYFO9/KWg920tk1Sg5Ou4Kd5zV26gVljgA2JburmutcdeET7Fa/TsU+2Lle2hrsuMYOaZg/yGj9UH3NfEmGYcjhrlhUxTBILH3YVReWy7yNz3kcsVupLPP4IT1naDdgti/YuQthzw2gwgId362lLGs92FmuwpMOoFkHO9dLnsFueY67YpGK3NfU8hBfNX6/3XIdT3e927u6ceFAWRyxzHw/hNavD15j51dZjirTDzaWOGqgKDuDnfCKXgjzMuUe7Wa/w9VaF+vBzlYI2jFEx1hlfmkt2JkKU3qm4LNYBDsAKE3BgwbaR/erG8EOAIrBQXLkR7CrG8EOAIqhfAMykAHBrm4EOwAAgEYQ7AAAABpBsAMAAGgEwQ4AAKARBDsAAIBGEOwAAAAaQbADAABoBMEOAACgEQQ7AACARhDsAAAAGkGwAwAAaATBDgAAoBEEOwAAgEYQ7AAAABpBsAMAAGgEwS6Lx+3Sddf7vpn//Pdv/uea7n7tZpfbw/XyWlNWFgWc0+ZCFot3vYRl9+umdd2vYmU/bpfOiHrGJocGr4Ff/98yQknd2zqxuYJaHN0IdllED3am3biwPsPL1tasLAo4r3TBbixD73WNoxXBDoElCXZbRijTxIZebaygVkc3gl0W22pDKQOPUeF+lfv2sJefVzg8nF9XHhoWJTZWWRZwXiuFrGe4fcFuOabgV3fLcOWYXt1JAD7CBzu9HLaMUOqINHR+ZWJLBTU7uhHs7B63y9A/pk4x9pW5j6gdzXnoWOpZem245t0e7Ox/y7wuqfHSy/YZxbbW3/dxEvkK+T1AsJvHHu+yGxvlnp5Y17RiBq+BewjTymHLCDVMLb2k9Hx7BbU7uhHs7IZ3/XoVjtVebnf50K1yEEyydBj1eO+wUL95+xDBTioUU2+2fgqzdn3GBFQiZyG/vH0fDnbj8j1HnHky1/SUcOOKGbwG24LdphFqaqH9iJ29gtod3Qh2dlOPFnvF8p5L3Uc9fT8+lj4lzX1FWa526l/pp4eC3VyVQk81lYj9A5F2FHzcDtV3fZxFvkJ+6rru19dqsPOLel7BTjh6YZ9e28OgNcUMXoMdwc57hBL/vom5b+sV0ezoRrCzU/uoEuWFnmbqzcvUepcUp7d+Ovn9Lvf791/aR6PJ08ePZUgw3fBj2uE7yka7vPRyu6nHvoFy5Svkj4+XuZDnWl4tZMuf4RHspJNS1unbOCYBl8yDl9rtHH3+cnuowW7jCGU6bGisE8NyWx3dCHZ2aldSu8WyFzX25fllW3canjIcyh6MJ3H6zRcoGAmFvf3zkNDIy+2xNjFQlDIKOfoRu7U/03c5aEDmPq+ODjGP2GmnXo13T5g2gvh0Y6Mbwc6ujPGgDxPshL9m4xUMOkYG1KSMQo4d7KzfZ6LUOtV7BvUHuyNXgVvCmffFDPXXB8HObmNtbDiaLS7K0mH1fr/hVOxKX992z5Fxs9T/kQanka+QjRku0qlYz2DXxriFFVkHL922U7GHv7dhf7BrZXQj2Nn514bX9afzq8brT8WuNDyjXmO35YidXqxSqSiXjGpXkIpMF2e00PNxGvkK+e2506+xi3Uq1mv6Nk40YU3mwUvtdluHsM0jlH4q1nbOVb95osXRjWBnt6E21u76Vl5dv2NcOgi351Ssuz2GT/fKH6bVuWVJQOnKKOQowW71K01M59g4YNe+rH1eHx82D2EHRyjzGOVz80QjoxvBzm5TbfRKH9H2ncurV8N3PMr963o3nIrdFuxW2yNVpOFAvLV1jAqoS5JCHmvw+/NJKOTXvysZzhbsTJOFCHYcsDuJrIOXbtcQtmGEktvgvHnW8cc1NLoR7AoVINgBSOVAYgv2y2NAmTYNYXT74wh2haIMgIoQ7AAbgl1iBLtCUQZARQh2gA3BLjGCXaEoA6AiBDvAhmCXGMGuUJQBUBGCHTA7kuTo9scR7ApFGQBl8glnBDucGcEuL4JdoSgDoEwEO8CNYJcXwa5QlAFQJoId4Eawy4tgVyjKACgTwQ5wI9jlRbArFFUBlIlgB7gR7PIi2JXC1rkdL1EVQHoEO8AtYLCjCnYg2JWCYAdUgWAHuBHs8iLYlYJgB1SBYAe4EezyItiVgmAHVCFXsCPtoRYEu7wIdqUg2AFVINgBbgS7vAh2pSDYAQU6kroIdjgngl1eBLtSEOyAAhHsgFXxkhwdfgeCXSkIdkCBCHbAKoJdUQh2pSDYAQUi2AGrCHZFIdiVgmAHFIhgB6wi2BWFYFcKgh1QIIIdsIpgVxSCXSkIdkCBCHbAKoJdUQh2pSDYAQUi2AGrCHZFIdhl49mb/aekDIAjgqcugh1OgmBXFIJdNgQ7oCgEO2Afgl1RCHbZEOyAohDsgH0IdkUh2GVDsAOKQrAD9iHYFYVglw3BDigKwQ7Yh2BXFIJdNmmCHVUBeCLYAfskC3b0fx8Eu2wIdkBRCHbAPgS7ohDssiHYAUUh2AH7EOyKQrDLhmAHFKXGYEfUQwkIdkUh2GVDsAOKQrAD9iHYFYVglw3BDsgldjgj2OFUCHZFIdhlQ7ADciHYAQER7IpCsMuGYAfkQrADAiLYFYVglw3BDsiFYAcckSvJ0dt9EOyy8ey+/lNSBoCn2OHsyLwHVwEkkCvJ0dt9EOyy8ey+/lNSBoCn2OHsyLwHVwEkkCvJ0dt9EOyy8ey+/lNSBoCn2OHsyLwHVwEkkCvJMcD5INhl49lf/aek3wOeYoezI/MeXAWQQK7oxgDng2CXjWd/9Z+Sfg94ih3OjswbYxVAWMmyGgPcDgS7bDz7q/+UlAHgKXY4OzJvjFUAYcUblRjRjiPYZePZX/2npAwAT7HD2ZF5Y6wCCCveqMSIdhzBLhvP/uo/JWUAeIodzo7MG2MVQFjxRiVGtOMIdtl49lf/KSkDwFPscHZk3hirAMKKNyoxoh1HsMvGs7/6T0kZAJ5ih7Mj88ZYBRBWvFGJEe04gl06+/rr7hkpA8AmeHIKOG+MVQAHxRuGGNGCI9ils6+/7p6RMgBsgiengPPGWAVwULxhiBEtOIJdOvv66+4ZKQPAJnhyCjhvjFUAB8UbhhjRgiPYpbOvv+6ekTIAbIInp4DzxlgFcFC8YYgRLTiCXTr7+uvuGSkDwCZ4cgo4b4xVAAfFG4YY0YIj2KWzr7/unpEyAAYJklPAeWOsAjgo3jDEiBYcwS6dff1194yUATBIkJwCzhtjFcBB8YYhRrTgCHbp7Ouvu2ekDIBBguQUcN4YqwAOijcMMaIFR7BLZ19/3T0jZQAMEiSngPPGWAVwULxhiBEtOIJdOvv66+4ZKQNgkCA5BZw3xiqAg+INQ4xowRHs0tnXX3fPSFUAgwTJKeC8MVYBHJRsGGJEO45gl86+Drp7RsoAGCRITgHnjbEK4KBkwxAj2nEEu3T2ddDdM1IGwCBBcgo4b4xVAAclG4YY0Y4j2EUUpIOGWg5lgNNKkJwCzhtjFcBByYYhRrTjCHYRBemgoZZDGeC0EiSngPPGWAVwULJhiBHtOIJdREE6aKjlUAY4rQTJKeC8MVYBHJRsGGJEO45gF1GQDhpqOZQBTitBcgo4b4xVAAclG4YY0Y4j2EUUpIOGWg5lgNNKkJwCzhtjFcBByYYhRrTjCHYRBemgoZZDGeAMjqSfQuZNswpgk2TDUMCHp0WwiyhIBw21HMoAZxAw/eSaN80qALdk4068h6dFsIsoSAcNtRyqAmcQMP3kmjfNKgC3XONOwIenRbCLKEgHDbUcqgJnEDD95Jo3zSoAt1zjTsCHp0WwiyhIBw21HKoCZxAw/eSaN80qALdc407Ah6dFsIsoSAcNtRyqAmcQMP3kmjfNKgC3XONOwIenRbCLKEgHDbUcqgJnEDD95Jo3zSoAt1zjTsCHp0WwiyhIBw21HKoCZxAw/eSaN80qALdc407Ah6dFsIsoSAcNtRyqAmcQMP3kmjfNKgC3XONOwIenRbCLKEgHDbUcqgJnEDD95Jo3zSoAt1zjTsCHp0WwiyhIBw21HKoCZxAw/eSaN80qALdc407Ah6dFsIsoSAcNtRyqAmcQMP3kmjfNKgC3XONOvIfnQbCL437tZs+fH3Ine//VKS63x/Tqz9vz+OTrX6l3fny8dM+fD2fffdwuXff73fDq16uwlo+PF7UFXdd113dpyY/bpesut0feLQls4Zt+/v7eVKFPHz9zXSgVOq9iqNAP53o/Pl7EChUm+3qd1jJNppNm/PPfz+3SiWs0rtTWEsDtz3//9DIRe5StTPq+Fwey610bpOxl0s+l1P1+N7x6fxWK8XG7OMpkmvfn7bkbhs7zINiFZ+ptv+9LB116/GyOXO+/xi745+/vrnt5+5b2+69/Vz6UHAt2ylhFsEN9fGLNjgqdx5L7dRrhpgqdVjFWqDtOHQt2UoUeDHbGyYCZqROKmclaJn3fDwPZVCaXoUyGfn51lkkfINiNZdIT7BDOeLBu+ZgyfOi53seH359P2oeYXt6599M4Mf1/ORjQBwt24kj2789//6aGz2PVY6iHD3b6KJVnXpGf+boqCclUocrxcnHeaVzRothaxuq3BTttsvHwycs0TBLsEJIeoTqxEIbu9+trfGgvEzl+iSXz7+PjpbvcHD223xzstPFuKpPbg2CHUITjXHN3fP/VLSc6//4W+pwj2P28PXfdr69ePhjQxwt2U+lOowvBDqXbE+y+Py/C2GCrUOFQnC3YjRU6jzficYgowW4atOYGEOwQkNQZHrcn6QqE+dzrlLfsZSLHr5/bZSyToZNbDmqEC3ZTmUzjHcEOxwzDxpNwpE0YNsbD0WOoWi7CE8+3/rw9G47Yvf+aCyNysJtq4INghxoEDXZShb4vVxe9CCeSfm4XwxG7+3Wp0LjBTqhQgh3C8gt244DlKJNh7NCP2A0DmbqiGMFuKJPx0CDBDsfMO2jpmNzY9cdhQ7/gVDxva7jG7nF7kmpmNdhZeQQ78TMZwQ6l2xPsxlOxQlzzqND5aJzhGrvvz4t0FdFqsLPyCHZLhRLsEJYeoaRjcsKVAP1wuMFSJr3xGrvvzyfpKMZqsFspE2ewE4/BE+xwjBCbRC9vH/MFpMMFp+KpWPGzxT/9rtj5YMCyWOHYQE+ww4ntCnZ6mRgqVDjH9CMXgnpX7Puv8fo8pUKNzSPYoVhrfXUokyGcPRxlsjzsum46bDGfd1qqz1QmPcEuBIJdSEsHlQ9Qf3y8yHcGSV3Q9ep8MGC5SVY6Ft0bgt2RU7EEO9RkX7Dr+361QrVhxnLV3ffn0/A5bbnkSLrtyTRYHjkVS7BDLIYespTJZSoT61E316vT4TrhJllzmQiLOnIqlmCHcGy7TunS7A3B7ufteUxj85HtXjzKzTV2OLf9wU5+Rq9Q44ilzbtcFCtW5XJBRfhgxzV2iMXRQ3r95okNwW655O7913J6ylgmfZhgxzV2CMdcGN+fT93UyUy3Vnx8vJhj39/fc/xKFOy4KxalOpJXVp4xVahxmFHn/ft7Hn4SBTvuikU0jh7Sj0fdpmHIdGuFNY3dr3NnThTsLHfFnqe3E+xCmnfQQjcaLzWQOpkUqr5eO+E7tJYOKt4hmyTYjTfqzgMYwQ4FOZJX5Ge+rp2Yh8wVKsS4r9fOeFv6z3wcok8T7MbzYvPzBDuEpEco/dJS6WO/pUzkRYk3kicJdlOZ3MUhmGCH3YTYJBN7lfgjLdNVPoavtRMOBkxzhb3Gzkz85QmCHcoRLtiZ7jFaq1DDbelThQoFG/Yau5UKJdghLFNfVcvkQ+rw5jKRFvX3t3Q+Kvw1dtYy6Ql2CGLuNFJJ6B9ivj+F/mjtvobfihUWGCHYqb8VS7BDOQIGu14ZtEwV+iRXqDGHGX4r1rbAMMEu5G/FEuygsAQsqUyk/mMpk14edwy/FWsvkz5AsLP+Vux5ejvB7hDHvjLUwzRr2fEQSOlIXtk3Wa55M64CpxJvdKjiYcMIdock6IJp1kJVoHBi9zP2SeMzRybLNW/GVeBU4o0OVTxsGMHukARdMM1aqAoUTux+xj5pfObIZLnmzbgKnEq80aGKhw0j2B2SoAumWQtVgcKJ3c/YJ43PHJks17wZV4FTiTc6VPGwYQQ7jXxB5vVum+5+dd7HEOphpMVueCje/aTfSzF+SUrXdV13ua18CeSmiVEl8T22V89IKjZtaulV4W5QV7eUpjmWfnLNm2IVj5tw2bt6n5Yy78p76n4PcZRfQVlvEDXsZ+/XcUmWTmL7Bq6qH34Jm9H8HcvmTe4ap+5X/T0pY4wj2MnEd8X57owTth7sTHfnjVXRm3cm1l3PpolRJcN7vG0oWqY2vTpnO0e3FPuwsWMbn/F/MsG80Vdh+KIK8/ez/Pnv38p7aniVT2zh+BeUd7AbJ7QHu/df43pSDjRxH0p37467C/37+e0b0rzNxwzgroZMYxzBTjS8L0Il6G+dMF3XNR/svubvqBxeHWp++Iogw9YZnnBUgbZpSXYtGd7T+U1WHhonFjrA0kPEnjb1RvHrUr9ep473R0wqv776QpLTsXkjr2L4BmbhNwzVL3SV/rnfU+VN1HagOGRTQRlnViZeDlxYgt3ykamhYDfsTOZv+1IezhOLW2h1UFsygDJpGWMcwU7wuF3MlSC/NY/bpeuu99OcihUfij+JJvz4kmtzbdnaqJk+kjjeY9OwMxE/URh6o+l3+YZfsTRGk157Un/G/8kE88Zdxffnkzaq2bfel/s9vV+VVzfsALBmU0GtzjtHluU9UjvJ9+dT173+bexUrPJbUP/kX5JdJrZt4JUM4Ojt+cY4gt0K+57rjMHu/Ve3/MyzNdjt3PGgbqa92LgD1Cd27hKHY0jiRXW9V7BTf8HF2JONz/g/mWDe9Kuw/Qba8FXqjvd0OG7BEbsoNhWUzHy0bpzPFuzu02+CtRXsvH+f3bbNtY0pZ4CVy/DzVAPBzk3NdcLjUwW7n+nAs/ijMcMZsTHnzZtkrSvPx7DZ/7fFNOZYP7JOhSRelTLPO/88l3gtnZDzxo63xLjxZCLBbt9k1t9AG4Kd+z1VL0vmcF0omwpKm8j4Rizlo123cL/Klzq0FeyEXxj792cMdsLVCMuxCe10tzaoqRnAsJ3zj3EEOwetsKSYd6pgJ95SNB2xM1zDfrnd1COcmjLuGkJwpr2cbRwanr9YLrz/+HjpupenZ/XV+aeEtI738vbxm2C3a7Kft2frb6ANidn5nuq3m5HsAtlSUNp85mmGnbkh2C0/4fqvtWA3Hf4XX7UHO/0GCHlQ0zOAob/nH+MIdjaP28V9ecOpgt38ULteYbnDbj4A49mVORnbli0HGKa9p/qxaXhmym3ir0/eX+VnhI43H9sj2G2d7OfteYnLhnlXjthpp145FxvQziN2rvODyz5cDHbqZZdtBbttR+yWTWgY1IwZwPVBJtsYR7AzsVaUkXh2MnwHjbTY3Q+NI+j86nwVtvdmZhBoxeZr7MZp5b3ty+2xnIoVO9jwqvLk/H9untg82TTgueZ1X2Nnfccp6hD2XWPnvnr1v396sDN9eZA0tPmPDiU+9L7GTplX36L2DGDt8LnKgWCn0W6XGJwx2Mm3RwglYQl2w9cFPX8aujFjQPu23MQnDFDy3vbl9jBnDlewEzqe2C2NHdv4jP+TCeZNsYr7VdmYlnm/rp10e/J4ua1wwI6ijmbPXbHug0jDe3quYOd9V6wyr7jJhQN29mBXUjkQ7GSWVDdTd3mNn4odvvJK/OpXsUiUb7kbJh4/G2n0C3o5FduaLV+7tfSHuaCErrV0vKk3at9jt+yppY4n9mFjxzY+4/9kgnmjr0K6oGpl3uFNnE/Xzt8BZn6HxW+/xWHbv8dODRKm99d880SvVGIzp2K9v8du+duXYxNjh7ZscyVFFzTGEexEpt+d6JYBpj9dsDN+Sb3j5gnp24zVK68NW5dY1xb3F+Ubd4Oy5VIYw5fFO395Qux4a3nFc7Jc80Zexfj1zvpe7i6U7fLV0PKPjynvqfEYBlUdzKaCMj1l6ganC3YrvzwxjVN/bPsWa4fW3oBixjiCncCW684c7NSqUDeFVAbT1uiNwa7nZyXPwP7TlrZxaKRf0iBGEPUIk6vjrWUdz8lyzRt3FabfExtK2xzs5mPzwmQSqagJdcFtKSjtzJ+pG5wv2P33z/VbsXKw+2Pat7jeGtfF+NnGOILdiuw9soQ2HGw/4CD2HL0vbXoy7GS55i1/FaiF463k4aaHdSHYrcjeq0pow8H2AwOxkxh7jvEZ/yfLWUXhzYuxChTI8cbxcNPDuhDsVIV0o6LacLD9wEDsJMaeY3zG/8lyVlF482KsAlk43hoeBnxYF4KdqpBuVFQborYfTRLfZeNbb3sywbw0L8YqkIXjreFhvIeFI9ipfN5dy9XHp+NTBp6btJaCgSfxbTW+17YnbZPl7uxnsfsNQha700nujlaclnIewW5PYeTugaWIEewKLxh4IthVimBXuCO7U4KdQ0vjFMFuj9w9sBS53wecSO7Ofha532ckkrujFSf3GxJSU39MMrl7YClyvw84kdyd/Sxyv89IJHdHK07uNySkpv6YZHL3wFLkfh9wIrk7+1nkfp+RSO6OVpzcb0hITf0xAAAAZ0awAwAAaMTpgl3uw70oRe6eiL6nHnPL/f6fUe73HBHl7lyjUtqRTKRNX9FiK2pqdYvFVnnfiDOvvYQGnFOuzZ5lvfyxWZTSjmTqCgoEu7oWi63OHK2yd8LsDTgnsk6r6y2noEppRzJ1BQWCXV2LxVZnjlbZO2H2BpwTWafV9ZZTUKW0I5m6ggLBrq7FYqszR6vsnTB7A86JrNPqesspqFLakUxdQYFgV9disdWZo1X2Tpi9AedE1ml1veUUVCntSKauoECwq2ux2OrM0Sp7J8zegHMi67S63nIKqpR2JFNXUCDY1bVYbHXmaJW9E2ZvwDmRdVpdbzkFVUo7kqkrKBDs6lostjpztMreCbM34JzIOq2ut5yCKqUdAAAAOIhgBwAA0AiCHQAAQCMIdgAAAI0g2AEAADSCYAcAANCI1oPd43bpRte7/qzw3PHFCk92l9tjx2LvV9P8x1prXObxpi7LEVt1eMMaF3u8teYlBGotvKxt7VB9cqUJtvc6fmdwrD32327eDRqnoBYCWt2qyxsT9G13rncZEcKu2r/Aw3Zx7/UGr6y8+xOntoPd/Tpv1/t1fmPnZ4WXAyy2v18PvYmP20VclLiG3a21LPNoUyf3q9RxD29Y82IDtNa0hFCthY/1rR2oT7qbYHun43cG19oj/+2W/ZVpCmohoNWturwb4o46/noNk0Zf6fL043YJ2Ml81rts45B9O+/+ZEXTwe5+VQ7PTJt6enbfW21erPz0UfPCDrfWsMxATb1fu8tFaFKgpqqLDdBa0xICblisWt/aQctHNXx8vlyvljc6bmdYW3vcv922vzJPQS0Es7pVI73t/u9mwNSxtlLpyYB/eJ715t2f+Gg62EmmThz4nZ5r43G7BKvS5fNbkJG8MwAABlhJREFUuNYKnwmDNPVxu1zvYvPCNFVbbIDWmpYQa0cDk/WtHbJ8DOu/312f2SN3hpW1R/7bZaaxnFqIYXWrRhryvd/NkEcJvQpcOGycbL2Rgl3W/YmPswS7GFFJWux0oOnoyfzx1Lwx7+9trbLMEE0dx6LQwc6w2ACtNS2hgNo7kfWtHah8NjRjU/Nirj3N3z43wTnmUgvBrG7V+7W73O7BL8XyfTeD5kqvlU6XnQXsX147FvEUcNAoTbDLSzqpH/YYmLRYMTkd6z/iccBQXURa5sGm3q96Sg7QVNtij7XWuIQCau9EvD5Yhysfr2ZsaV70tcf/2x3XNlELMfgEO+m674D3MHi8m2GPF/oFrGXvHmrVPn/skifvgY+REuwyetwu2g2hATa6ulhJuKvCQnYR65GSrctcBp+gwc682MOtNS+hgNo7kW1bO0u0ynjEThRr7a79FbUQg98RO+e1j3HWG3h9fiuN1Mfy7lgIdrmYPhwEuLBx7TNHwMv9o9w8sfbk6nJUIe5KsS32YGttS8h/feuZbNraeU6GRu8M+YLdhv0VtRDMpuv6A253r3czdD8r9eYJaxtCyLs/cWo62FmuDr0fvBXZuFhxWfsuhrYs4VBrjcs83tSZ3GuPbljjYsvcsNhoZWsH7JMOznEuemewrT323+5xkTy1EMPqVl3G/6Db3ePdDB831lYqX+We7o+NtI37vs+9P3FqOdjph4DEsLD7ilXbYqdT+fsv0LQs4VBrjcs83lRhQcGaaltsmRsWGxm2tvjJPVifdNA+hwh5J35nsK896t9u21+l/vPPaKXPi+980O2+ut4YB4bXVrr8rWFXvbbe+fUoH5ny7U+cWg52AAAAp0KwAwAAaATBDgAAoBEEOwBIzvsb6yJ+tR1QJOGqUy773INgBwApbb4vkJtVcR5DqhPuZaXzb0awA1CvZDv++9W1HuneU+ftdzu/acLjS0uABqiFFvbrUU6CYAegXomC3RjcjOuRzhu5w92hQSrg704BhdJKxPkjTzAj2AGoV4pgtxyOM38DqvyS/bfOjw5RDHE4H/ehchgR7ADUK3KwU47GmX8zQ0tb41zy1MOTx1p6vxoTI9Am6Xo7eCPYAaiXMdjJccwyLAiXxV3vpiFE+lkAS4A05jpjsrNMuY0xMQJt4s6JvQh2AOql7fv1n9AyHGvTL4u7XK/GYCf+tLA+xlgPKGgvBMl1JDucBcfqjiDYAaiXHLjmwGZIevqvVs7PLGHQOpCYg531kIIawBw3X1hO0Vqedt3FAZTH8knL8pFLnIdUtxvBDkC9pGhlO54lpyHTVPY7HkzrWXlWf2VcvL50e6R0n+Ql2aESa8HOUot08CMIdgDqJQYo+3lKMSVZzvGsHCQ4FuzGwc28dPOLlltgVxMoUBTnhXL6IXZu/Q6BYAegXuKw4chmwmSWgWblmp5DwW4tjZkabvlGYoId6rJ2B4T8ccz4pZB0960IdgDqpQc74xCSOtiZR6uVpat3bjj+EkY6VGL11lbufQ2PYAegXnuP2FkuXtsY7Dzvil0/zKY26n61jXQEO1RlLbdxO1AEBLvTcF7DyjCBOuW9xs59+E+7KdZRZtLq7bGOU7GojDPYzaMSuS4ogt0prN6YxDCBOgW6K3Y1etmGJ79fnvC4l3WZ53G72OuRu2JRF4+7Yhl+QiPYnYD+RV72K3iAquT9HjthXvdvxXqcb5q+EuXiiHWcuEJtVoIdPTkGgl3zXEcx+KSEykX85Qn3etSXVg5CeHzlqs/BOH55ApUxVA7XE8RGsGveyjcyMESgZsF+K3bfNXaGhTlv31hPdu6C5A5CVMZ12xF9ORKCXfPuV9t4wjd8A5P4kcleiSkXASTlfxEDgiHYNc8R3xgmcD7ms5kpfp3y6Jfq86X8qI/HRQxEu9AIdu2zJruVb+4CmmQYTpJ9Odz9un8Y46pYVMh1KJwTspEQ7M7AfGSOXIdzyvqzRZafClt1JBIC2bivceA+iigIduegn37iKmycmRLu0o4r9+vGFW6eASiE52/F0sFDItidiPK1DIQ6IBf/b5LkUB1qtn4IgUvtgiPYAQAANIJgBwAA0AiCHQAAQCMIdgAAAI0g2AEAADSCYAcAANAIgh0AAEAjCHYAAACNINgBAAA0gmAHAADQCIIdAABAI/4HUYUxWjZDiqAAAAAASUVORK5CYII=)
where \(\nu\) is the normality (a.k.a. df) parameter of the \(t\) distribution.
Now for the issue of interest here: What is the predicted SAT score for a hypothetical state that spends, say, 9 thousand dollars per student and has 10% of the students take the exam? The answer, using the method described above, is shown below:
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAIAAAAAVb93AAAYpklEQVR4nO3dubHbyBaAYWaiDOSNDMbwvPGnZDEChaCq67IUhVw5zEGOEmAICoLPIJZG9zm9YD3d+L8aZ0gQRJPCf0Es5OUFADDpcvQCAABkBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0ABgFIEGAKMIdKMet0vgen/OmNPzfvXmc3usvrhreA+5aOHCsS1+uZZ43q/D8kcXbeZCvue5/7gwG4FujVTm2X2NzMxepdcO9N6Rni4/gcaLQLfF7WkYKneVz8lYP723Qg+zsZbo2YEWi6WMfkOJ5Z8xPB+Brg6BbkZmUbLD000oFaG7y9iavm6ghxnu9oeIQCNAoBsRxsTZnL49XpO1831XYkWNBcGP9zits5muzN/77O4/wTgr9/OANC//E8Hqgfbm+Lh10/qvrDAsbylSgwr2ZwijyBheMBvxw8/kRvlvbeI9wm4IdBP8NS/cd3x7TFbw/LU9a+XsZndPHU9Udqy6U3V/PK7hhJN5CSO87RLo221csv6B4rDcuaYGtUKg1cMFziP84fZPO5lpxnukvI46Cj8XgW7BsG33er2EnRjSVu37tty1LrWCjXUYnqC/SdimD5Lh3jTOyts+lyYKJtlgF0cwIm1zM/j7OGtQ83ZxdHObLllw43S44h6crPcoQKC3QqDrN22t+JlV2PH8vF/TuyOFzTLxIcJqHSyJtkXuRTKSGm8aeYhrBbofUpBdf2pxXF5Ncwa1JNDKPit/0ZzhRndtpN4j7IZA1+9xc9Yp5cCWsOZNHhYjbx4JOxyCFdhtir6OZ2zXTRdffraSXTLxoSl/jhY8bcagVj5IKJ7QM7zS0aCn3yPshkBXb7J/Qzv8F65gGfs4lOcSGq3Ew33WZAtj9ZUCHTzb2udBy4cv/de24IBrYnN1aaDTe6G94Wpb3PkvCbZGoKvnbgmrm3Ph2j3dbz3nOS/+dns00MkLaLIDrQ1y/bM4cuaf9aybB1q+3FPexXEZD3V6S5T7HmU8fd7jkEKgazfZl6wFR9mAjpQptWE4vT870OkWHrAFvUKgj92C1va/K4F2zrrM2lOVRqC3QqBrN9lVoWxcBh1Khym1a1XacRzdPZu7izg70Ovtg14U6IJ90JsFWgurv+s7tbd/zkuITRHo6gnHCBPnN2Sth/GzIrxZiGcpeNOolx9KG+Pxlu1xFkco+jEhfgDwmEAHJ6LIp8sHf7iT7xF2Q6CrNz3aJ+8ddm5SThYTqOcWq2eghbkRz4MOT40IHpdoWXDa4CbnQXu0QobnMAZ/LAoCnXnVjLgEwsnL0UALQc56j7AbAl0/73yM4EjP9X6/+bfkbQlln+TQxeMWPrV0ykN0XlktE4+KeVcSqluDwXwXBjrjRMSivzrKUudcqOK/+g/pQ4x0qbfwZzQyHOyGQDcgWMnGVb1bGZ2VrvxjarjGRo7PBVvs2tJqy5MZaG/Bxqc+ItDhsMQlLfirUxpobQGmTywOVzpHO3U+HvZCoJuQv9/Ce8hKq96MEyh28bgRF9SMQDdi2ObJqaS0o3ERm4F+3q/2FgooQKDbMflcGj35YvXPrSYD/bxf2XxG3Qh0W5IX6+ZvZpcwGWigegS6TVqoN2oogQa2QKABwCgCDQBGEWgAMIpAA4BRBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWg0ZzprxWs+bu46k8STH75nB/awloINJqi/JTM0h976QIszkZ8SiKNNRBoNGTYkB1L2udzQTHHzWMh0P2dw10rPCHQIdBoRp9GP6Pd7XO2or2t43AW7z57MV7whICLQKMZYiydewqDOcb59tDmoDwlhcY6CDSaoWdYT3fE834dHyPP+x1iYb7qHUAJAo1mpLagl2zRyoFW/ySwCY1VEGg0o6tiWOjgQF65wkDP26kCeAg0GtKX2G20c4oygUZlCDSaMrliZAzz4l4SaByBQKM5k0jfHq81jtmxDxpHINA4geXbs5zFgSMQaLRv1ll20iyU86CDm9mAxjoINNoR28xduDmrbYNzJSG2RKDRkG7vs9PL4XLAhbVUd5LwXRzYEIFGS5Qvs5O3cEsSGtmLLZ44Qp6xBgKN1iS/4GjlQI93E2esjEDjnB43QgrzCDTO6Hm/chAP9hFonNDzfmXzGRUg0ABgFIEGAKMINAAYRaABwCgCDQBGEWgAMIpAA4BRBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDSA47x/YT3x85DeD7Unf01y8iPrGdOXzn8/BBrAUfqSRorox7an/Kik39pkdQvnvzMCDeAQThq1QI+5HaeQbgtm6eR1nD5sbun8d0egAezO225NbN36d/cN9Yqr3Ow015tR4fwPQKAB7EnapyAGOpbJbibTu8Qb9VkVz/8ABBpYrlvX3+uztxvUjc80Tvrq7+9JjYUiL3iT4niz3zFDzsJe789hQcRARyMp3VkY6OL5H4BAA8sNgX6Ix6jcFuWmVJ5L5qTh5ENw7vJjdirRezGGQccCHY23eG/ZLo7y+R+AQAPLTerr9MErqHjoapIA6SiXc6sUGL9H43Oq2+7OHbFDaK6PP3+9/wpeHlUsg4ltWPFu6TXRDvrNmf/uCDSwnH7kfwxjsKqHCYj0KrgrY1r3GeVsq5MHdg90aqnUgMofK9StagINtC6ytuun+gZ1ijehuzfnM7cwo34xIkfE9FnrJxf7CvcIrB9odVG1CQk00LpIaNJHrvoHTY40Fj1HOJEa6PK9vW9VbEHLHxPkvTgEGjiLVQIdPeQ3mj5HbONWCLTYm0YCnfMhYbyPQANncUygpQdc78/YLo56Al18EC+xo0abnkADrVsz0JlJiJx+ETlIWFGgC0+DSw2i8AGcZgc0Y9Ut6KwmxD6htxHowgtJigPNhSrASawS6KJCRwIiniC9KNBHnMVRfCl26S4OLvUGTmKdQMsXpHhTdzPSn1H+kriyQIcby+IW9OJt6sSfhrIvMyo7SFg+/yMQaGC5lQKtXtcinD8mXnOof0ncToEuTHbufgnx60Bj16ikT7ObM//9EWhgudUCHd+dIF6rHBi/ESQ8B62yQEdGGT+2l/+IsvnvjUADy60Y6MldyV4oX/YRfkavNtDSy1E4+eo/qbUfAg3At1egkUCgAfgItBEEGoCPQBtBoIGzm1FeAr0PAg2cHYE2i0ADZ0egzSLQwNkRaLMINHB2BNosAg2cHYE2i0ADZ0egzSLQwNkRaLMINHB2BNosAg2cHYE2i0ADZ0egzSLQwNkRaLMINHB2mwaaZC9BoIGzI9BmEWjg7Ai0WQQaODsCbRaBBs6OQJtFoIFzWaWzBHofBBo4FwJdEQINnAuBrgiBBs6FQFeEQAPnQqArQqCBcyHQFSHQwLkQ6IoQaOBcCHRFCDRwLgS6IgQaOBcCXRECDZwLga4IgQbOhUBXhEAD50KgK0KggXMh0BUh0MC5HB5okp2PQAPnQqArQqCBcyHQFSHQwLkQ6IoQaOBcCHRFCDRwLgS6IgQaaNlGVSXQ+yDQQMsIdNUINNAyAl01Ag20jEBXjUADLasi0CRbQ6CBlhHoqhFooGUEumoEGmgZga4agQZaRqCrRqCBlhHoqhFooGUEumoEGmgZga4agQZaRqCrRqCBlhHoqhFooGUEumoEGmgZga4agQZaRqCrRqCBlhHoqhFooGUEumoEGmjHbg0l0Psg0IBNj9tl4vbIn3bq359+DX99/+TM+KtWzB//XS6Xy+fvz7JA//7f5+559Wn+fvz5+/Xfy+VyuVzv38YJ+sfqMl6H6/255jtxIAINmKP3Vq7T834tCPQ7uxP//O9X2NCfX/q7XpYCHUswgQawKae2To7HW8NG9/fdHhkNfXYF/Pz925+/H39+d4/t/nd81Ldv/wxxD+ezeaCV7WTnxZEqTKABbGgokBCZfsPaj1d3++2RsXv3ef803WQWbnm93M1nU4GevEThRAQawHb6+MQ/wk/D1D3men9mBLqbw39fx2keXy6Xy+Xy5cf4KHfz2Vygx9cheJkINIDtJAoz7swIHnO9h4fyXsEt/QyEQH/69nu6Te1OYyzQ2h8yAg1gM4kNaLFAY5+9Y4vS2Rc5W9DvSYZeHxNo58Dmh3TinfxKEWgA28ncgnY2MBNncLid/cjZB/3r+3VacJuBlnfIE2gA28nbB+1madxsvt6fbtS6/cjOpvFHxlkc72j6WTcYaPGVItAANiRsIw+cXRjDvdPDZV4NpT3O0fOg33cNxXSndM7Dyw10JgKtI9CALZNNYulmJd+vl3pS3eQMjY8/f/uz6Lp5fZ3G8X263rAB7kjs93gR6LURaMAc7UrC22Nyol3GPoen17vYoyYXdv90T+0YDiQKV40fvQ+aQAPYnXfwr2uOc87zyoHuyvjlh7sbpNtkDm+xE2gOEgIwYnKlyoqB9i7s7v633+/8er3671fyLjg8LNDypSoEGsBhpgEaG/f5+ze5mMFFKHIx/Qu73dmaDLRyKSGBBrChaGK8Oz+GHcfKd84p37PhFzO8sNv6FjSXegM4QORE6GCv68efv+P5GFINu0fEj+yJF3Zb3gfNlyUBOMpwCoebGTFK73gN58N5F6p0+UudG6dc2G30LA7n/Ba+bhTAAfSrt6ebjEO/pHOWu4q538QvVFW/sHuF86D5wv7FCDRgUsYvXk3b99N7wKdvv4Np/KomLuxeeCXhJoEe97roLxqBBnC0eHw3vWXnpxNvOQMCDdSKQDePQAO1ItDNI9BArQh08wg0UCsC3TwCDaS8T3vL/JW8yFkXa09/9kBnvS/BS1j6QifOCSl9H8sQaCBO/No0eRKftmqvNf301LfXuQKd8b6Uvo+Rnw8TH1X6PpYj0ECE9CMmHukiv8jVyCtPn3VpX3uB/pl+X0rfR/kLPsZX2q9u6fs4C4EGFN72kbLKaRtm2pdqrD19zrcXtRVo7ye7kilc9j5qd5W+j/MQaCAkfXYVV+zY6ihd1jZ3eu9iv48/f7O/SrSlQLu/1JUT6OL3UZ5ZGOPS93EuAg14vC/kia67Jd8OumR67ys13s0SvrW55UA7df78PfG+FL6PyXllvFWZd5Yh0IDnvX4N6+r8FTu8d/b0X6U+9t9n5HwpaPuBHgabGeic97G4z6Xv42wEGoiLrW6JbaXg7tnTfxP76Hxb/wkC7d3i/5rXR+Lk6Mj76L3u6Z3Wpe/jbAQaiMvY8spdUedPT6B3CPT4o+mB6TtW+j7OR6CBOALdeqCHl/kqTDGcOee8aQQasIJAnyXQyrscvm0EGrDigED352z0315PoHcKtHZYz3sjCTRgxQEHCYeT6qTfHCHQawda/ZFw9bEcJASMOOA0u+Gkuv5EOk6z2zTQyav//Mdymh1gxAEXqozby902MheqbBro5CZv99jMj0JcqALsJhrobS71Hkv06/uny+XCpd7ZgQ4fmPM+DnuhSz8Kcak3cKzEB9blX370Tkm/v+Kf+3NSmW4zud+PMfx33i9L2iDQ8d3QfFkSYFZqj2L0ayeF1TSY3qnz5dO33y+vMr++99P/F+ybnmxZ+w8k0AWBdq4fnL5nw81aiXPf91kINBCXcchHu/4s2QLPvz/lykSnJ9DrBDryviR2Nme/7+UINBCXeUx+6U9YffkRz1N/TrQz/wP72GagX69XWN3ExjA/eQU0pI0+Hr4ASxapIgQa2FUbfTx8AZYsUkUINLCrNvp4+AIsWaSKEGhgV2308fAFWLJIFSHQwK7a6OPhC7DuIplFoIENGYxRGwuw7iKZRaCBDRmMURsLsPUiGUGggQ1VEaMaF2DrRTKCQAOrWaUO2sVsKNVAsgk0MNNGm29HZ60d+2/Ur45AA7YcnbV2HP1OrqCFMQAtOTpr7Tj6nVxBC2MAWnJ01tpx9Du5ghbGAABNItAAYBSBbtbRny+B0dFrQ6144ZrV6lrBuKrT8NC2xgvXrFbXCsZVnYaHtjVeuGa1ulYwruo0PLSt8cI1q9W1gnFVp+GhbY0XrlmtrhWMqzoND21rvHDNanWtYFzVaXhoW+OFa1arawXjqk7DQ9saL1yzWl0rGFd1Gh7a1njhmtXqWsG4qtPw0LbGC9esVtcKxlWdhoe2NV64ZrW6VjCu6jQ8tK3xwjWr1bWCcVWn4aFtjRcOAIwi0ABgFIEGAKMINAAYRaABwCgCDQBGEWgAMIpAA4BRBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEujKP22V0e4QTPO/XhRMcIb1QtY3rcdtsIEePVB1acgLzQzOHQFdk8q+3c70/nSkm+Zb+lScnOIKwUKWLbWxc3TsVLsLygRw9UnVoyQnMD80iAl2P7t/v8G+2WxPGRPsT+P+fnuAI/Wq5YLFtjWv8O+ovwPKBHDxSfWjJCawPzSgCXY3H7eJvML9Xh/6f8PIJjvBeD6dLUO+4+r8216vw/HW/g9GhJScwPTTDCHQ1hH+uk5uEf+HTKZITHEF+/qLFtjOu8Umlp18+kANHGh9acgLLQzONQNfD37k3/X/lH/P4Dz85wSGSga5zXKm/pqOCgdgYabKam4z9rAh0VYLDKOM/6VrXAen53b3rdY6LQBPoVRDomsQOc1e7DnSDGhdhONBUU7Y8BJpAr4JAVyM8qj25pd51IDx78Hq/15YtD4Em0Ksg0LWIHEWpafVWOJ8Nbo/+/9kHXTjBLgj0ngh0JWLHzt//gls6UO4uapXjau4sjoLn4yyOFRHoSjjbytqtdZ5q6vyJGUwXtMZxic9f93nQ2U+40djPiUDXot9TG+6DnpbMmaKOi7XeCxFcD+msqxWOS25L9VcSvl5zA13F0Awi0PUID6Z5IROnEPYFRiY4grRQ0qZUTeOKXs2xZCDHj3RmoGsYmkEEui7Tf8TS4RP3VDzx8EpygiOkF6qucUUqtnwgx450bqBfL/NDM4hAA4BRBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0ABgFIEGAKMINAAYRaABwCgCDQBGEWgAMIpAA4BRBBoAjCLQUD3v10vS9f48ejlnq2SA/WIevyTYG4GGKqtfb7fHzsv2uK3wrJYHKCwmgT4fAg1VQb92rcewXDsG+sg8EujzItBQZYTBTdxu/Vg90NYGqCwDgT4fAg1VXhjGhO0VkD0D/TpigMoCEOjzIdBQ5YbhcbusE8zC5dop0PsP0EOgz4tAQ7VKv4Y7Mxon7BP2JvfntixcCwbY3XR7vKK7QYIRJRI/nb6b2eJA98uvzoA/AVYRaKiWBlqsqZYpfWKhi4YC7SV4nEI/Aik/nTL99f5cns9EoemzWQQaqtJdtJPsxoIbTCxug7pz6Kc+ItDiAPtA36YL1E+ROj8keMbIy3W93Rb3MzpS+mwXgYYqveJOOuQmV4qrfnvkifoHuDPZcR901gCDh6v7PJw7Jrc7c9NuXxjQyFDps2EEGqqi04TFgkoJFU6KSO4jVeaw63nQ0ycb0xkuRHT4kf0liZdrWUHVF62/48hLcaAh0FDl98tLR+YuT3mHQDoUBwRa3yUR6bOyeMEWa+LlGp5q4Sau9GGEPhtHoKHK6ZcUjeRn5nACeRdsonC7BDpRzcy9MuLTdhMkX67CDxgqacGSC4tDEWio5u6czA60m4VYLL167H4etCAj0Lntzw764p3E4TPRZ+MINFSbB7qgbu6ktgNdvN9kv0AH+zO6p+bwoFkEGqojAu3wWj1ObDvQxVulK71cRc/1Xjb6bB6Bhmp2GEoPEmYtxTi18UAXL17uy7VGSN1C02f7CDRUs/sV/Xqh8M7481QX6NRpF+HiRx8wfpBYpaRjoelzBQg0VCv0y8+ofLtzq99c8XzjtU5r2CjQsS/nkHOb83KtltJNZoptEGioFu37LLrUO33mQ+QCu/G+0p0BWwU650hh9veWaJd6z9/zof05gD0EGqqFB6cikYpfMZee3stZcMHL4YFODEjMonouuLyYC3ZNxy6DhCkEGqo1zh4o+77N9PeNilP6O7NNBFoZUCKJk0pHv250ybFDTn+uBYEGzoY+V4NAA+fCt9dVhEADZ7LqWdXYGoEGTiA4AsnujSoQaOAElPNeYByBBk5gre/9x74INAAYRaABwCgCDQBGEWgAMIpAA4BRBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0ABgFIEGAKMINAAY9X9SLBqg2CIRuQAAAABJRU5ErkJggg==)
(Please note that the numerical annotation in these figures only shows the first three significant digits, so you'll need to examine the actual MCMC chain for more digits!) As another example, what is the predicted SAT score for a hypothetical state that spends,
say, 9 thousand dollars per student and has 80% of the students take the
exam? Answer:
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAIAAAAAVb93AAAZN0lEQVR4nO3dvZHbyBZAYWaiDNbbNRjD89ZXyWIEE4Kq5LI2CrnrMAc5SmBC2CD4DAJEo/ve/gEawO3G+WqdJUEQzRHPYBoAeXkCAEy6HL0BAAAZgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0ABgFIEGAKMINAAYRaABwCgC3anH7RK43j8XrOnzfvXWc3tU39waXkMu2rhwbKtfrjU+79f39kc3beFGvta5/7iwGIHujVTmxX2NrMxepWsHeu9Iz7efQONJoPvi9jQMlfuWz8nYuLz3hn6vxlqiFwdaLJYy+g0ltn/B8HwEujkEuhuZRckOz7CgVIThLmPv9LqBfq9wt19EBBoBAt2JMCbO7vTt8Zy9O193Jd6osSD48Z6WdXbTlfV7f7v7TzCtyv17QFqX/xdB9UB7a3zchmX9V1YYlrcVqUEF8xnCKDKGF6xG/ONndqP8uzbxM8JuCHQX/HdeOHd8e8ze4Pnv9qw357C6e+p4ojKx6i41/PK4hgvO1iWM8LZLoG+3acvGB4rDcteaGlSFQKuHC5xH+MMdn3a20oyfkfI66ij8UgS6B+99u+fzKUxiSHu1r9ty33WpN9hUh/cTjDcJ+/RBMtybplV5++fSQsEiG0xxBCPSdjeD34+LBrVsimNY23zLghvnwxVncLJ+RgECvRUC3b55a8W/WYWJ58/7NT0dKeyWiQ8R3tbBlmh75F4kI6nxlpGHWCvQ45CC7PpLi+PyapozqDWBVuas/E1zhhud2kj9jLAbAt2+x815TykHtoR33uxhMfLukTDhELyB3abo7/GM/br55svPVjIlEx+a8utoxdNmDKryQULxhJ73Kx0NevpnhN0Q6ObN5je0w3/hGyxjjkN5LqHRSjzcZ022MFZfKdDBs9U+D1o+fOm/tgUHXBO7q2sDnZ6F9oar7XHnvyTYGoFunrsnrO7Ohe/u+bz1kue8+Pvt0UAnL6DJDrQ2yPpnceSsP+tZNw+0fLmnPMVxmQ51eluU+zPKePq8xyGFQLduNpesBUfZgY6UKbVjOL8/O9DpFh6wB10h0MfuQWvz70qgnbMus2aq0gj0Vgh062ZTFcrOZdChdJhSU6vSxHF0ejZ3ijg70PXmoFcFumAOerNAa2H1p75Ts/1LXkJsikA3TzhGmDi/Iet9GD8rwluFeJaCt4x6+aG0Mx5v2R5ncYSifybEDwAeE+jgRBT5dPngF3fyZ4TdEOjmzY/2ybPDzk3KyWIC9dxi9Qy0MDfiedDhqRHB4xItC04b3OQ8aI9WyPAcxuCXRUGgM6+aEbdAOHk5GmghyFk/I+yGQLfPOx8jONJzvd9v/i15e0LZJzkM8biFTy2d8hBdV1bLxKNi3pWE6t5gsN6Vgc44EbHot46y1TkXqviv/kP6I0a61Fv4NRoZDnZDoDsQvMmmt/rwZnTedOV/pobv2MjxuWCPXdtabXsyA+1t2PTURwQ6HJa4pQW/dUoDrW3A/InF4UrnaKfOx8NeCHQX8uctvIdUeustOIFiF48bcUHLCHQn3vs8OZWUJhpXsRnoz/vV3kYBBQh0P2Z/l0ZPvqj+d6vJQH/er+w+o20Eui/Ji3Xzd7NLmAw00DwC3Sct1Bs1lEADWyDQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0OjO/NsK1n3t1fwrzdWvJJgtxhdtoRYCja4oXyWz6Mte5nFWVyU+JZFGDQQaHXkndaromM/iYo7rmlYV3iLduPgJgQCBRjfGNPr7uMPtZXvRr/B6lQ3XlLkYsAiBRjfEWDr3FAVTfoj/FMpTUmjUQaDRDT3DerpVr8Ymdo3FhaJ3ACUINLqR2oMu26MNJ0yGtThPoP5KYBcaVRBodGOoYlho8ehe/gr1szNS++wEGusQaHRkLLGbUedkucJeyqfsuSsh0NgWgUZXxHOXb4/yXkp73cFtBBrbItDozizSt8dzwTE7dQ55Ps/NHDS2RaBxAsX7s5knhHAWB7ZFoNG/8rPscs/YU5ZjBxp1EGj0Q86lenJHROqEEG8WmisJsQ0CjY6EZyq/z8RYegaHdEKIdBufxYENEGj0RPkwO3kPN5FQ+cPstMu/488ILEKg0Rsv0vo5FhkR9YOv7YfzedDYBIHGOT1uhBTmEWic0ef9ykE82EegcUKf9yu7z2gAgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0ABgFIEGAKMINAAYRaDRpdnXbGtfx+0vpBkfXLo8sBKBRm8+79fMcBJoGEeg0Zd3RJ1vhZ2SXfBVseOKch9SujyQRqDRkzHFQSXfjc7cux1zm7szXLo8kINAoydDJ4W9WDXdkdWQZxyMQDdkaMwrMd5Mq5uG+VSpXiR/sjbWLmn6NczRrI/e6nf5079KoJncgBUEuiHvQD/Ew2DX+6d2gExNqbyWzEXDxd99vMuP2TxhFaY4yDPsINANmdXX6YFXUPHo2KxL0oE051Z32fcavP5Mz6nuuzt3TNvx5ePXj9//vf+r8bLMSRusvAzRh2fOVpQuD5Qg0A3ROzOFMdiPC//oH1cjJCW4K2NZ9xnlbM8W/+P7xzzQG/Ra3ufP2cVl9xmmEOiGRKZR9aNUQWL1aVrn3pz9QWFFsV6N9339pgT64+NPqauC+Nap50GnK8rBQdhCoBsS2Z1NHxwbHzQ70lj0HOFCaqClR48PUgNdZQ9a3ofPOxGa6Q1YQ6AbUiXQmVfDzZ9DvzhPDLRYwT0CnbMHrya66ES8BcsDxQh0Q44JtPSA6/0zNsVxVKATMw7xyR36DHsIdENqBjozK5G5gchBwoMCnZqeid9f9sKULw+UI9ANqboHnTVvGttL3CjQ7v8WJntVoOkzDCLQDakS6KJCRyokniC9KtDrz+JYMcWRdXB0xfLAEgS6IXUCLV+Q4i09rEh/Rmdm2rnz2D3oNQcJ2YGGRQS6IZUCrV7XIpyiJl5z6B02jATaTe24g7xloBefZscONEwi0A2pFuj4eXPzJbWzPqZPBAkvOj8w0CVDy3r96iwPLEKgG1Ix0LO7kgFTPuwjPExoIdBlQ0u+SFWWB5Yh0NhQWNv4LeICwGkRaGyIQANrEGhsiEADaxBobIhAA2sQaGyIQANrEGhsiEADaxBobIhAA2sQaGyIQANrEGhsiEADaxBobIhAA2sQaGyIQANrEGjUVJTjzECTbJwWgUZNBBqoiECjJgINVESgUROBBioi0KiJQAMVEWjURKCBigg0aiLQQEUEGjURaKAiAo2aCDRQEYFGTQQaqIhAo6Z9Ak2ycRIEGjURaKAiAo2aCDRQEYFGTQQaqIhAoyYCDVREoFETgQYqItCoiUADFRFoLFe9tgseQqDRMQKN5Qg0sCkCjeUINLApAo3lCDSwKQKN5Qg0sCkCjeUINLApAo3lCDSwKQKN5Qg0sCkCjeUINLApAo3lCDSwKQKN5Qg0sCkCjeUINLApAo3lCDSwKQKN5Qg0sCkCjeXsBJpeo0sEGssRaGBTBBrLEWhgUwQayxFoYFMEGssRaGBTBBrLEWhgUwQauXao7YKHEGh0jEAjF4EGdkagkYtAAzsj0MhFoIGdEWjkItDAzgg0chFoYGcEGrkINLAzAo1cBBrYGYFGLgIN7IxAIxeBBnZGoJGLQAM7I9DIRaCBnRHo83jcLjO3R/6yc3//9LP47/cv091fv2klfa32j+8fZYH+9b8/lOedP+rb35fL5XK53j+nkXzer7Gh5L0OszUCuyHQp6D3Vq5TompeKP/5Gizx5//+DUv68zbeVbgHvW2gYwkm0DgUge6fkygnx9OtYaPH++b3KFkc6/nH9083psP/Tst/fPz5juwBgVb2k50XR6owgcahCHTv3gUSIjPuWPvxGm73bpazOExuOLvM4y33T/chP/+6XC6Xq7v7bCHQs5coXIhA41AEunNjfOJ/ws/DNDzGf4icxWF+4+u3KZSvFl9uj+khw+7z7aG19chAT69DMGYCjUMR6M4lCiPmKziUF8niUF4p0Nf754/ZIcSvj4VncRQH2tny8bHxQGu/yAg0DkWg+5bYgRYL9Lrpy8cv/+ifmOyMPehXN798/HraDbTyShFoHIpAdy5zD9rpV+K0hy8fv8rmoMfd52/Rth4eaHlCnkDjUAS6c3lz0G6WplPy3i1+RXCczbj89Y+bxcRZHO/d53hbjw+0+EoRaByKQPdO2Ed+c06Pft873ubWOTrjHD0P+nXXODcyOx/7j+8fpYHORKDRCwLdv6mKUnv8QI/kvdTxisG//vFiOsw7v/P94/d/z+fnK46vPe73DrjDP2yoPC+BxkkR6FPQriS8PZRznpOXpfz9U19m/N/Z2SBDwccd8zHo8xPvlHXuOgdNoGEHgT4N7+Df0BzlnOcKgf71WvNr93mcBrl9ey8Q3mIj0BwkhB0E+uS0Hei1gfYu7B7+103nMFtyjXw0x56Bli9VIdA4FIE+tyBA7xs+5AhOMxXPWKD9C7uFdNoKtHIpIYHGoQh076KJUS9T0T40IzzrWQpleGH38XvQwYc0SS8El3rDFgLdu8iJ0PJnJY23Sh87NxQwHkrxwu7D56AjgebDkmAVge7f+xQONzORKL3v8i5UGfIX+Tz++ZUp8wWOPotDCbRzfgsfNwpzCPQZ6Fdvpz8l2TN93P5TDKV+YXeF86DXBDqBD+yHRQT6NEq+8er5fM6vPblclGsLvVsSF3avvJJwk0BHXwcCjUMRaAwyQll2yz4PWbwSwD4CjQGBBqwh0BgQaMAaAo0BgQasIdAYEGjAGgJtUOH5FsFJcbkfOzFf/qyBLnm1tU8FzHvdhx9UcjFgRKBtKT5jWUuGcmJYZP0ZV/T1Fmi9t/Krnfg2sMyP+yDQyEagLZGu+ZuiEDZXuhwwcolgYv3Od8I+TxBo58KZr9MrFXv5xK9AL/3REmgUINB2qJ+aoX1rlbZLpqwovf7hs5vPEOjP+xfn19L89dBeV/WjWRO8HXUCjWwE2g79qjUxrbGvgxVXlV6//43d/QZ6HPFw5bry+s1Tqn63QYQ0iUKgkY1A21EY6NLPESXQzn/Da/H+NlvlBRG+XaWgz94HMa2ZIcFJEWg7yqY44m936d7Z+t1avWdjTzPFMXwdlxpo/YOyr/dPf7dYTfZrufePgECjGIG25F1i8SCe+LHNWh3Eu531v48HTsfKnA/kfHYe6Nw9aPHgqyhjx5pAoxiBtkY+9Uvdqy4KdHT9h4Ry/UOWrSRvDtqN6fS6Kcdes09XJ9DIRqBtUffT/NIuDHRk/R9HhHL9Qxau5H0Wx/gRpuFrN4up9p1Ys1c1nl4CjWIE2pApDNLUsfyNKCWBdtc/1Wr4iP1L+sP4ewr09AUxl8sf38UXqSCmWe0l0ChGoM2I7aOF95UHer6OebB+3sZUfeweyqMC/eP9Af+B26P0nOec+BJoFCPQViQuBA6CW3SQ0ImR8GVUP36/v+tkmpM9Q6B//Hb/gHB/Rf0qPOeZQGMTBNqI1Ls3uL/oNLv0t0aNc7KvM+2EBXoNtLzAz8KLBgk0NkGgjSgOdNGFKgS67CH/fC37eyXruwsJNIoRaCtKpziKLvX+wRRHyXjfXz7rX6aS+JDA+IwIgUYxAm1G2UFC99b0hyU5e4WXLx+/nhwkHP9iEL6AfHyhxAuDxLzmHVEk0ChGoA0pOc3Ou0u44s1derZjeOE0u9lpdu5vrOm6Svcl+v3fj+naFv/HkP0xzwQaxQi0KbHrieX3tfaZ8/J12+NMtOT9QRwnCXTs1fj7p/gQ/ceTc7oHgUYxAm1PmIHMv52l1Ap5+ve7v/6gR/5D+gz0fz9+TxMao6/f9Ic8n8/wV2LVU/GAGQLdoUPa2mqgSx4C7IxAd4hAE2j0gUB3iEATaPSBQPfAQuPOEGiSjZ0R6B5YaByBBqoj0D2w0DgCDVRHoNtjs3HnDDS9xqYIdHtsNs5UKO08L7AGgbauYmuUq+AQQ6BxIAJtznY7g0e3rklb77kDEQT6RI5uXZOO/qHh1Pj3dyJHt65JR//QcGr8+zuRo1vXpKN/aDg1/v0BgFEEGgCMItDdOnpuAJgc/W5oFS9ct3p9VzCu5nQ8tK3xwnWr13cF42pOx0PbGi9ct3p9VzCu5nQ8tK3xwnWr13cF42pOx0PbGi9ct3p9VzCu5nQ8tK3xwnWr13cF42pOx0PbGi9ct3p9VzCu5nQ8tK3xwnWr13cF42pOx0PbGi9ct3p9VzCu5nQ8tK3xwnWr13cF42pOx0PbGi9ct3p9VzCu5nQ8tK3xwnWr13cF42pOx0PbGi8cABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0ABgFIEGAKMIdGMet8vk9ggX+LxfVy5whPRGtTaux22zgRw9UnVoyQXMD80cAt2Q2b/ewfX+6Swxy7f0rzy5wBGEjSrdbGPjGn5S4SasH8jRI1WHllzA/NAsItDtGP79vv/NDu+EKdH+Av7/pxc4wvi2XLHZtsY1/R71N2D9QA4eqT605ALWh2YUgW7G43bxd5hfb4fxn/D6BY7weh/Ot6DdcY2/ba5X4fnb/glGh5ZcwPTQDCPQzRD+uc5uEv6Fz5dILnAE+fmLNtvOuKYnlZ5+/UAOHGl8aMkFLA/NNALdDn9yb/7/yj/m6R9+coFDJAPd5rhSv00nBQOxMdJkNTcZ+1kR6KYEh1Gmf9Ktvgek53dn19scF4Em0FUQ6JbEDnM3+x4YBjVtwvtAU0vZ8hBoAl0FgW5GeFR7dku774Hw7MHr/d5atjwEmkBXQaBbETmK0tLbW+H8bXB7jP/PHHThArsg0Hsi0I2IHTt//Qvu6UC5u6lNjqu7szgKno+zOCoi0I1w9pW1W9s81dT5FfM239AWxyU+f9vnQWc/4UZjPycC3Ypxpjacg56XzFmijYu1XhsRXA/pvFcbHJfcluavJHw+lwa6iaEZRKDbER5M80ImLiHMBUYWOIK0UdKuVEvjil7NsWYgx490YaBbGJpBBLot83/E0uET91Q88fBKcoEjpDeqrXFFKrZ+IMeOdGmgn0/zQzOIQAOAUQQaAIwi0ABgFIEGAKMINAAYRaABwCgCDQBGEWgAMIpAA4BRBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0FB93q+XpOv98+jtXKyRAY6befyWYG8EGqqsfr3cHjtv2+NW4VktD1DYTAJ9PgQaqoJ+7VqP93btGOgj80igz4tAQ5URBjdxu/WjeqCtDVDZBgJ9PgQaqrwwTAnbKyB7Bvp5xACVDSDQ50OgocoNw+N2qRPMwu3aKdD7D9BDoM+LQENVpV/vOzMaJ8wJe4v7a1sXrhUDHG66PZ7RaZBgRInEz5cfVrY60OP2qyvgV4BVBBqqtYEWa6plSl9Y6KKhQHsJnpbQj0DKT6csf71/rs9notD02SwCDVXpFO0su7HgBguL+6DuGsaljwi0OMAx0Lf5Bo1LpM4PCZ4x8nJdb7fV/YyOlD7bRaChSr9xZx1ykyvFVb898kTjA9yV7DgHnTXA4OHqnIdzx+x2Z23a7SsDGhkqfTaMQENVdJqwWFApocJJEck5UmUNu54HPX+yKZ3hRkSHH5kvSbxc6wqqvmjjHUdeigMNgYYqv19eOjKnPOUJgXQoDgi0PiUR6bOyecEea+Llej/Vyl1c6Y8R+mwcgYYqp19SNJJ/M4cLyFOwicLtEuhENTNnZcSnHRZIvlyFf2CopA1LbiwORaChWjo5mR1oNwuxWHr12P08aEFGoHPbnx301ZPE4TPRZ+MINFSbB7qgbu6itgNdPG+yX6CD+YzhqTk8aBaBhuqIQDu8Vk8L2w508V5ppZer6Lle20afzSPQUC0OQ+lBwqytmJY2Hujizct9uWqE1C00fbaPQEO1uF/RjxcK74w/T3OBTp12EW5+9AHTHxJVSjoVmj43gEBDVaFffkbl251b/eaK5xvXOq1ho0DHPpxDzm3Oy1UtpZusFNsg0FCtmvssutQ7feZD5AK76b7SyYCtAp1zpDD7c0u0S72Xz3xovw5gD4GGauXBqUik4lfMpZf3chZc8HJ4oBMDErOongsub+aKqenYZZAwhUBDVePsgbLP20x/3qi4pD+ZbSLQyoASSZxVOvpxo2uOHXL6cysINHA29LkZBBo4Fz69riEEGjiTqmdVY2sEGjiB4Agk0xtNINDACSjnvcA4Ag2cQK3P/ce+CDQAGEWgAcAoAg0ARhFoADCKQAOAUQQaAIwi0ABgFIEGAKMINAAYRaABwCgCDQBGEWgAMIpAA4BRBBoAjCLQAGAUgQYAowg0ABhFoAHAKAINAEYRaAAwikADgFEEGgCMItAAYBSBBgCjCDQAGPV/naA56DKcMqkAAAAASUVORK5CYII=)
A couple more examples of predictions:
Now for
the R code I used to generate those results. I modified the scripts supplied with DBDA2E, named
Jags-Ymet-XmetMulti-Mrobust.R and
Jags-Ymet-XmetMulti-Mrobust-Example.R. Some of the key changes are highlighted below.
The to-be-probed \(x\) values are put in a matrix called
xProbe., which has columns that correspond to the \(x\) predictors, and one row for each probe. The number of probed points (i.e., the number of rows of
xProbe), is denoted
Nprobe. The number of predictors (i.e., the number of columns of
xProbe), is denoted
Nx. Then, in the new low-level script, called
Jags-Ymet-XmetMulti-MrobustPredict.R, the Jags model specification looks like this:
# Standardize the data:
data {
ym <- mean(y)
ysd <- sd(y)
for ( i in 1:Ntotal ) {
zy[i] <- ( y[i] - ym ) / ysd
}
for ( j in 1:Nx ) {
xm[j] <- mean(x[,j])
xsd[j] <- sd(x[,j])
for ( i in 1:Ntotal ) {
zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
}
# standardize the probe values:
for ( i in 1:Nprobe ) {
zxProbe[i,j] <- ( xProbe[i,j] - xm[j] ) / xsd[j]
}
}
}
# Specify the model for standardized data:
model {
for ( i in 1:Ntotal ) {
zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigma^2 , nu )
}
# Priors vague on standardized scale:
zbeta0 ~ dnorm( 0 , 1/2^2 )
for ( j in 1:Nx ) {
zbeta[j] ~ dnorm( 0 , 1/2^2 )
}
zsigma ~ dunif( 1.0E-5 , 1.0E+1 )
nu ~ dexp(1/30.0)
# Transform to original scale:
beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd
beta0 <- zbeta0*ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd
sigma <- zsigma*ysd
# Predicted y values at xProbe:
for ( i in 1:Nprobe ) {
zyP[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zxProbe[i,1:Nx] ) ,
1/zsigma^2 , nu )
yP[i] <- zyP[i] * ysd + ym
}
}
The changes noted above are analogous to those I used for simple linear regression in
the previous post. The MCMC chain monitors the values of
yP[i], and subsequently we can examine the posterior distribution.
I hope this helps. Here is complete R code for the high-level and low-level scripts:
High-level script:
#===== Begin high-level script ================
# Example for Jags-Ymet-XmetMulti-MrobustPredict.R
#-------------------------------------------------------------------------------
# Optional generic preliminaries:
graphics.off() # This closes all of R's graphics windows.
rm(list=ls()) # Careful! This clears all of R's memory!
#.............................................................................
# # Two predictors:
myData = read.csv( file="Guber1999data.csv" )
yName = "SATT"
xName = c("Spend","PrcntTake")
xProbe = matrix( c( 4 , 10 , # Spend, PrcntTake
9 , 10 ,
4 , 80 ,
9 , 80 ) , nrow=4 , byrow=TRUE )
fileNameRoot = "Guber1999data-Predict-"
numSavedSteps=15000 ; thinSteps=5
graphFileType = "png"
#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Jags-Ymet-XmetMulti-MrobustPredict.R")
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
mcmcCoda = genMCMC( data=myData , xName=xName , yName=yName , xProbe=xProbe ,
numSavedSteps=numSavedSteps , thinSteps=thinSteps ,
saveName=fileNameRoot )
#-------------------------------------------------------------------------------
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda) # get all parameter names
for ( parName in parameterNames ) {
diagMCMC( codaObject=mcmcCoda , parName=parName ,
saveName=fileNameRoot , saveType=graphFileType )
}
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda , saveName=fileNameRoot )
show(summaryInfo)
# Display posterior information:
plotMCMC( mcmcCoda , data=myData , xName=xName , yName=yName ,
pairsPlot=TRUE , showCurve=FALSE ,
saveName=fileNameRoot , saveType=graphFileType )
#-------------------------------------------------------------------------------
# Plot posterior predicted y at xProbe:
mcmcMat = as.matrix(mcmcCoda)
xPcols = grep( "xProbe" , colnames(mcmcMat) , value=FALSE )
yPcols = grep( "yP" , colnames(mcmcMat) , value=FALSE )
xLim = quantile( mcmcMat[,yPcols] , probs=c(0.005,0.995) )
for ( i in 1:length(yPcols) ) {
openGraph(width=4,height=3)
xNameText = paste( "@" , paste( xName , collapse=", " ) , "=" )
xProbeValText = paste(mcmcMat[1,xPcols[seq(i,
by=length(yPcols),
length=length(xName))]],
collapse=", ")
plotPost( mcmcMat[,yPcols[i]] , xlab="Post. Pred. y" , xlim=xLim ,
cenTend="mean" ,
main= bquote(atop(.(xNameText),.(xProbeValText))) )
}
#-------------------------------------------------------------------------------
#===== End high-level script ================
Low-level script, named
Jags-Ymet-XmetMulti-MrobustPredict.R
and called by high-level script:
#===== Begin low-level script ================
# Jags-Ymet-XmetMulti-MrobustPredict.R
# Accompanies the book:
# Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
# A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
source("DBDA2E-utilities.R")
#===============================================================================
genMCMC = function( data , xName="x" , yName="y" , xProbe=NULL ,
numSavedSteps=10000 , thinSteps=1 , saveName=NULL ,
runjagsMethod=runjagsMethodDefault ,
nChains=nChainsDefault ) {
require(runjags)
#-----------------------------------------------------------------------------
# THE DATA.
y = data[,yName]
x = as.matrix(data[,xName],ncol=length(xName))
# Do some checking that data make sense:
if ( any( !is.finite(y) ) ) { stop("All y values must be finite.") }
if ( any( !is.finite(x) ) ) { stop("All x values must be finite.") }
cat("\nCORRELATION MATRIX OF PREDICTORS:\n ")
show( round(cor(x),3) )
cat("\n")
flush.console()
Nx = ncol(x) # number of x predictors
Ntotal = nrow(x) # number of data points
# Check the probe values:
if ( !is.null(xProbe) ) {
if ( any( !is.finite(xProbe) ) ) {
stop("All xProbe values must be finite.") }
if ( ncol(xProbe) != Nx ) {
stop("xProbe must have same number of columns as x.") }
} else { # fill in placeholder so JAGS doesn't balk
xProbe = matrix( 0 , ncol=Nx , nrow=3 )
for ( xIdx in 1:Nx ) {
xProbe[,xIdx] = quantile(x[,xIdx],probs=c(0.0,0.5,1.0))
}
}
# Specify the data in a list, for later shipment to JAGS:
dataList = list(
x = x ,
y = y ,
Nx = Nx ,
Ntotal = Ntotal ,
xProbe = xProbe ,
Nprobe = nrow(xProbe)
)
#-----------------------------------------------------------------------------
# THE MODEL.
modelString = "
# Standardize the data:
data {
ym <- mean(y)
ysd <- sd(y)
for ( i in 1:Ntotal ) {
zy[i] <- ( y[i] - ym ) / ysd
}
for ( j in 1:Nx ) {
xm[j] <- mean(x[,j])
xsd[j] <- sd(x[,j])
for ( i in 1:Ntotal ) {
zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
}
# standardize the probe values:
for ( i in 1:Nprobe ) {
zxProbe[i,j] <- ( xProbe[i,j] - xm[j] ) / xsd[j]
}
}
}
# Specify the model for standardized data:
model {
for ( i in 1:Ntotal ) {
zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigma^2 , nu )
}
# Priors vague on standardized scale:
zbeta0 ~ dnorm( 0 , 1/2^2 )
for ( j in 1:Nx ) {
zbeta[j] ~ dnorm( 0 , 1/2^2 )
}
zsigma ~ dunif( 1.0E-5 , 1.0E+1 )
nu ~ dexp(1/30.0)
# Transform to original scale:
beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd
beta0 <- zbeta0*ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd
sigma <- zsigma*ysd
# Predicted y values at xProbe:
for ( i in 1:Nprobe ) {
zyP[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zxProbe[i,1:Nx] ) ,
1/zsigma^2 , nu )
yP[i] <- zyP[i] * ysd + ym
}
}
" # close quote for modelString
# Write out modelString to a text file
writeLines( modelString , con="TEMPmodel.txt" )
#-----------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Let JAGS do it...
#-----------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "beta0" , "beta" , "sigma",
"zbeta0" , "zbeta" , "zsigma", "nu" , "xProbe" , "yP" )
adaptSteps = 500 # Number of steps to "tune" the samplers
burnInSteps = 1000
runJagsOut <- run.jags( method="parallel" ,
model="TEMPmodel.txt" ,
monitor=parameters ,
data=dataList ,
#inits=initsList ,
n.chains=nChains ,
adapt=adaptSteps ,
burnin=burnInSteps ,
sample=ceiling(numSavedSteps/nChains) ,
thin=thinSteps ,
summarise=FALSE ,
plots=FALSE )
codaSamples = as.mcmc.list( runJagsOut )
# resulting codaSamples object has these indices:
# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
if ( !is.null(saveName) ) {
save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") )
}
return( codaSamples )
} # end function
#===============================================================================
smryMCMC = function( codaSamples ,
saveName=NULL ) {
summaryInfo = NULL
mcmcMat = as.matrix(codaSamples,chains=TRUE)
paramName = colnames(mcmcMat)
for ( pName in paramName ) {
summaryInfo = rbind( summaryInfo , summarizePost( mcmcMat[,pName] ) )
}
rownames(summaryInfo) = paramName
summaryInfo = rbind( summaryInfo ,
"log10(nu)" = summarizePost( log10(mcmcMat[,"nu"]) ) )
if ( !is.null(saveName) ) {
write.csv( summaryInfo , file=paste(saveName,"SummaryInfo.csv",sep="") )
}
return( summaryInfo )
}
#===============================================================================
plotMCMC = function( codaSamples , data , xName="x" , yName="y" ,
showCurve=FALSE , pairsPlot=FALSE ,
saveName=NULL , saveType="jpg" ) {
# showCurve is TRUE or FALSE and indicates whether the posterior should
# be displayed as a histogram (by default) or by an approximate curve.
# pairsPlot is TRUE or FALSE and indicates whether scatterplots of pairs
# of parameters should be displayed.
#-----------------------------------------------------------------------------
y = data[,yName]
x = as.matrix(data[,xName])
mcmcMat = as.matrix(codaSamples,chains=TRUE)
chainLength = NROW( mcmcMat )
zbeta0 = mcmcMat[,"zbeta0"]
zbeta = mcmcMat[,grep("^zbeta$|^zbeta\\[",colnames(mcmcMat))]
if ( ncol(x)==1 ) { zbeta = matrix( zbeta , ncol=1 ) }
zsigma = mcmcMat[,"zsigma"]
beta0 = mcmcMat[,"beta0"]
beta = mcmcMat[,grep("^beta$|^beta\\[",colnames(mcmcMat))]
if ( ncol(x)==1 ) { beta = matrix( beta , ncol=1 ) }
sigma = mcmcMat[,"sigma"]
nu = mcmcMat[,"nu"]
log10nu = log10(nu)
#-----------------------------------------------------------------------------
# Compute R^2 for credible parameters:
YcorX = cor( y , x ) # correlation of y with each x predictor
Rsq = zbeta %*% matrix( YcorX , ncol=1 )
#-----------------------------------------------------------------------------
if ( pairsPlot ) {
# Plot the parameters pairwise, to see correlations:
openGraph()
nPtToPlot = 1000
plotIdx = floor(seq(1,chainLength,by=chainLength/nPtToPlot))
panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
usr = par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r = (cor(x, y))
txt = format(c(r, 0.123456789), digits=digits)[1]
txt = paste(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex=1.25 ) # was cex=cex.cor*r
}
pairs( cbind( beta0 , beta , sigma , log10nu )[plotIdx,] ,
labels=c( "beta[0]" ,
paste0("beta[",1:ncol(beta),"]\n",xName) ,
expression(sigma) , expression(log10(nu)) ) ,
lower.panel=panel.cor , col="skyblue" )
if ( !is.null(saveName) ) {
saveGraph( file=paste(saveName,"PostPairs",sep=""), type=saveType)
}
}
#-----------------------------------------------------------------------------
# Marginal histograms:
decideOpenGraph = function( panelCount , saveName , finished=FALSE ,
nRow=2 , nCol=3 ) {
# If finishing a set:
if ( finished==TRUE ) {
if ( !is.null(saveName) ) {
saveGraph( file=paste0(saveName,ceiling((panelCount-1)/(nRow*nCol))),
type=saveType)
}
panelCount = 1 # re-set panelCount
return(panelCount)
} else {
# If this is first panel of a graph:
if ( ( panelCount %% (nRow*nCol) ) == 1 ) {
# If previous graph was open, save previous one:
if ( panelCount>1 & !is.null(saveName) ) {
saveGraph( file=paste0(saveName,(panelCount%/%(nRow*nCol))),
type=saveType)
}
# Open new graph
openGraph(width=nCol*7.0/3,height=nRow*2.0)
layout( matrix( 1:(nRow*nCol) , nrow=nRow, byrow=TRUE ) )
par( mar=c(4,4,2.5,0.5) , mgp=c(2.5,0.7,0) )
}
# Increment and return panel count:
panelCount = panelCount+1
return(panelCount)
}
}
# Original scale:
panelCount = 1
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
histInfo = plotPost( beta0 , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(beta[0]) , main="Intercept" )
for ( bIdx in 1:ncol(beta) ) {
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
histInfo = plotPost( beta[,bIdx] , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(beta[.(bIdx)]) , main=xName[bIdx] )
}
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
histInfo = plotPost( sigma , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(sigma) , main=paste("Scale") )
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
histInfo = plotPost( log10nu , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(log10(nu)) , main=paste("Normality") )
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
histInfo = plotPost( Rsq , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(R^2) , main=paste("Prop Var Accntd") )
panelCount = decideOpenGraph( panelCount , finished=TRUE , saveName=paste0(saveName,"PostMarg") )
# Standardized scale:
panelCount = 1
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
histInfo = plotPost( zbeta0 , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(z*beta[0]) , main="Intercept" )
for ( bIdx in 1:ncol(beta) ) {
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
histInfo = plotPost( zbeta[,bIdx] , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(z*beta[.(bIdx)]) , main=xName[bIdx] )
}
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
histInfo = plotPost( zsigma , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(z*sigma) , main=paste("Scale") )
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
histInfo = plotPost( log10nu , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(log10(nu)) , main=paste("Normality") )
panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
histInfo = plotPost( Rsq , cex.lab = 1.75 , showCurve=showCurve ,
xlab=bquote(R^2) , main=paste("Prop Var Accntd") )
panelCount = decideOpenGraph( panelCount , finished=TRUE , saveName=paste0(saveName,"PostMargZ") )
#-----------------------------------------------------------------------------
}
#===============================================================================
#===== End low-level script ================