OpenTURNS
Note
Go to the end to download the full example code
Customize your Metropolis-Hastings algorithm¶
This simple example shows how you can build your own variant of the Metropolis-Hastings algorithm.
We want to sample from the distribution with support ![[-3, 3]^2](data:image/svg+xml;base64,PD94bWwgdmVyc2lvbj0nMS4wJyBlbmNvZGluZz0nVVRGLTgnPz4KPCEtLSBUaGlzIGZpbGUgd2FzIGdlbmVyYXRlZCBieSBkdmlzdmdtIDMuMS4yIC0tPgo8c3ZnIHZlcnNpb249JzEuMScgeG1sbnM9J2h0dHA6Ly93d3cudzMub3JnLzIwMDAvc3ZnJyB4bWxuczp4bGluaz0naHR0cDovL3d3dy53My5vcmcvMTk5OS94bGluaycgd2lkdGg9JzM2Ljk4NjE0MnB0JyBoZWlnaHQ9JzEyLjQ2MzUzMXB0JyB2aWV3Qm94PScwIC05LjQ3NDczOSAzNi45ODYxNDIgMTIuNDYzNTMxJz4KPGRlZnM+CjxwYXRoIGlkPSdnMS01OScgZD0nTTIuMzMxMjU4IC4wNDc4MjFDMi4zMzEyNTgtLjY0NTU3OSAyLjEwNDExLTEuMTU5NjUxIDEuNjEzOTQ4LTEuMTU5NjUxQzEuMjMxMzgyLTEuMTU5NjUxIDEuMDQwMS0uODQ4ODE3IDEuMDQwMS0uNTg1ODAzUzEuMjE5NDI3IDAgMS42MjU5MDMgMEMxLjc4MTMyIDAgMS45MTI4MjctLjA0NzgyMSAyLjAyMDQyMy0uMTU1NDE3QzIuMDQ0MzM0LS4xNzkzMjggMi4wNTYyODktLjE3OTMyOCAyLjA2ODI0NC0uMTc5MzI4QzIuMDkyMTU0LS4xNzkzMjggMi4wOTIxNTQtLjAxMTk1NSAyLjA5MjE1NCAuMDQ3ODIxQzIuMDkyMTU0IC40NDIzNDEgMi4wMjA0MjMgMS4yMTk0MjcgMS4zMjcwMjQgMS45OTY1MTNDMS4xOTU1MTcgMi4xMzk5NzUgMS4xOTU1MTcgMi4xNjM4ODUgMS4xOTU1MTcgMi4xODc3OTZDMS4xOTU1MTcgMi4yNDc1NzIgMS4yNTUyOTMgMi4zMDczNDcgMS4zMTUwNjggMi4zMDczNDdDMS40MTA3MSAyLjMwNzM0NyAyLjMzMTI1OCAxLjQyMjY2NSAyLjMzMTI1OCAuMDQ3ODIxWicvPgo8cGF0aCBpZD0nZzItNTAnIGQ9J00yLjI0NzU3Mi0xLjYyNTkwM0MyLjM3NTA5My0xLjc0NTQ1NSAyLjcwOTgzOC0yLjAwODQ2OCAyLjgzNzM2LTIuMTIwMDVDMy4zMzE1MDctMi41NzQzNDYgMy44MDE3NDMtMy4wMTI3MDIgMy44MDE3NDMtMy43Mzc5ODNDMy44MDE3NDMtNC42ODY0MjYgMy4wMDQ3MzItNS4zMDAxMjUgMi4wMDg0NjgtNS4zMDAxMjVDMS4wNTIwNTUtNS4zMDAxMjUgLjQyMjQxNi00LjU3NDg0NCAuNDIyNDE2LTMuODY1NTA0Qy40MjI0MTYtMy40NzQ5NjkgLjczMzI1LTMuNDE5MTc4IC44NDQ4MzItMy40MTkxNzhDMS4wMTIyMDQtMy40MTkxNzggMS4yNTkyNzgtMy41Mzg3MyAxLjI1OTI3OC0zLjg0MTU5NEMxLjI1OTI3OC00LjI1NjA0IC44NjA3NzItNC4yNTYwNCAuNzY1MTMxLTQuMjU2MDRDLjk5NjI2NC00LjgzNzg1OCAxLjUzMDI2Mi01LjAzNzExMSAxLjkyMDc5Ny01LjAzNzExMUMyLjY2MjAxNy01LjAzNzExMSAzLjA0NDU4My00LjQwNzQ3MiAzLjA0NDU4My0zLjczNzk4M0MzLjA0NDU4My0yLjkwOTA5MSAyLjQ2Mjc2NS0yLjMwMzM2MiAxLjUyMjI5MS0xLjMzODk3OUwuNTE4MDU3LS4zMDI4NjRDLjQyMjQxNi0uMjE1MTkzIC40MjI0MTYtLjE5OTI1MyAuNDIyNDE2IDBIMy41NzA2MUwzLjgwMTc0My0xLjQyNjY1SDMuNTU0NjdDMy41MzA3Ni0xLjI2NzI0OCAzLjQ2Njk5OS0uODY4NzQyIDMuMzcxMzU3LS43MTczMUMzLjMyMzUzNy0uNjUzNTQ5IDIuNzE3ODA4LS42NTM1NDkgMi41OTAyODYtLjY1MzU0OUgxLjE3MTYwNkwyLjI0NzU3Mi0xLjYyNTkwM1onLz4KPHBhdGggaWQ9J2cwLTAnIGQ9J003Ljg3ODQ1Ni0yLjc0OTY4OUM4LjA4MTY5NC0yLjc0OTY4OSA4LjI5Njg4Ny0yLjc0OTY4OSA4LjI5Njg4Ny0yLjk4ODc5MlM4LjA4MTY5NC0zLjIyNzg5NSA3Ljg3ODQ1Ni0zLjIyNzg5NUgxLjQxMDcxQzEuMjA3NDcyLTMuMjI3ODk1IC45OTIyNzktMy4yMjc4OTUgLjk5MjI3OS0yLjk4ODc5MlMxLjIwNzQ3Mi0yLjc0OTY4OSAxLjQxMDcxLTIuNzQ5Njg5SDcuODc4NDU2WicvPgo8cGF0aCBpZD0nZzMtNTEnIGQ9J00yLjE5OTc1MS00LjI5MTkwNUMxLjk5NjUxMy00LjI3OTk1IDEuOTQ4NjkyLTQuMjY3OTk1IDEuOTQ4NjkyLTQuMTYwMzk5QzEuOTQ4NjkyLTQuMDQwODQ3IDIuMDA4NDY4LTQuMDQwODQ3IDIuMjIzNjYxLTQuMDQwODQ3SDIuNzczNTk5QzMuNzg5Nzg4LTQuMDQwODQ3IDQuMjQ0MDg1LTMuMjAzOTg1IDQuMjQ0MDg1LTIuMDU2Mjg5QzQuMjQ0MDg1LS40OTAxNjIgMy40MzExMzMtLjA3MTczMSAyLjg0NTMzLS4wNzE3MzFDMi4yNzE0ODItLjA3MTczMSAxLjI5MTE1OC0uMzQ2NyAuOTQ0NDU4LTEuMTM1NzQxQzEuMzI3MDI0LTEuMDc1OTY1IDEuNjczNzI0LTEuMjkxMTU4IDEuNjczNzI0LTEuNzIxNTQ0QzEuNjczNzI0LTIuMDY4MjQ0IDEuNDIyNjY1LTIuMzA3MzQ3IDEuMDg3OTItMi4zMDczNDdDLjgwMDk5Ni0yLjMwNzM0NyAuNDkwMTYyLTIuMTM5OTc1IC40OTAxNjItMS42ODU2NzlDLjQ5MDE2Mi0uNjIxNjY5IDEuNTU0MTcyIC4yNTEwNTkgMi44ODExOTYgLjI1MTA1OUM0LjMwMzg2MSAuMjUxMDU5IDUuMzU1OTE1LS44MzY4NjIgNS4zNTU5MTUtMi4wNDQzMzRDNS4zNTU5MTUtMy4xNDQyMDkgNC40NzEyMzMtNC4wMDQ5ODEgMy4zMjM1MzctNC4yMDgyMTlDNC4zNjM2MzYtNC41MDcwOTggNS4wMzMxMjYtNS4zNzk4MjYgNS4wMzMxMjYtNi4zMTIzMjlDNS4wMzMxMjYtNy4yNTY3ODcgNC4wNTI4MDItNy45NTAxODcgMi44OTMxNTEtNy45NTAxODdDMS42OTc2MzQtNy45NTAxODcgLjgxMjk1MS03LjIyMDkyMiAuODEyOTUxLTYuMzQ4MTk0Qy44MTI5NTEtNS44Njk5ODggMS4xODM1NjItNS43NzQzNDYgMS4zNjI4ODktNS43NzQzNDZDMS42MTM5NDgtNS43NzQzNDYgMS45MDA4NzItNS45NTM2NzQgMS45MDA4NzItNi4zMTIzMjlDMS45MDA4NzItNi42OTQ4OTQgMS42MTM5NDgtNi44NjIyNjcgMS4zNTA5MzQtNi44NjIyNjdDMS4yNzkyMDMtNi44NjIyNjcgMS4yNTUyOTMtNi44NjIyNjcgMS4yMTk0MjctNi44NTAzMTFDMS42NzM3MjQtNy42NjMyNjMgMi43OTc1MDktNy42NjMyNjMgMi44NTcyODUtNy42NjMyNjNDMy4yNTE4MDYtNy42NjMyNjMgNC4wMjg4OTItNy40ODM5MzUgNC4wMjg4OTItNi4zMTIzMjlDNC4wMjg4OTItNi4wODUxODEgMy45OTMwMjYtNS40MTU2OTEgMy42NDYzMjYtNC45MDE2MTlDMy4yODc2NzEtNC4zNzU1OTIgMi44ODExOTYtNC4zMzk3MjYgMi41NTg0MDYtNC4zMjc3NzFMMi4xOTk3NTEtNC4yOTE5MDVaJy8+CjxwYXRoIGlkPSdnMy05MScgZD0nTTIuOTg4NzkyIDIuOTg4NzkyVjIuNTQ2NDUxSDEuODI5MTQxVi04LjUyNDAzNUgyLjk4ODc5MlYtOC45NjYzNzZIMS4zODY4VjIuOTg4NzkySDIuOTg4NzkyWicvPgo8cGF0aCBpZD0nZzMtOTMnIGQ9J00xLjg1MzA1MS04Ljk2NjM3NkguMjUxMDU5Vi04LjUyNDAzNUgxLjQxMDcxVjIuNTQ2NDUxSC4yNTEwNTlWMi45ODg3OTJIMS44NTMwNTFWLTguOTY2Mzc2WicvPgo8L2RlZnM+CjxnIGlkPSdwYWdlMSc+Cjx1c2UgeD0nMCcgeT0nMCcgeGxpbms6aHJlZj0nI2czLTkxJy8+Cjx1c2UgeD0nMy4yNTE2NjEnIHk9JzAnIHhsaW5rOmhyZWY9JyNnMC0wJy8+Cjx1c2UgeD0nMTIuNTUwMTU4JyB5PScwJyB4bGluazpocmVmPScjZzMtNTEnLz4KPHVzZSB4PScxOC40MDMxNDknIHk9JzAnIHhsaW5rOmhyZWY9JyNnMS01OScvPgo8dXNlIHg9JzIzLjY0NzMwNycgeT0nMCcgeGxpbms6aHJlZj0nI2czLTUxJy8+Cjx1c2UgeD0nMjkuNTAwMjk4JyB5PScwJyB4bGluazpocmVmPScjZzMtOTMnLz4KPHVzZSB4PSczMi43NTE5NTknIHk9Jy00LjMzODQzNycgeGxpbms6aHJlZj0nI2cyLTUwJy8+CjwvZz4KPC9zdmc+CjwhLS0gREVQVEg9NCAtLT4=)
whose PDF 
is proportional to the Ackley function to the tenth power:
![\forall \vect{x} \in [-3, 3]^2 \quad f(\vect{x}) \propto a(\vect{x})^{10},](data:image/svg+xml;base64,PD94bWwgdmVyc2lvbj0nMS4wJyBlbmNvZGluZz0nVVRGLTgnPz4KPCEtLSBUaGlzIGZpbGUgd2FzIGdlbmVyYXRlZCBieSBkdmlzdmdtIDMuMS4yIC0tPgo8c3ZnIHZlcnNpb249JzEuMScgeG1sbnM9J2h0dHA6Ly93d3cudzMub3JnLzIwMDAvc3ZnJyB4bWxuczp4bGluaz0naHR0cDovL3d3dy53My5vcmcvMTk5OS94bGluaycgd2lkdGg9JzE1My42MzkxN3B0JyBoZWlnaHQ9JzEzLjA2MTI4cHQnIHZpZXdCb3g9JzExNy40NTE5MDUgLTE0LjMwMTU2MyAxNTMuNjM5MTcgMTMuMDYxMjgnPgo8ZGVmcz4KPHBhdGggaWQ9J2czLTQ4JyBkPSdNMy44OTczODUtMi41NDI0NjZDMy44OTczODUtMy4zOTUyNjggMy44MDk3MTQtMy45MTMzMjUgMy41NDY3LTQuNDIzNDEyQzMuMTk2MDE1LTUuMTI0NzgyIDIuNTUwNDM2LTUuMzAwMTI1IDIuMTEyMDgtNS4zMDAxMjVDMS4xMDc4NDYtNS4zMDAxMjUgLjc0MTIyLTQuNTUwOTM0IC42Mjk2MzktNC4zMjc3NzFDLjM0MjcxNS0zLjc0NTk1MyAuMzI2Nzc1LTIuOTU2OTEyIC4zMjY3NzUtMi41NDI0NjZDLjMyNjc3NS0yLjAxNjQzOCAuMzUwNjg1LTEuMjExNDU3IC43MzMyNS0uNTczODQ4QzEuMDk5ODc1IC4wMTU5NCAxLjY4OTY2NCAuMTY3MzcyIDIuMTEyMDggLjE2NzM3MkMyLjQ5NDY0NSAuMTY3MzcyIDMuMTgwMDc1IC4wNDc4MjEgMy41Nzg1OC0uNzQxMjJDMy44NzM0NzQtMS4zMTUwNjggMy44OTczODUtMi4wMjQ0MDggMy44OTczODUtMi41NDI0NjZaTTIuMTEyMDgtLjA1NTc5MUMxLjg0MTA5Ni0uMDU1NzkxIDEuMjkxMTU4LS4xODMzMTMgMS4xMjM3ODYtMS4wMjAxNzRDMS4wMzYxMTUtMS40NzQ0NzEgMS4wMzYxMTUtMi4yMjM2NjEgMS4wMzYxMTUtMi42MzgxMDdDMS4wMzYxMTUtMy4xODgwNDUgMS4wMzYxMTUtMy43NDU5NTMgMS4xMjM3ODYtNC4xODQzMDlDMS4yOTExNTgtNC45OTcyNiAxLjkxMjgyNy01LjA3Njk2MSAyLjExMjA4LTUuMDc2OTYxQzIuMzgzMDY0LTUuMDc2OTYxIDIuOTMzMDAxLTQuOTQxNDY5IDMuMDkyNDAzLTQuMjE2MTg5QzMuMTg4MDQ1LTMuNzc3ODMzIDMuMTg4MDQ1LTMuMTgwMDc1IDMuMTg4MDQ1LTIuNjM4MTA3QzMuMTg4MDQ1LTIuMTY3ODcgMy4xODgwNDUtMS40NTA1NiAzLjA5MjQwMy0xLjAwNDIzNEMyLjkyNTAzMS0uMTY3MzcyIDIuMzc1MDkzLS4wNTU3OTEgMi4xMTIwOC0uMDU1NzkxWicvPgo8cGF0aCBpZD0nZzMtNDknIGQ9J00yLjUwMjYxNS01LjA3Njk2MUMyLjUwMjYxNS01LjI5MjE1NCAyLjQ4NjY3NS01LjMwMDEyNSAyLjI3MTQ4Mi01LjMwMDEyNUMxLjk0NDcwNy00Ljk4MTMyIDEuNTIyMjkxLTQuNzkwMDM3IC43NjUxMzEtNC43OTAwMzdWLTQuNTI3MDI0Qy45ODAzMjQtNC41MjcwMjQgMS40MTA3MS00LjUyNzAyNCAxLjg3Mjk3Ni00Ljc0MjIxN1YtLjY1MzU0OUMxLjg3Mjk3Ni0uMzU4NjU1IDEuODQ5MDY2LS4yNjMwMTQgMS4wOTE5MDUtLjI2MzAxNEguODEyOTUxVjBDMS4xMzk3MjYtLjAyMzkxIDEuODI1MTU2LS4wMjM5MSAyLjE4MzgxMS0uMDIzOTFTMy4yMzU4NjYtLjAyMzkxIDMuNTYyNjQgMFYtLjI2MzAxNEgzLjI4MzY4NkMyLjUyNjUyNi0uMjYzMDE0IDIuNTAyNjE1LS4zNTg2NTUgMi41MDI2MTUtLjY1MzU0OVYtNS4wNzY5NjFaJy8+CjxwYXRoIGlkPSdnMy01MCcgZD0nTTIuMjQ3NTcyLTEuNjI1OTAzQzIuMzc1MDkzLTEuNzQ1NDU1IDIuNzA5ODM4LTIuMDA4NDY4IDIuODM3MzYtMi4xMjAwNUMzLjMzMTUwNy0yLjU3NDM0NiAzLjgwMTc0My0zLjAxMjcwMiAzLjgwMTc0My0zLjczNzk4M0MzLjgwMTc0My00LjY4NjQyNiAzLjAwNDczMi01LjMwMDEyNSAyLjAwODQ2OC01LjMwMDEyNUMxLjA1MjA1NS01LjMwMDEyNSAuNDIyNDE2LTQuNTc0ODQ0IC40MjI0MTYtMy44NjU1MDRDLjQyMjQxNi0zLjQ3NDk2OSAuNzMzMjUtMy40MTkxNzggLjg0NDgzMi0zLjQxOTE3OEMxLjAxMjIwNC0zLjQxOTE3OCAxLjI1OTI3OC0zLjUzODczIDEuMjU5Mjc4LTMuODQxNTk0QzEuMjU5Mjc4LTQuMjU2MDQgLjg2MDc3Mi00LjI1NjA0IC43NjUxMzEtNC4yNTYwNEMuOTk2MjY0LTQuODM3ODU4IDEuNTMwMjYyLTUuMDM3MTExIDEuOTIwNzk3LTUuMDM3MTExQzIuNjYyMDE3LTUuMDM3MTExIDMuMDQ0NTgzLTQuNDA3NDcyIDMuMDQ0NTgzLTMuNzM3OTgzQzMuMDQ0NTgzLTIuOTA5MDkxIDIuNDYyNzY1LTIuMzAzMzYyIDEuNTIyMjkxLTEuMzM4OTc5TC41MTgwNTctLjMwMjg2NEMuNDIyNDE2LS4yMTUxOTMgLjQyMjQxNi0uMTk5MjUzIC40MjI0MTYgMEgzLjU3MDYxTDMuODAxNzQzLTEuNDI2NjVIMy41NTQ2N0MzLjUzMDc2LTEuMjY3MjQ4IDMuNDY2OTk5LS44Njg3NDIgMy4zNzEzNTctLjcxNzMxQzMuMzIzNTM3LS42NTM1NDkgMi43MTc4MDgtLjY1MzU0OSAyLjU5MDI4Ni0uNjUzNTQ5SDEuMTcxNjA2TDIuMjQ3NTcyLTEuNjI1OTAzWicvPgo8cGF0aCBpZD0nZzItNTknIGQ9J00yLjMzMTI1OCAuMDQ3ODIxQzIuMzMxMjU4LS42NDU1NzkgMi4xMDQxMS0xLjE1OTY1MSAxLjYxMzk0OC0xLjE1OTY1MUMxLjIzMTM4Mi0xLjE1OTY1MSAxLjA0MDEtLjg0ODgxNyAxLjA0MDEtLjU4NTgwM1MxLjIxOTQyNyAwIDEuNjI1OTAzIDBDMS43ODEzMiAwIDEuOTEyODI3LS4wNDc4MjEgMi4wMjA0MjMtLjE1NTQxN0MyLjA0NDMzNC0uMTc5MzI4IDIuMDU2Mjg5LS4xNzkzMjggMi4wNjgyNDQtLjE3OTMyOEMyLjA5MjE1NC0uMTc5MzI4IDIuMDkyMTU0LS4wMTE5NTUgMi4wOTIxNTQgLjA0NzgyMUMyLjA5MjE1NCAuNDQyMzQxIDIuMDIwNDIzIDEuMjE5NDI3IDEuMzI3MDI0IDEuOTk2NTEzQzEuMTk1NTE3IDIuMTM5OTc1IDEuMTk1NTE3IDIuMTYzODg1IDEuMTk1NTE3IDIuMTg3Nzk2QzEuMTk1NTE3IDIuMjQ3NTcyIDEuMjU1MjkzIDIuMzA3MzQ3IDEuMzE1MDY4IDIuMzA3MzQ3QzEuNDEwNzEgMi4zMDczNDcgMi4zMzEyNTggMS40MjI2NjUgMi4zMzEyNTggLjA0NzgyMVonLz4KPHBhdGggaWQ9J2cyLTk3JyBkPSdNMy41OTg1MDYtMS40MjI2NjVDMy41Mzg3My0xLjIxOTQyNyAzLjUzODczLTEuMTk1NTE3IDMuMzcxMzU3LS45NjgzNjlDMy4xMDgzNDQtLjYzMzYyNCAyLjU4MjMxNi0uMTE5NTUyIDIuMDIwNDIzLS4xMTk1NTJDMS41MzAyNjItLjExOTU1MiAxLjI1NTI5My0uNTYxODkzIDEuMjU1MjkzLTEuMjY3MjQ4QzEuMjU1MjkzLTEuOTI0NzgyIDEuNjI1OTAzLTMuMjYzNzYxIDEuODUzMDUxLTMuNzY1ODc4QzIuMjU5NTI3LTQuNjAyNzQgMi44MjE0Mi01LjAzMzEyNiAzLjI4NzY3MS01LjAzMzEyNkM0LjA3NjcxMi01LjAzMzEyNiA0LjIzMjEzLTQuMDUyODAyIDQuMjMyMTMtMy45NTcxNjFDNC4yMzIxMy0zLjk0NTIwNSA0LjE5NjI2NC0zLjc4OTc4OCA0LjE4NDMwOS0zLjc2NTg3OEwzLjU5ODUwNi0xLjQyMjY2NVpNNC4zNjM2MzYtNC40ODMxODhDNC4yMzIxMy00Ljc5NDAyMiAzLjkwOTM0LTUuMjcyMjI5IDMuMjg3NjcxLTUuMjcyMjI5QzEuOTM2NzM3LTUuMjcyMjI5IC40NzgyMDctMy41MjY3NzUgLjQ3ODIwNy0xLjc1NzQxQy40NzgyMDctLjU3Mzg0OCAxLjE3MTYwNiAuMTE5NTUyIDEuOTg0NTU4IC4xMTk1NTJDMi42NDIwOTIgLjExOTU1MiAzLjIwMzk4NS0uMzk0NTIxIDMuNTM4NzMtLjc4OTA0MUMzLjY1ODI4MS0uMDgzNjg2IDQuMjIwMTc0IC4xMTk1NTIgNC41Nzg4MjkgLjExOTU1MlM1LjIyNDQwOC0uMDk1NjQxIDUuNDM5NjAxLS41MjYwMjdDNS42MzA4ODQtLjkzMjUwMyA1Ljc5ODI1Ny0xLjY2MTc2OCA1Ljc5ODI1Ny0xLjcwOTU4OUM1Ljc5ODI1Ny0xLjc2OTM2NSA1Ljc1MDQzNi0xLjgxNzE4NiA1LjY3ODcwNS0xLjgxNzE4NkM1LjU3MTEwOC0xLjgxNzE4NiA1LjU1OTE1My0xLjc1NzQxIDUuNTExMzMzLTEuNTc4MDgyQzUuMzMyMDA1LS44NzI3MjcgNS4xMDQ4NTctLjExOTU1MiA0LjYxNDY5NS0uMTE5NTUyQzQuMjY3OTk1LS4xMTk1NTIgNC4yNDQwODUtLjQzMDM4NiA0LjI0NDA4NS0uNjY5NDg5QzQuMjQ0MDg1LS45NDQ0NTggNC4yNzk5NS0xLjA3NTk2NSA0LjM4NzU0Ny0xLjU0MjIxN0M0LjQ3MTIzMy0xLjg0MTA5NiA0LjUzMTAwOS0yLjEwNDExIDQuNjI2NjUtMi40NTA4MDlDNS4wNjg5OTEtNC4yNDQwODUgNS4xNzY1ODgtNC42NzQ0NzEgNS4xNzY1ODgtNC43NDYyMDJDNS4xNzY1ODgtNC45MTM1NzQgNS4wNDUwODEtNS4wNDUwODEgNC44NjU3NTMtNS4wNDUwODFDNC40ODMxODgtNS4wNDUwODEgNC4zODc1NDctNC42MjY2NSA0LjM2MzYzNi00LjQ4MzE4OFonLz4KPHBhdGggaWQ9J2cyLTEwMicgZD0nTTUuMzMyMDA1LTQuODA1OTc4QzUuNTcxMTA4LTQuODA1OTc4IDUuNjY2NzUtNC44MDU5NzggNS42NjY3NS01LjAzMzEyNkM1LjY2Njc1LTUuMTUyNjc3IDUuNTcxMTA4LTUuMTUyNjc3IDUuMzU1OTE1LTUuMTUyNjc3SDQuMzg3NTQ3QzQuNjE0Njk1LTYuMzg0MDYgNC43ODIwNjctNy4yMzI4NzcgNC44Nzc3MDktNy42MTU0NDJDNC45NDk0NC03LjkwMjM2NiA1LjIwMDQ5OC04LjE3NzMzNSA1LjUxMTMzMy04LjE3NzMzNUM1Ljc2MjM5MS04LjE3NzMzNSA2LjAxMzQ1LTguMDY5NzM4IDYuMTMzMDAxLTcuOTYyMTQyQzUuNjY2NzUtNy45MTQzMjEgNS41MjMyODgtNy41Njc2MjEgNS41MjMyODgtNy4zNjQzODRDNS41MjMyODgtNy4xMjUyOCA1LjcwMjYxNS02Ljk4MTgxOCA1LjkyOTc2My02Ljk4MTgxOEM2LjE2ODg2Ny02Ljk4MTgxOCA2LjUyNzUyMi03LjE4NTA1NiA2LjUyNzUyMi03LjYzOTM1MkM2LjUyNzUyMi04LjE0MTQ2OSA2LjAyNTQwNS04LjQxNjQzOCA1LjQ5OTM3Ny04LjQxNjQzOEM0Ljk4NTMwNS04LjQxNjQzOCA0LjQ4MzE4OC04LjAzMzg3MyA0LjI0NDA4NS03LjU2NzYyMUM0LjAyODg5Mi03LjE0OTE5MSAzLjkwOTM0LTYuNzE4ODA0IDMuNjM0MzcxLTUuMTUyNjc3SDIuODMzMzc1QzIuNjA2MjI3LTUuMTUyNjc3IDIuNDg2Njc1LTUuMTUyNjc3IDIuNDg2Njc1LTQuOTM3NDg0QzIuNDg2Njc1LTQuODA1OTc4IDIuNTU4NDA2LTQuODA1OTc4IDIuNzk3NTA5LTQuODA1OTc4SDMuNTYyNjRDMy4zNDc0NDctMy42OTQxNDcgMi44NTcyODUtLjk5MjI3OSAyLjU4MjMxNiAuMjg2OTI0QzIuMzc5MDc4IDEuMzI3MDI0IDIuMTk5NzUxIDIuMTk5NzUxIDEuNjAxOTkzIDIuMTk5NzUxQzEuNTY2MTI3IDIuMTk5NzUxIDEuMjE5NDI3IDIuMTk5NzUxIDEuMDA0MjM0IDEuOTcyNjAzQzEuNjEzOTQ4IDEuOTI0NzgyIDEuNjEzOTQ4IDEuMzk4NzU1IDEuNjEzOTQ4IDEuMzg2OEMxLjYxMzk0OCAxLjE0NzY5NiAxLjQzNDYyIDEuMDA0MjM0IDEuMjA3NDcyIDEuMDA0MjM0Qy45NjgzNjkgMS4wMDQyMzQgLjYwOTcxNCAxLjIwNzQ3MiAuNjA5NzE0IDEuNjYxNzY4Qy42MDk3MTQgMi4xNzU4NDEgMS4xMzU3NDEgMi40Mzg4NTQgMS42MDE5OTMgMi40Mzg4NTRDMi44MjE0MiAyLjQzODg1NCAzLjMyMzUzNyAuMjUxMDU5IDMuNDU1MDQ0LS4zNDY3QzMuNjcwMjM3LTEuMjY3MjQ4IDQuMjU2MDQtNC40NDczMjMgNC4zMTU4MTYtNC44MDU5NzhINS4zMzIwMDVaJy8+CjxwYXRoIGlkPSdnNC00MCcgZD0nTTMuODg1NDMgMi45MDUxMDZDMy44ODU0MyAyLjg2OTI0IDMuODg1NDMgMi44NDUzMyAzLjY4MjE5MiAyLjY0MjA5MkMyLjQ4NjY3NSAxLjQzNDYyIDEuODE3MTg2LS41Mzc5ODMgMS44MTcxODYtMi45NzY4MzdDMS44MTcxODYtNS4yOTYxMzkgMi4zNzkwNzgtNy4yOTI2NTMgMy43NjU4NzgtOC43MDMzNjJDMy44ODU0My04LjgxMDk1OSAzLjg4NTQzLTguODM0ODY5IDMuODg1NDMtOC44NzA3MzVDMy44ODU0My04Ljk0MjQ2NiAzLjgyNTY1NC04Ljk2NjM3NiAzLjc3NzgzMy04Ljk2NjM3NkMzLjYyMjQxNi04Ljk2NjM3NiAyLjY0MjA5Mi04LjEwNTYwNCAyLjA1NjI4OS02LjkzMzk5OEMxLjQ0NjU3NS01LjcyNjUyNiAxLjE3MTYwNi00LjQ0NzMyMyAxLjE3MTYwNi0yLjk3NjgzN0MxLjE3MTYwNi0xLjkxMjgyNyAxLjMzODk3OS0uNDkwMTYyIDEuOTYwNjQ4IC43ODkwNDFDMi42NjYwMDIgMi4yMjM2NjEgMy42NDYzMjYgMy4wMDA3NDcgMy43Nzc4MzMgMy4wMDA3NDdDMy44MjU2NTQgMy4wMDA3NDcgMy44ODU0MyAyLjk3NjgzNyAzLjg4NTQzIDIuOTA1MTA2WicvPgo8cGF0aCBpZD0nZzQtNDEnIGQ9J00zLjM3MTM1Ny0yLjk3NjgzN0MzLjM3MTM1Ny0zLjg4NTQzIDMuMjUxODA2LTUuMzY3ODcgMi41ODIzMTYtNi43NTQ2N0MxLjg3Njk2MS04LjE4OTI5IC44OTY2MzgtOC45NjYzNzYgLjc2NTEzMS04Ljk2NjM3NkMuNzE3MzEtOC45NjYzNzYgLjY1NzUzNC04Ljk0MjQ2NiAuNjU3NTM0LTguODcwNzM1Qy42NTc1MzQtOC44MzQ4NjkgLjY1NzUzNC04LjgxMDk1OSAuODYwNzcyLTguNjA3NzIxQzIuMDU2Mjg5LTcuNDAwMjQ5IDIuNzI1Nzc4LTUuNDI3NjQ2IDIuNzI1Nzc4LTIuOTg4NzkyQzIuNzI1Nzc4LS42Njk0ODkgMi4xNjM4ODUgMS4zMjcwMjQgLjc3NzA4NiAyLjczNzczM0MuNjU3NTM0IDIuODQ1MzMgLjY1NzUzNCAyLjg2OTI0IC42NTc1MzQgMi45MDUxMDZDLjY1NzUzNCAyLjk3NjgzNyAuNzE3MzEgMy4wMDA3NDcgLjc2NTEzMSAzLjAwMDc0N0MuOTIwNTQ4IDMuMDAwNzQ3IDEuOTAwODcyIDIuMTM5OTc1IDIuNDg2Njc1IC45NjgzNjlDMy4wOTYzODktLjI1MTA1OSAzLjM3MTM1Ny0xLjU0MjIxNyAzLjM3MTM1Ny0yLjk3NjgzN1onLz4KPHBhdGggaWQ9J2c0LTUxJyBkPSdNMi4xOTk3NTEtNC4yOTE5MDVDMS45OTY1MTMtNC4yNzk5NSAxLjk0ODY5Mi00LjI2Nzk5NSAxLjk0ODY5Mi00LjE2MDM5OUMxLjk0ODY5Mi00LjA0MDg0NyAyLjAwODQ2OC00LjA0MDg0NyAyLjIyMzY2MS00LjA0MDg0N0gyLjc3MzU5OUMzLjc4OTc4OC00LjA0MDg0NyA0LjI0NDA4NS0zLjIwMzk4NSA0LjI0NDA4NS0yLjA1NjI4OUM0LjI0NDA4NS0uNDkwMTYyIDMuNDMxMTMzLS4wNzE3MzEgMi44NDUzMy0uMDcxNzMxQzIuMjcxNDgyLS4wNzE3MzEgMS4yOTExNTgtLjM0NjcgLjk0NDQ1OC0xLjEzNTc0MUMxLjMyNzAyNC0xLjA3NTk2NSAxLjY3MzcyNC0xLjI5MTE1OCAxLjY3MzcyNC0xLjcyMTU0NEMxLjY3MzcyNC0yLjA2ODI0NCAxLjQyMjY2NS0yLjMwNzM0NyAxLjA4NzkyLTIuMzA3MzQ3Qy44MDA5OTYtMi4zMDczNDcgLjQ5MDE2Mi0yLjEzOTk3NSAuNDkwMTYyLTEuNjg1Njc5Qy40OTAxNjItLjYyMTY2OSAxLjU1NDE3MiAuMjUxMDU5IDIuODgxMTk2IC4yNTEwNTlDNC4zMDM4NjEgLjI1MTA1OSA1LjM1NTkxNS0uODM2ODYyIDUuMzU1OTE1LTIuMDQ0MzM0QzUuMzU1OTE1LTMuMTQ0MjA5IDQuNDcxMjMzLTQuMDA0OTgxIDMuMzIzNTM3LTQuMjA4MjE5QzQuMzYzNjM2LTQuNTA3MDk4IDUuMDMzMTI2LTUuMzc5ODI2IDUuMDMzMTI2LTYuMzEyMzI5QzUuMDMzMTI2LTcuMjU2Nzg3IDQuMDUyODAyLTcuOTUwMTg3IDIuODkzMTUxLTcuOTUwMTg3QzEuNjk3NjM0LTcuOTUwMTg3IC44MTI5NTEtNy4yMjA5MjIgLjgxMjk1MS02LjM0ODE5NEMuODEyOTUxLTUuODY5OTg4IDEuMTgzNTYyLTUuNzc0MzQ2IDEuMzYyODg5LTUuNzc0MzQ2QzEuNjEzOTQ4LTUuNzc0MzQ2IDEuOTAwODcyLTUuOTUzNjc0IDEuOTAwODcyLTYuMzEyMzI5QzEuOTAwODcyLTYuNjk0ODk0IDEuNjEzOTQ4LTYuODYyMjY3IDEuMzUwOTM0LTYuODYyMjY3QzEuMjc5MjAzLTYuODYyMjY3IDEuMjU1MjkzLTYuODYyMjY3IDEuMjE5NDI3LTYuODUwMzExQzEuNjczNzI0LTcuNjYzMjYzIDIuNzk3NTA5LTcuNjYzMjYzIDIuODU3Mjg1LTcuNjYzMjYzQzMuMjUxODA2LTcuNjYzMjYzIDQuMDI4ODkyLTcuNDgzOTM1IDQuMDI4ODkyLTYuMzEyMzI5QzQuMDI4ODkyLTYuMDg1MTgxIDMuOTkzMDI2LTUuNDE1NjkxIDMuNjQ2MzI2LTQuOTAxNjE5QzMuMjg3NjcxLTQuMzc1NTkyIDIuODgxMTk2LTQuMzM5NzI2IDIuNTU4NDA2LTQuMzI3NzcxTDIuMTk5NzUxLTQuMjkxOTA1WicvPgo8cGF0aCBpZD0nZzQtOTEnIGQ9J00yLjk4ODc5MiAyLjk4ODc5MlYyLjU0NjQ1MUgxLjgyOTE0MVYtOC41MjQwMzVIMi45ODg3OTJWLTguOTY2Mzc2SDEuMzg2OFYyLjk4ODc5MkgyLjk4ODc5MlonLz4KPHBhdGggaWQ9J2c0LTkzJyBkPSdNMS44NTMwNTEtOC45NjYzNzZILjI1MTA1OVYtOC41MjQwMzVIMS40MTA3MVYyLjU0NjQ1MUguMjUxMDU5VjIuOTg4NzkySDEuODUzMDUxVi04Ljk2NjM3NlonLz4KPHBhdGggaWQ9J2cwLTEyMCcgZD0nTTYuNDA3OTctNC43OTQwMjJDNS45Nzc1ODQtNC42NzQ0NzEgNS43NjIzOTEtNC4yNjc5OTUgNS43NjIzOTEtMy45NjkxMTZDNS43NjIzOTEtMy43MDYxMDIgNS45NjU2MjktMy40MTkxNzggNi4zNjAxNDktMy40MTkxNzhDNi43Nzg1OC0zLjQxOTE3OCA3LjIyMDkyMi0zLjc2NTg3OCA3LjIyMDkyMi00LjM1MTY4MUM3LjIyMDkyMi00Ljk4NTMwNSA2LjU4NzI5OC01LjQwMzczNiA1Ljg1ODAzMi01LjQwMzczNkM1LjE3NjU4OC01LjQwMzczNiA0LjczNDI0Ny00Ljg4OTY2NCA0LjU3ODgyOS00LjY3NDQ3MUM0LjI3OTk1LTUuMTc2NTg4IDMuNjEwNDYxLTUuNDAzNzM2IDIuOTI5MDE2LTUuNDAzNzM2QzEuNDIyNjY1LTUuNDAzNzM2IC42MDk3MTQtMy45MzMyNSAuNjA5NzE0LTMuNTM4NzNDLjYwOTcxNC0zLjM3MTM1NyAuNzg5MDQxLTMuMzcxMzU3IC44OTY2MzgtMy4zNzEzNTdDMS4wNDAxLTMuMzcxMzU3IDEuMTIzNzg2LTMuMzcxMzU3IDEuMTcxNjA2LTMuNTI2Nzc1QzEuNTE4MzA2LTQuNjE0Njk1IDIuMzc5MDc4LTQuOTczMzUgMi44NjkyNC00Ljk3MzM1QzMuMzIzNTM3LTQuOTczMzUgMy41Mzg3My00Ljc1ODE1NyAzLjUzODczLTQuMzc1NTkyQzMuNTM4NzMtNC4xNDg0NDMgMy4zNzEzNTctMy40OTA5MDkgMy4yNjM3NjEtMy4wNjA1MjNMMi44NTcyODUtMS40MjI2NjVDMi42Nzc5NTgtLjY5MzQgMi4yNDc1NzItLjMzNDc0NSAxLjg0MTA5Ni0uMzM0NzQ1QzEuNzgxMzItLjMzNDc0NSAxLjUwNjM1MS0uMzM0NzQ1IDEuMjY3MjQ4LS41MTQwNzJDMS42OTc2MzQtLjYzMzYyNCAxLjkxMjgyNy0xLjA0MDEgMS45MTI4MjctMS4zMzg5NzlDMS45MTI4MjctMS42MDE5OTMgMS43MDk1ODktMS44ODg5MTcgMS4zMTUwNjgtMS44ODg5MTdDLjg5NjYzOC0xLjg4ODkxNyAuNDU0Mjk2LTEuNTQyMjE3IC40NTQyOTYtLjk1NjQxM0MuNDU0Mjk2LS4zMjI3OSAxLjA4NzkyIC4wOTU2NDEgMS44MTcxODYgLjA5NTY0MUMyLjQ5ODYzIC4wOTU2NDEgMi45NDA5NzEtLjQxODQzMSAzLjA5NjM4OS0uNjMzNjI0QzMuMzk1MjY4LS4xMzE1MDcgNC4wNjQ3NTcgLjA5NTY0MSA0Ljc0NjIwMiAuMDk1NjQxQzYuMjUyNTUzIC4wOTU2NDEgNy4wNjU1MDQtMS4zNzQ4NDQgNy4wNjU1MDQtMS43NjkzNjVDNy4wNjU1MDQtMS45MzY3MzcgNi44ODYxNzctMS45MzY3MzcgNi43Nzg1OC0xLjkzNjczN0M2LjYzNTExOC0xLjkzNjczNyA2LjU1MTQzMi0xLjkzNjczNyA2LjUwMzYxMS0xLjc4MTMyQzYuMTU2OTEyLS42OTM0IDUuMjk2MTM5LS4zMzQ3NDUgNC44MDU5NzgtLjMzNDc0NUM0LjM1MTY4MS0uMzM0NzQ1IDQuMTM2NDg4LS41NDk5MzggNC4xMzY0ODgtLjkzMjUwM0M0LjEzNjQ4OC0xLjE4MzU2MiA0LjI5MTkwNS0xLjgxNzE4NiA0LjM5OTUwMi0yLjI1OTUyN0M0LjQ4MzE4OC0yLjU3MDM2MSA0Ljc1ODE1Ny0zLjY5NDE0NyA0LjgxNzkzMy0zLjg4NTQzQzQuOTk3MjYtNC42MDI3NCA1LjQxNTY5MS00Ljk3MzM1IDUuODM0MTIyLTQuOTczMzVDNS44OTM4OTgtNC45NzMzNSA2LjE2ODg2Ny00Ljk3MzM1IDYuNDA3OTctNC43OTQwMjJaJy8+CjxwYXRoIGlkPSdnMS0wJyBkPSdNNy44Nzg0NTYtMi43NDk2ODlDOC4wODE2OTQtMi43NDk2ODkgOC4yOTY4ODctMi43NDk2ODkgOC4yOTY4ODctMi45ODg3OTJTOC4wODE2OTQtMy4yMjc4OTUgNy44Nzg0NTYtMy4yMjc4OTVIMS40MTA3MUMxLjIwNzQ3Mi0zLjIyNzg5NSAuOTkyMjc5LTMuMjI3ODk1IC45OTIyNzktMi45ODg3OTJTMS4yMDc0NzItMi43NDk2ODkgMS40MTA3MS0yLjc0OTY4OUg3Ljg3ODQ1NlonLz4KPHBhdGggaWQ9J2cxLTQ3JyBkPSdNOC42MzE2MzEtLjM4MjU2NUM4LjU4MzgxMS0uMzgyNTY1IDguMzkyNTI4LS4zNTg2NTUgOC4zNTY2NjMtLjM1ODY1NUM3LjQxMjIwNC0uMzU4NjU1IDYuNzc4NTgtMS4yOTExNTggNi4zNzIxMDUtMS45MTI4MjdDNi4yNTI1NTMtMi4xMTYwNjUgNS45Mjk3NjMtMi42MDYyMjcgNS43OTgyNTctMi44MDk0NjVDNi4wODUxODEtMy40NTUwNDQgNi44ODYxNzctNC45Mzc0ODQgOC4zMDg4NDItNC45Mzc0ODRDOC4zOTI1MjgtNC45Mzc0ODQgOC41MDAxMjUtNC45Mzc0ODQgOC42MzE2MzEtNC45MDE2MTlDOC42MzE2MzEtNS4yMTI0NTMgOC42MTk2NzYtNS4yMjQ0MDggOC41OTU3NjYtNS4yNDgzMTlDOC41MTIwOC01LjI3MjIyOSA4LjMzMjc1Mi01LjI4NDE4NCA4LjIyNTE1Ni01LjI4NDE4NEM2LjcwNjg0OS01LjI4NDE4NCA1LjgyMjE2Ny0zLjgyNTY1NCA1LjUzNTI0My0zLjIyNzg5NUM1LjA1NzAzNi0zLjk1NzE2MSA0Ljg4OTY2NC00LjIyMDE3NCA0LjQ0NzMyMy00LjU5MDc4NUMzLjcxODA1Ny01LjIyNDQwOCAzLjA2MDUyMy01LjI4NDE4NCAyLjc2MTY0NC01LjI4NDE4NEMxLjQzNDYyLTUuMjg0MTg0IC42Njk0ODktMy45NDUyMDUgLjY2OTQ4OS0yLjU3MDM2MUMuNjY5NDg5LTEuMjMxMzgyIDEuNDEwNzEgLjEzMTUwNyAyLjcyNTc3OCAuMTMxNTA3QzQuMjQ0MDg1IC4xMzE1MDcgNS4xMjg3NjctMS4zMjcwMjQgNS40MTU2OTEtMS45MjQ3ODJDNS44OTM4OTgtMS4xOTU1MTcgNi4wNjEyNy0uOTMyNTAzIDYuNTAzNjExLS41NjE4OTNDNy4yMzI4NzcgLjA3MTczMSA3Ljg5MDQxMSAuMTMxNTA3IDguMTg5MjkgLjEzMTUwN0M4LjMzMjc1MiAuMTMxNTA3IDguNTI0MDM1IC4xMDc1OTcgOC42MzE2MzEgLjA4MzY4NlYtLjM4MjU2NVpNNS4xNTI2NzctMi4zNDMyMTNDNC44NjU3NTMtMS42OTc2MzQgNC4wNjQ3NTctLjIxNTE5MyAyLjY0MjA5Mi0uMjE1MTkzQzEuNDgyNDQxLS4yMTUxOTMgLjkzMjUwMy0xLjQ4MjQ0MSAuOTMyNTAzLTIuNTcwMzYxQy45MzI1MDMtMy43NTM5MjMgMS42MDE5OTMtNC43OTQwMjIgMi41OTQyNzEtNC43OTQwMjJDMy41Mzg3My00Ljc5NDAyMiA0LjE3MjM1NC0zLjg2MTUxOSA0LjU3ODgyOS0zLjIzOTg1MUM0LjY5ODM4MS0zLjAzNjYxMyA1LjAyMTE3MS0yLjU0NjQ1MSA1LjE1MjY3Ny0yLjM0MzIxM1onLz4KPHBhdGggaWQ9J2cxLTUwJyBkPSdNNi41NTE0MzItMi43NDk2ODlDNi43NTQ2Ny0yLjc0OTY4OSA2Ljk2OTg2My0yLjc0OTY4OSA2Ljk2OTg2My0yLjk4ODc5MlM2Ljc1NDY3LTMuMjI3ODk1IDYuNTUxNDMyLTMuMjI3ODk1SDEuNDgyNDQxQzEuNjI1OTAzLTQuODI5ODg4IDMuMDAwNzQ3LTUuOTc3NTg0IDQuNjg2NDI2LTUuOTc3NTg0SDYuNTUxNDMyQzYuNzU0NjctNS45Nzc1ODQgNi45Njk4NjMtNS45Nzc1ODQgNi45Njk4NjMtNi4yMTY2ODdTNi43NTQ2Ny02LjQ1NTc5MSA2LjU1MTQzMi02LjQ1NTc5MUg0LjY2MjUxNkMyLjYxODE4Mi02LjQ1NTc5MSAuOTkyMjc5LTQuOTAxNjE5IC45OTIyNzktMi45ODg3OTJTMi42MTgxODIgLjQ3ODIwNyA0LjY2MjUxNiAuNDc4MjA3SDYuNTUxNDMyQzYuNzU0NjcgLjQ3ODIwNyA2Ljk2OTg2MyAuNDc4MjA3IDYuOTY5ODYzIC4yMzkxMDNTNi43NTQ2NyAwIDYuNTUxNDMyIDBINC42ODY0MjZDMy4wMDA3NDcgMCAxLjYyNTkwMy0xLjE0NzY5NiAxLjQ4MjQ0MS0yLjc0OTY4OUg2LjU1MTQzMlonLz4KPHBhdGggaWQ9J2cxLTU2JyBkPSdNNi41ODcyOTgtNy44NDI1OUM2LjY0NzA3My03Ljk3NDA5NyA2LjY0NzA3My03Ljk5ODAwNyA2LjY0NzA3My04LjA1Nzc4M0M2LjY0NzA3My04LjE3NzMzNSA2LjU1MTQzMi04LjI5Njg4NyA2LjQwNzk3LTguMjk2ODg3QzYuMjUyNTUzLTguMjk2ODg3IDYuMTgwODIyLTguMTUzNDI1IDYuMTMzMDAxLTguMDIxOTE4TDUuMTQwNzIyLTUuMzkxNzgxSDEuNTA2MzUxTC41MTQwNzItOC4wMjE5MThDLjQ1NDI5Ni04LjE4OTI5IC4zOTQ1MjEtOC4yOTY4ODcgLjIzOTEwMy04LjI5Njg4N0MuMTE5NTUyLTguMjk2ODg3IDAtOC4xNzczMzUgMC04LjA1Nzc4M0MwLTguMDMzODczIDAtOC4wMDk5NjMgLjA3MTczMS03Ljg0MjU5TDMuMDQ4NTY4LS4wMTE5NTVDMy4xMDgzNDQgLjE1NTQxNyAzLjE2ODEyIC4yNjMwMTQgMy4zMjM1MzcgLjI2MzAxNEMzLjQ5MDkwOSAuMjYzMDE0IDMuNTM4NzMgLjEzMTUwNyAzLjU4NjU1IC4wMTE5NTVMNi41ODcyOTgtNy44NDI1OVpNMS42OTc2MzQtNC45MTM1NzRINC45NDk0NEwzLjMyMzUzNy0uNjU3NTM0TDEuNjk3NjM0LTQuOTEzNTc0WicvPgo8L2RlZnM+CjxnIGlkPSdwYWdlMSc+Cjx1c2UgeD0nMTE3LjQ1MTkwNScgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzEtNTYnLz4KPHVzZSB4PScxMjQuMDkzNjg1JyB5PSctNC4yMjkwNzUnIHhsaW5rOmhyZWY9JyNnMC0xMjAnLz4KPHVzZSB4PScxMzUuMjkzMjk0JyB5PSctNC4yMjkwNzUnIHhsaW5rOmhyZWY9JyNnMS01MCcvPgo8dXNlIHg9JzE0Ni41ODQyNjInIHk9Jy00LjIyOTA3NScgeGxpbms6aHJlZj0nI2c0LTkxJy8+Cjx1c2UgeD0nMTQ5LjgzNTkyMycgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzEtMCcvPgo8dXNlIHg9JzE1OS4xMzQ0MicgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzQtNTEnLz4KPHVzZSB4PScxNjQuOTg3NDExJyB5PSctNC4yMjkwNzUnIHhsaW5rOmhyZWY9JyNnMi01OScvPgo8dXNlIHg9JzE3MC4yMzE1NycgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzQtNTEnLz4KPHVzZSB4PScxNzYuMDg0NTYnIHk9Jy00LjIyOTA3NScgeGxpbms6aHJlZj0nI2c0LTkzJy8+Cjx1c2UgeD0nMTc5LjMzNjIyMScgeT0nLTkuMTY1MjYxJyB4bGluazpocmVmPScjZzMtNTAnLz4KPHVzZSB4PScxOTUuNzc0NTE2JyB5PSctNC4yMjkwNzUnIHhsaW5rOmhyZWY9JyNnMi0xMDInLz4KPHVzZSB4PScyMDIuODIwOTUzJyB5PSctNC4yMjkwNzUnIHhsaW5rOmhyZWY9JyNnNC00MCcvPgo8dXNlIHg9JzIwNy4zNzMyNzknIHk9Jy00LjIyOTA3NScgeGxpbms6aHJlZj0nI2cwLTEyMCcvPgo8dXNlIHg9JzIxNS4yNTIwNTgnIHk9Jy00LjIyOTA3NScgeGxpbms6aHJlZj0nI2c0LTQxJy8+Cjx1c2UgeD0nMjIzLjEyNTIxNCcgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzEtNDcnLz4KPHVzZSB4PScyMzUuNzQ0NTQnIHk9Jy00LjIyOTA3NScgeGxpbms6aHJlZj0nI2cyLTk3Jy8+Cjx1c2UgeD0nMjQxLjg4OTQ4NCcgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzQtNDAnLz4KPHVzZSB4PScyNDYuNDQxODEnIHk9Jy00LjIyOTA3NScgeGxpbms6aHJlZj0nI2cwLTEyMCcvPgo8dXNlIHg9JzI1NC4zMjA1OScgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzQtNDEnLz4KPHVzZSB4PScyNTguODcyOTE2JyB5PSctOS4xNjUyNjEnIHhsaW5rOmhyZWY9JyNnMy00OScvPgo8dXNlIHg9JzI2My4xMDcwOTgnIHk9Jy05LjE2NTI2MScgeGxpbms6aHJlZj0nI2czLTQ4Jy8+Cjx1c2UgeD0nMjY3LjgzOTQxMycgeT0nLTQuMjI5MDc1JyB4bGluazpocmVmPScjZzItNTknLz4KPC9nPgo8L3N2Zz4KPCEtLSBERVBUSD0wIC0tPg==)
where 
is the Ackey function defined in The Ackley test case page.
In the following we call it the “Ackley distribution”.
import openturns as ot
import openturns.experimental as otexp
from openturns.viewer import View
from openturns.usecases import ackley_function
from numpy import exp, format_float_scientific
ot.RandomGenerator.SetSeed(100)
Prepare the Metropolis-Hastings algorithm¶
Define the Ackley distribution support and density (up to a constant factor).
am = ackley_function.AckleyModel()
ackley = am.model
power10 = ot.SymbolicFunction("x", "x^10")
ackley_pdf = ot.ComposedFunction(power10, ackley)
logarithm = ot.SymbolicFunction("x", "10 * log(x)")
ackley_logpdf = ot.ComposedFunction(logarithm, ackley)
lb = -3.0
ub = 3.0
support = ot.Interval([lb] * 2, [ub] * 2)
Define the proposal distribution as a Histogram.
Its ticks (on the X axis of the PDF of the histogram) will remain the same,
but its frequencies (on the Y axis) will be updated
during the course of the Metropolis-Hastings algorithm.
n_bins = 50
myticks = ot.RegularGrid(lb, (ub - lb) / n_bins, n_bins + 1).getValues()
frequencies = [1.0] * (myticks.getSize() - 1)
proposal = ot.Histogram(myticks, frequencies)
The state of the Markov chain must be converted to an acceptable set
of parameters for the Histogram distribution.
This is the job of the link function,
which we construct with the OpenTURNSPythonFunction class.
It takes a state as input and outputs the parameters (ticks and frequencies)
of the proposal distribution.
In our case, the ticks will not depend on the inputs,
but the frequencies will be outputs of the Ackley function.
parameter_dim = proposal.getParameter().getSize()
parameter_desc = proposal.getParameterDescription()
class ConditionalAckley(ot.OpenTURNSPythonFunction):
"""
When executed, this function returns the parameters of a Histogram
which approximates the conditional Ackley distribution obtained
when one of the 2 coordinates is fixed.
To compute the frequencies of the Histogram,
this OpenTURNSPythonFunction computes the values of the Ackley function
on a regular grid on the line parallel to
either the (1, 0) vector (if the second coordinate is fixed)
or the (0, 1) vector (if the first coordinate is fixed)
containing the point passed as input.
The regular grid covers the part of this line which is contained in
the support of the Ackley distribution, implicitly defined as the
smallest square that contains the cartesian product of the regular grid
with itself. For example, if the regular grid covers the interval [-3, 3],
then the support is the square [-3, 3] x [-3, 3].
Parameters
----------
marginal : int
The marginal whose value is *not* fixed.
If 0, then the line of the regular grid is parallel to the (1, 0) vector.
If 1, then the line of the regular grid is parallel to the (0, 1) vector.
ticks : RegularGrid
Ticks of the Histogram distribution.
"""
def __init__(self, marginal, ticks):
super().__init__(2, parameter_dim)
self.setInputDescription(["X0", "X1"]) # input: X0 and X1 coordinates
self.setOutputDescription(parameter_desc) # output: parameters of the Histogram
self._marginal = marginal # parameter which does not vary after initialization
offset = (ticks[1] - ticks[0]) / 2
self._marginal_inputs = ot.Sample.BuildFromPoint(ticks)[0:-1] + offset
# _marginal_inputs contains the varying coordinate of the points in the regular grid
self._size = self._marginal_inputs.getSize()
self._ticks = ticks
def _exec(self, X):
"""
Execute the function on a point X = (X0, X1).
Parameters
----------
X : list of 2 floats
Point through which the line containing the regular grid passes.
Returns
-------
parameters : :class:`~openturns.Point`
Parameters of the :class:`~openturns.Histogram`.
"""
inputs = ot.Sample(self._size, X) # sample of inputs for the Ackley function
# All input points are initialized at the point X passed as argument.
# Replace the varying coordinate with the values of the regular grid.
inputs[:, self._marginal] = self._marginal_inputs
# Compute the Ackley function on these inputs.
outputs = exp(ackley_logpdf(inputs).asPoint())
# The outputs are the unnormalized frequencies of the Histogram
# proposal distribution, but the Histogram.setParameter() method
# expects a full set of parameters.
# The easiest way to provide it is to construct a new Histogram object
# with the adequate frequencies and call its getParameter() method.
return ot.Histogram(self._ticks, outputs).getParameter()
The 2 components of the state of the Markov chain will be updated
one after the other, not simultaneously.
We define 2 UserDefinedMetropolisHastings algorithms
encapsulated within a Gibbs algorithm,
so we need 2 link functions, each corresponding to one of the marginals
of the Ackley distribution.
Note that thanks to the OpenTURNSPythonFunction class,
we were able to only code one template to be used by two different functions
instead of directly coding two functions with the PythonFunction class.
link_function_0 = ot.Function(ConditionalAckley(0, myticks))
link_function_1 = ot.Function(ConditionalAckley(1, myticks))
Let us illustrate the first of these functions.
We can start by evaluating it at 
.
Let 
be a bidimensional random variable
following the Ackley distribution.
The output is the set of parameters for a Histogram
distribution that approximates the distribution of 
conditional on 
.
Let us update the Histogram we created before with this set of parameters.
par = link_function_0([0.5, 1.5])
proposal.setParameter(par)
print(par)
[-3,0.12,0.574885,0.12,0.590611,0.12,0.616354,0.12,0.575301,0.12,0.474915,0.12,0.345447,0.12,0.217746,0.12,0.127608,0.12,0.0892587,0.12,0.086088,0.12,0.093021,0.12,0.0909145,0.12,0.0769473,0.12,0.0565498,0.12,0.0353399,0.12,0.0193349,0.12,0.0118315,0.12,0.01083,0.12,0.0129736,0.12,0.0148585,0.12,0.0146287,0.12,0.0124841,0.12,0.0091527,0.12,0.00577608,0.12,0.00381056,0.12,0.00381056,0.12,0.00577608,0.12,0.0091527,0.12,0.0124841,0.12,0.0146287,0.12,0.0148585,0.12,0.0129736,0.12,0.01083,0.12,0.0118315,0.12,0.0193349,0.12,0.0353399,0.12,0.0565498,0.12,0.0769473,0.12,0.0909145,0.12,0.093021,0.12,0.086088,0.12,0.0892587,0.12,0.127608,0.12,0.217746,0.12,0.345447,0.12,0.474915,0.12,0.575301,0.12,0.616354,0.12,0.590611,0.12,0.574885]#101
The link_function_0 function computes histogram frequencies proportional
to the values of the unnormalized Ackley distribution PDF
along the line.
title = "Ackley PDF (up to a constant factor) and $X_1 = 1.5$ cross-section"
graph = ot.Graph(title, "$X_0$", "$X_1$", True)
line = ot.Curve([[-3, 1.5], [3, 1.5]], "black", "dashed", 2)
graph.add(line)
graph.setLegendPosition("topright")
contour = ackley_pdf.draw([lb] * 2, [ub] * 2)
reversed_colors = [color for color in reversed(contour.getColors())]
contour.setColors(reversed_colors)
legends = contour.getLegends()
legends = [format_float_scientific(float(v), precision=1) for v in legends[:-1]]
contour.setLegends(legends)
graph.add(contour)
view = View(graph, legend_kw={"bbox_to_anchor": (1, 1), "loc": "upper left"})
view.getFigure().tight_layout()
Let us compare the histogram to the unnormalized PDF (we still need to rescale it to make it appear properly on the graph).
scaling = ot.SymbolicFunction("x", "1.3e-10 * x")
scaled_ackley_pdf = ot.ComposedFunction(scaling, ackley_pdf)
graph = proposal.drawPDF()
graph.setXTitle("")
graph.add(scaled_ackley_pdf.draw(0, 0, [0.0, 1.5], -3.0, 3.0, 100))
graph.setLegends(["Histogram PDF", "Rescaled unnormalized PDF"])
graph.setLegendPosition("top")
graph.setTitle("Conditional distribution of $X_0$ given $X_1 = 1.5$")
_ = View(graph)
Let us now do the same think with link_function_1.
This time, the relevant cross-section is along the line 
.
par = link_function_1([0.5, 1.5])
proposal.setParameter(par)
print(par)
title = "Ackley PDF (up to a constant factor) and $X_0 = 0.5$ cross-section"
graph = ot.Graph(title, "$X_0$", "$X_1$", True)
line = ot.Curve([[0.5, -3], [0.5, 3]], "black", "dashed", 2)
graph.add(line)
graph.setLegendPosition("topright")
contour = ackley_pdf.draw([lb] * 2, [ub] * 2)
reversed_colors = [color for color in reversed(contour.getColors())]
contour.setColors(reversed_colors)
legends = contour.getLegends()
legends = [format_float_scientific(float(v), precision=1) for v in legends[:-1]]
contour.setLegends(legends)
graph.add(contour)
view = View(graph, legend_kw={"bbox_to_anchor": (1, 1), "loc": "upper left"})
view.getFigure().tight_layout()
scaling = ot.SymbolicFunction("x", "3.1e-10 * x")
scaled_ackley_pdf = ot.ComposedFunction(scaling, ackley_pdf)
graph = proposal.drawPDF()
graph.setXTitle("")
graph.add(scaled_ackley_pdf.draw(1, 0, [0.5, 0.0], -3.0, 3.0, 100))
graph.setLegends(["Histogram PDF", "Rescaled unnormalized PDF"])
graph.setLegendPosition("top")
graph.setTitle("Conditional distribution of $X_1$ given $X_0 = 0.5$")
_ = View(graph)
[-3,0.12,0.714193,0.12,0.709153,0.12,0.714195,0.12,0.635879,0.12,0.49337,0.12,0.331163,0.12,0.187923,0.12,0.0968651,0.12,0.0598784,0.12,0.0526439,0.12,0.0522747,0.12,0.0457978,0.12,0.0334477,0.12,0.0202285,0.12,0.00973413,0.12,0.00378105,0.12,0.00159728,0.12,0.0011016,0.12,0.001075,0.12,0.00097469,0.12,0.00070403,0.12,0.00040715,0.12,0.0001854,0.12,6.70262e-05,0.12,2.77594e-05,0.12,2.77594e-05,0.12,6.70262e-05,0.12,0.0001854,0.12,0.00040715,0.12,0.00070403,0.12,0.00097469,0.12,0.001075,0.12,0.0011016,0.12,0.00159728,0.12,0.00378105,0.12,0.00973413,0.12,0.0202285,0.12,0.0334477,0.12,0.0457978,0.12,0.0522747,0.12,0.0526439,0.12,0.0598784,0.12,0.0968651,0.12,0.187923,0.12,0.331163,0.12,0.49337,0.12,0.635879,0.12,0.714195,0.12,0.709153,0.12,0.714193]#101
We now choose the initial state of the Markov chain.
initialState = [0.1] * 2
Sample from the Ackley distribution¶
We can finally define the two Metropolis-Hastings algorithms and the Gibbs algorithm which encapsulates them.
gmh_0 = otexp.UserDefinedMetropolisHastings(
ackley_logpdf, support, initialState, proposal, link_function_0, [0]
)
gmh_1 = otexp.UserDefinedMetropolisHastings(
ackley_logpdf, support, initialState, proposal, link_function_1, [1]
)
gibbs = ot.Gibbs([gmh_0, gmh_1])
sample = gibbs.getSample(100)
Print the acceptance rates of the two Metropolis-Hastings samplers
mhlist = gibbs.getMetropolisHastingsCollection()
rate_gmh_0 = mhlist[0].getAcceptanceRate()
rate_gmh_1 = mhlist[1].getAcceptanceRate()
print("Acceptance rates: {} and {}".format(rate_gmh_0, rate_gmh_1))
Acceptance rates: 0.99 and 0.91
Draw the contour plot of the Ackley distribution PDF (up to a multiplicative constant) and the MCMC sample: we can see the sample is accurate.
title = "Ackley PDF (up to a constant factor) and Metropolis-Hastings sample"
graph = ot.Graph(title, "$X_0$", "$X_1$", True)
graph.add(ackley_pdf.draw([lb] * 2, [ub] * 2))
reversed_colors = [color for color in reversed(graph.getColors())]
graph.setColors(reversed_colors)
graph.add(ot.Cloud(sample, "black", "plus"))
legends = graph.getLegends()
legends = [format_float_scientific(float(v), precision=1) for v in legends[:-1]]
graph.setLegends(legends)
graph.setLegendPosition("topright")
view = View(graph, legend_kw={"bbox_to_anchor": (1, 1), "loc": "upper left"})
view.getFigure().tight_layout()
