Welcome back everybody. We're going to hit the major transition point in the course today where we switch algorithms. Our algorithms for approximating the posterior distribution. Before we get into the details, let me remind you really what the heart of Bayesian inference is. Bayesian inference is not about how you compute the posterior distribution. It's just about the posterior distribution, however you get it. There are lots of valid ways to get an approximation. Remember the globe tossing from the beginning? In this case, you can calculate this exact posterior a huge number of ways and they're all valid. That's fine. I emphasize this because often in applied statistics, especially in biology, Bayes is thought of as using Markov chains. They're the same thing. You use the Markov chain as Bayes. Actually, they're not the same thing. You can use a Markov chain and not have a Bayesian model. Actually, the Markov chain strategy I'm going to show you is essentially frequentist. There's this weird thing about it, but it's used to calculate a posterior distribution. You can use a Markov chain to do lots of different things. The reason in this course that Markov chains come so late is partly to emphasize that to you that there's an essential thing about Bayesian inference, but it's not how you compute the posterior. It's just that you're after the posterior, however you want to get it. The second thing was to make it easier. There's a bunch of additional machinery you need to fool around with once you start running Markov chains. Today, I'm going to introduce you to that machinery and start to train you to be expert users so that you can make good decisions and recognize when the machine malfunctions. Remind you, in this course we've discussed four different, so far only three, but four different ways to compute a posterior distribution. There are others. First, there's this analytical approach, which I did not show you, but I used to make the graph on the previous slide. For very simple models, you can just get a closed form analytical expression for the posterior distribution. It's great. You take MAS stats courses, you'll learn that sort of thing. But for interesting models and applied stats, this is nearly always impossible, so we don't do it. So you need some other numerical way to do the integration. This is all about doing numerical integration. And I taught you grid approximation in the beginning, partly because I wanted it to reinforce the concept of what the posterior distribution is. It's just a thing that's proportional to the product of the probability of the data times the prior probability. That's all it is. That's the posterior. And then you can use in quadratic approximation with this funny tool named Quap for months now. Yeah, you've been loving it. Natalia's eyes rolled very significantly there. Sorry, I didn't mean to single you out. Lots of eyes rolled. But quadratic approximation is in this course instead of starting with Markov chains because I also want to... Well, first of all, it's very effective. It's unreasonably effective for a large number of very simple models. And it's also numerically the same thing as the way lots of non-Basian models are fit, estimated. If you're doing maximum likelihood estimation and getting some sort of symmetric standard error around each estimate, then you're doing the same mathematical steps. Even though people tell you it's philosophically completely different for some reason, but the numbers will be the same. Yeah. And I want you to see that connection. That there's a connection there between lots of standard tools and very basic Basian models. But now we're going to get into things where that connection gets a lot blurrier. So before getting into the details of what Markov chains are, let's have a parable, okay?  King Markov and you rule an island kingdom like all Polynesian kings did. And it's good to be king, but it also comes with obligations. There are certain things you must do to avoid being killed by your people. It's like all kings. Royalty is fragile as an institution, yeah. And so you rule over a kingdom called the Metropolis Archipelago and it has ten islands. There aren't ten here, but there'll be ten in a moment. And the islands are numbered and in order and they have different population sizes on them, different population densities. And the contract you have with your people is they'll let you continue to tax them mercilessly as long as you visit them. They love you. Well, you're a descendant of a god, right? The Hawaiian kings. And so as long as you visit them and they can bask in your godly glow, then they'll be happy. So you must visit them in proportion to their population density. So this is your goal. You're not numerate, not highly educated. You're king after all. No, sorry.  But a Polynesian king. And so you need some simple algorithm to do this that doesn't involve having a reliable census and a schedule and calendar for the whole year.  And here it is.  And now you're ready to move to another island. And so the first thing you do is you flip a coin. And this coin determines whether you go left or right. Well, it doesn't matter if you go left or right. This generates what's called a proposal. The proposal island, you're going to consider going to one of the islands, either to the left or to the right. You get this randomly from an unbiased coin toss or a shell toss or whatever you want. Think of it as. Then step two, you send a servant in a boat across to the proposal island. And they take a rough census of the population density on that island. You have another servant do the same for the island you're currently on. So you get the population of the proposal island from this. And we're going to call that number of the population of the proposal island little p sub five for the fifth island. You're also going to get the population of the current island. So this is p sub four. You want to compare these two numbers in a particular way. What you want to do is construct the ratio of them, p five or p four. And then this should be your probability of accepting the proposal of moving from island four to island five in this case is p five over p four. Now you're saying that p five over p four can be greater than one. Yes, it can. In that case, you'll move. Right, for sure. But you might stay. It depends. If p five is less than p four, then you might stay. I'm not guaranteed to stay. So imagine that p five is half of p four. Then you've got some chance, but you've got a good chance of moving and a good chance of staying. So you move except move to the proposal island with that probability. And then finally, then you get either you stay and then you do another tour of the island you're currently on.  And then you just repeat this starting over with stage one.  And this is this is the life of the politician king. The thing about this is that this is a valid way to fill the contract. This algorithm guarantees this is provable mathematically guarantees that in the long run, you will visit each island in proportion to its relative population size and fulfill your royal contract and not be assassinated by your people. At least not for this reason. This is an example of Markov chain Monte Carlo. In fact, as I'll explain in a second, it's the most famous and most primitive algorithm of Markov chain Monte Carlo. And what I want to do today in this lecture is help you understand the general strategies of these algorithms. Why we would use them because it seems mad, right? I hope you felt that that parable on the previous slide was crazy. It's just sort of madness. You can think of lots of better ways to do something like that. But it has a huge advantage if you don't know the distribution. So if you don't know the distribution of population sizes on the islands, you don't need to actually visit each of the islands in proportion to its population size. And those are the kinds of problems that we have when we're trying to find posterior distributions. We don't know the posterior distribution and yet we can visit each part of it in proportion to its relative probability using an algorithm like this. And that's the magic. We can sample from a distribution that we don't know. That's what we're doing. And this is just a tool for doing numerical integration. So I also want to introduce you to the interfaces to Markov chain Monte Carlo that I'm going to use in the rest of the course. The engine is something called STAN. I'll tell you what it is when we get there. STAN, it's not an acronym, it doesn't stand for anything. It's not in all caps. It's a person's name, Stanislav Ulam. I'll show you a picture of it later. And the tool that I wrote in rethinking called Ulam, which is a simplified input. It uses formulas exactly like the ones you've been using so far, but it'll make a custom Markov chain out of it for you. So I'm going to show you how that works too. Then I need to show you how to sample responsibly from this and diagnose when things go wrong. And you'll continue to practice this and develop lifelong skills and a love for Markov chains. Okay, here is the R script version of King Monte's Royal Tour. This is a very short and heavily commented chain. It's called the metropolis algorithm. Metropolis was also a person. I'll show you a picture of it in a minute. I don't usually step through R code in this course, but if you'll indulge me, I'd like to stay on this life and run and walk through each line for you to have a sense of how simple this algorithm is and that this code instantiates that weird story that I gave you on the previous slides.  It really is super simple sort of thing to do. Okay, at the top, we define how many weeks King Monte is going to do a tour. And in this case, it's past his lifespan probably. But 10 to the fifth. Then we do some record keeping. Positions is an empty vector full of zeros with the length and number of weeks. We're just going to store which island the King is on in this vector. In R, you always want to initialize all of your storage space at the beginning instead of growing it as you go because this is much faster this way. And then I'm going to just put him on an island on Island 10 to start. You could randomize that if you want. That's what current is. It's just the current position. And then we loop over the number of weeks. And the first step is we just record in the record keeping position specter where the King is now. That's that first line of positions. I is current. And now we're going to flip the coin to generate a proposal. The proposal is the current plus a sample from minus one or one, size one. So this is just a clever our way of adding minus one or one to work in your current position. That gives us a proposal. Island to go. And now we've got to handle the looping. So in the book I've got a picture of this archipelago or badly drawn picture of this archipelago. And it's in a ring. Imagine it's a collapsed volcano or something like that. Super volcano. And so it's in a ring and so you can loop from island 10 to island one the other way. So these two lines here are not necessarily essential but they just make it a loop, a donut of islands. And then finally there's the action, the magic part. We construct this ratio of the proposal to the current and I've got the island's number such that their population sizes, relative population size are their numbers. That's what I'm asserting here. When you do this with a stat model you have probabilities of data in there. That's what changes. And then I do a random number if that random number is less than the probability of the moving then we go to the proposal otherwise we stay in current.  This is what it looks like. You run that code and you should it's the safety of your own home. Run that code and plot it out. You get this sort of thing and what we're looking at is time moving from left to right and then the king's position on the vertical and you'll see it zigzags around and it loops from top to bottom in various places. There are places where it stays on the same island for a while. Because this is possible too. You get stuck on very densely populated islands like island 10. You need to stay there for weeks and weeks and weeks and kiss the same babies over and over again. There's lots of babies, right? Dense islands. You need to stay there longer. That's the whole point. But then you can take a tour. He'll flip over the small islands and only spend a week there. Those sorts of things. And in the long run when we construct a density out of these data points it's in the right proportions. So the way we usually think of chains is you think about connecting these with lines. So it's not just a bunch of points and you get this zigzag sort of pattern of movement in space. So in the long run let me give you an idea about the evolution of the frequency properties here. These bar graphs, these histograms show you different points moving left to right on this graph. On the left it's only after 100 weeks. You can see it's starting to emerge. He is visiting island 10 more than the others. But it's not. What we want is a perfect staircase. That would be the right. That's the Royal Contract. So in the middle after 500 weeks it's looking better. After 2000 weeks we're almost there. That's a lot of weeks though. I don't know how many weeks are in a year. Not that many. This is a lot of years. So this is King Monty's like great grandson has already taken over the throne and picked up the Royal Contract again. But it's still a valid Markov chain. It keeps going. This is an important thing about Markov chains. They guaranteed to work but only in the long run. And what the long run means is a very sensitive issue. And this will be something the reason we're not going to use this algorithm to get estimates of the posterior distribution we're going to use a better one that is much more efficient than this. You can do a lot better than this. But this works. Guaranteed. It's just that you might be dead by the time it delivers on its guarantee. So there's my summary side of that. This algorithm, the metropolis algorithm it will converge to the correct proportions in the long run. But that's a kind of vague promise. That's the kind of promise you get from math. It's like, yeah, I guarantee it, eventually it'll happen before the universe implodes. It doesn't matter where you start. It's not sensitive to initial conditions. And for this exact algorithm you need symmetric proposals. Your probability of going one direction needs to be the same as going the other. Well, the proposals. Your probability of going one direction and going the other is actually not the same because of the relative population sizes. But the proposal time needs to be fair. There are other algorithms where that's not required. You can extend this to something called metropolis hastings. I mentioned this in the book notes, where the proposals can have any distribution you like.  You can make more intelligent proposals. So this is called the metropolis algorithm. And we're not going to actually use it, as you might imagine, to tour archipelagos. In our applications, the islands are parameter values. And there can be an infinite number of them. They don't have to even be discrete. You can just have a continuous smear of islands and you can jump among them. And the population size is the posterior probability. That's what we want to sample from. Because there's this distribution, which is the posterior. And we need to visit parts of it in proportion to their height, local height, to the posterior probability in that location. This works no matter how many parameters you have in the long run. This is going to be a big snag in a moment. If you get a lot of parameters, the long run is really long. But we have a solution. Nope, Fred. And it works for continuous and discrete parameters. So the reason it's called metropolis is because Nick Metropolis was a Greek-American mathematician who headed a research group at Los Alamos that developed the first computer application of this algorithm on a computer called Maniac. Here's the famous paper where they describe it, a really terrible title, right? Equation of state calculations by fast computing machines. Fast computing has got some irony if you read this now, right? Every one of you has a phone in your pocket that is vastly more powerful than this computer was. But this was really an extraordinary achievement at the time. And some famous people here, including famous physicists, Ariana and Marshall Rosenbluth, Gusta Teller and Edward Teller. So I think it was Ariana Rosenbluth who did most of the programming. And I think there's a story where Gusta Teller figured out the memory. They got a bicycle wheel and a bunch of magnetic tape and they spun the bicycle wheel to feed the tape in and out of the machine. There was no, like, store where you go buy a disc or a memory stick or something, right? They did this all from scratch. And so there's a great book. I forget who wrote it right now. I'm sorry, where they tell all these stories about how they invented all this stuff. It made this thing work. It's really an amazing thing. But I should say, though, this was all part of the Manhattan Project. It was this post-war, but this was making fusion bombs, right? They were running the invented disc algorithms so they could simulate fusion of, you know, a bomb explosion. That's what it was for. This is Teller, Ed Teller. Some of you know this history. So there's some monstrous stuff here. We're not going to do that. We're not going to do simulated fusion explosions. But you could if you wanted. But please don't. All right. Here's Nick Metropolis in the foreground smoking because everybody smoked. Everywhere. And that's maniac in the background, part of it. He was distributed across rooms and there's all these things, you know, making tubes and everything. Let me give you a sense of it. Maniac is an acronym that Nick Metropolis came up with because already acronyms have gotten ridiculous by then and he was like, okay, I think I've got to put it into this. Let's make a ridiculous acronym. And it didn't do any good, did it? Acronyms are still bad. But so it's the mathematical analyzer, numerical integrator, and computer. And I think there's also an irony here because many of the members of the Manhattan Project had moral doubts about what they were doing. Right? So this encrypted stuff that they did like this. So Maniac was a thousand pounds and had five k bytes of memory, which is basically nothing. That's not one meme. Right? You can't fit one cat meme in five k bytes of memory. Right? And 70 k multiplications per second, which sounds fast compared to what I could do with a pencil. But compared to your laptop or your mobile phone, your laptop, a really cheap Barbie-based laptop, would be four to seven pounds, two to eight million kilobytes of storage, and could do billions of multiplications per second. That's super cheap. And your phone as well could do that easily. So we've come a long way in computing power. But the fundamental algorithms and the mathematics that justifies those algorithms and guarantees they work is the same. It was developed during this time that computing was really difficult. So what are Markov chains? So the Metropolis algorithm is the simplest version of what's called Markov chain Monte Carlo. What's a chain? A chain is the sequence of moves. That's why it's called a chain. There's a sequential simulation that's being run. Then we call that a chain. In particular, this is a Markov chain, which is named after this fellow, Andrey Andreevich Markov, who is a very famous Russian mathematician. He did a lot of foundational work in probability theory. And he didn't call these things Markov chains himself. Other people did, but he did all this work on probability systems that have no memory. So what makes something a Markov chain is that it's a sequential simulation where the probability you go to any new place in the chain depends only on where you are now. It doesn't depend upon where you've been in the past. So that's why we say they're memoryless systems. What matters is the current state and not past states. Does that make sense? And that's so in King Monte's tour all that matters is where he is now. It doesn't matter where he was last week. And that's why it's a Markov chain. Markov chains are great for computing because you don't have to store a bunch of numbers. So if you can translate a system into a memoryless version of itself, then you can call it a Markov chain and then there are a bunch of great things, tricks you can use to calculate things about it. And this fellow, Andre, did the most foundational work in this. And then the Monte Carlo is named after the city where people gamble. You may have heard of it because it's random. It's stochastic. That's not what the Monte Carlo card is. There wasn't a person named Monte Carlo. It's really there was at some point. Nothing to do with this story. Okay. So why do we use these crazy things? I said a little bit about this before. Sometimes you get a model that's too complicated to have an analytical expression from the integrated posterior year. Nobody can do it. Integrals are hard. Those of you who've done much integral calculus know that with integrals you're usually guessing. There are nice ones where there are techniques that work. You can integrate by parts or something like that.  And if you could solve them, you'd be famous. It's not like derivatives. Derivatives, you just give it to your computer. Those are the derivative. They're just a rule. But going the other direction is a lot harder. And so our goal here, you get the posterior distribution by integrating. And integration is sometimes just not practical. So we do numerical integration instead. We figure out the integral through numerical techniques. And Mark Auté Monte Carlo is just one of those. And it's a very, very capable one that works for a very large class of problems. So we're going to be interested in problems like multi-level models if you're interested in networks. Network models have intractable likelihoods as well. Phylogenies. Mark Auté Monte Carlo's very big in phylogenetic inference for this reason as well. And spatial models also have lots of really difficult things about them. Optimization, which is what CRAP does, is in general a non-starter in high dimensions. Once you've got a lot of parameters, which we're going to call dimensions, optimization actually targets the wrong part of the distribution. I say more about this in the text. I don't think I'll have time today to talk about that more. There's this phenomenon I described in the textbook called concentration of measure, which guarantees that in a high-dimensional distribution the mode is very, very far from where all the probability mass is. And it's just a super weird thing. I have a whole section in the book about it. I might have time today to talk about it, but I don't think I do. You just have to realize this is why everybody gives up on optimization as soon as you get 100 parameters, basically. Optimization is giving you the wrong answer. It's giving you the right answer. You just ask the wrong question. The mode is now not relevant. Bottom line, the bold on this slide. Mark Auté Monte Carlo is not fancy. It's old. It's from the 1950s. This is middle 20th century technology. It's not some fancy thing. I say this because in some fields people say, oh, you're running a markup chain. Here's your fancy statistician now. No, this is from the 50s. This is old stuff. This is before people knew that smoking caused cancer. This is old technology. It's robust. It's used for a ton of applications, routine applications all the time. There's nothing fancy or shown off about markup chains. I want you to see it. It's just bread and butter sort of stuff in applied statistics. And it is. So there's a lot of text on this slide. I apologize. I just want to give you an overview before I get into some more interesting visualizations.  The thing I showed you is Metropolis. You can extend Metropolis very directly to something called Metropolis Hastings. Hastings was also a person. And Hastings published a very clear and famous paper that showed that you don't have to have symmetric proposals. You can choose more intelligent asymmetric proposal distributions. And you can get more efficient tours that converge faster to the target distribution. And that's a big achievement. And Metropolis Hastings has been customized to a bunch of really clever strategies, including the most famous of which is Gibbs sampling. Gibbs was a physicist who didn't work on this, but it was named in honor of him. Gibbs worked on stochastic distributions, statistical mechanics. And so they named this thing after it posthumously. Gibbs sampling in low dimensional models at least is extremely efficient because it uses a very clever proposal distribution. But it's basically a Metropolis algorithm. Same thing. What all these algorithms have in common is they guess and check. They make a proposal of someplace to move to. They check the posterior probability at that location and they compare it to where they are now. This is the guess and check strategy. And all of the Metropolis descended algorithms do that. As a consequence, the quality of your proposals is the most important thing because if you make dumb proposals, you'll reject them all.  The goal is to accept every proposal. If you could constantly be moving, you would tour the space as fast as possible if you could make really smart proposals.  So it's very hard to say how you're going to make smart proposals so that they're all accepted if you don't even know the distribution. Well, happy news. I'm going to show you a way today to do that. But you need a completely different strategy. And the strategy is known as Hamiltonian Monte Carlo. It does not guess and check. It's super cool. And this is what we're going to use in this class. And this is the new hotness. Can I still use that phrase? That's sort of old now. It's not even new. But now everybody is switching to using Hamiltonian Monte Carlo because it's vastly more efficient. And you can do really hard to mention models with it. So as I mentioned occasionally, wrapping up a project where there are 27,000 parameters and we just run it as a Hamiltonian chain and there's no problem. It just samples them all. No problem. I mean, it takes two days to run the chain, but that's no big deal, right? It took 15 years to collect the data. So it's no big deal. Okay. And now I say there are lots of new methods being developed all the time. But all of them have in common with Hamiltonian Monte Carlo is they use something called a gradient. They don't guess and check. The error of guess and check is kind of over. You can still use all those algorithms for simple problems. That's great. There's no problem. But you get a lot of efficiency out of avoiding guess and check. So let's see the problem of guess and check. So I want to show you a quick movie here. This is going to be the Metropolis Hastings Algorithm. In a two-dimensional Gaussian posterior hill that we're looking down on. And this is a simulation of the Markov chain running. So this is like King Monty, but now it's on a hill and he needs to tour this hill in proportion to its height. He needs to spend more time in the middle where it's high and less time in the bottom. So we're in the middle there. The arrows are proposals, green arrows are acceptances, red arrows are rejections. So when there's a red arrow, it stays. When there's a green arrow, it moves. We tour around every time and we record our position at every point and we stack up on these histograms and you see it's building the distribution. It's a two-dimensional Gaussian. You see this? Yeah. So these are up on my blog, by the way. If you want to play with these, you can change the animations and play around with it. I got it from this library. I'll give you the link at the bottom. There's a fellow who's written a giant library of these interactive animations of Markov chains. They're fun. So this is the Metropolis Algorithm working. Works really well, right? It's still an amazing achievement in mathematics that this sort of thing functions at all. The problem with the Metropolis Hastings Algorithm comes when the distribution's not so nice. So here's a distribution, which is not a nice Gaussian hill. It's a donut. And you're thinking, like, when am I going to have a donut distribution? Oh, you will. Oh, you will. You won't realize it. But in high-dimensional spaces, basically, as I said, the probability gets concentrated into a thin shell, into a thin region. It's very far from the middle, from the mode. And so it becomes this hyper-donut, if you will. That's probably the first time that word has ever been used. But picture a hyper-donut. And high-dimensional donut. 27,000-dimensional donut. And you must sample from that 27,000-dimensional donut. The Metropolis Algorithm is really bad at this because it makes proposals outside the high probability region everywhere. So let me show you how this animation goes. So you see Metropolis. It's making proposals. Most of them are rejected. It's very hard for it to move because it's making proposals into dead space all the time. And it doesn't know there's a donut here. It's got to map it out, right? You see the donut, but it doesn't, right? And your real problems, you don't know what the shape looks like either. So in the long run, this will still work. But the long run now is very long indeed, right? You need to finish your PhD.  So I'll let this run for a little bit. I think at some point I sped up this, yeah, speed it up. And you see, eventually it's going to work, but it hasn't even been to the south side yet. It's getting there. All right, yay, visit the south side. And it just gets stuck in a little narrow region because it's making bad proposals all the time. The proposals don't know anything about the target distribution. And that's the basic problem with Metropolis Hastings. Gibb sampling has the same issues. In a high-dimensional space, gibb sampling, it makes better proposals, but it still gets stuck because there's this concentration of the measure of the distribution. So let me show you the contrast here. You can tune the Metropolis family algorithm so they do better in a particular case, but there's this really tight trade-off. You're always given up on something. So I'll show you, this is the Metropolis algorithm. This is from this figure 9.3 in the book. And on the left, I'm showing you a distribution. This is two-dimensional Gaussian, but now there's a correlation in the posterior here. So this is negative .9 is the correlation between these two parameters. And we start in the upper left with the Markov chain. Field circles are accepted proposals. Open circles are rejected ones. And I made the step size really small. What's the step size? The step size is how wide the proposal distribution is. So if you have a very small proposal distribution, then you're only going to consider points near you. That means that a lot of them will be accepted because if you're in valid space, and you are, then you're going to make proposals in valid space. But as a consequence, you move very, very slow through the posterior distribution. So the acceptance rate is about 60%. We're rejecting 40%. We're going to do a lot better than this in a moment. Hang on. So we're moving. We're accepting most of them at least, but it takes us a very long time to get all the way down this valley, this pickle, whatever you want to call that shape. And on the right, it's the same posterior distribution, the same algorithm, but I increased the step size because we want to take a faster tour. And we do in a sense, but now we're rejecting the majority of proposals because they're outside the narrow space that we need to stay in. This is the fundamental trade-off, excuse me, in Metropolis. And really, you can't win. You've got to pick your poison on this. And again, Gibbs can do a little bit better, but as soon as the dimension of the posterior gets high, this problem always appears again. The problem is guessing and checking. So what do we do instead? Hamiltonian Monte Carlo has a different approach entirely. It does not guess and check. Instead, it runs a physics simulation. This is crazy, but bear with me. It really does. I'm not joking. This is not a metaphor. It actually runs a physics simulation. And what we're going to do is we're going to represent our parameter state. That is, there's a vector of parameter values in your model. And we're going to represent that as a coordinate in some high-dimensional space. So in a two-dimensional posterior distribution, it's just your x and y coordinates. In a three-dimensional one, it's your x, y, and z. And in more dimensions, you've got some hyper space that you have a coordinate position in. And now you're going to say you're a particle in this space and you have that position.  We're going to play God. We're going to flick this particle with some random momentum in a random direction in this high-dimensional space. And it's going to cruise on a surface. What's the surface? The surface is the posterior distribution. And it's going to cruise around on the posterior distribution. We're going to let the simulation run for a little while. And then we're going to stop the particle and we're going to record where it is. Then we're going to flick it in a new random direction where it is. And because it's following the curvature, I'm going to show you movies of this, so hang on. Because it's following the curvature, it always makes good proposals because it doesn't go into bad areas. Why not? Because you won't go into those areas because of the shape of the surface you're on. So there's no more guessing and checking. All proposals are good proposals. It's like dogs, Brent, are all good dogs. And we make really good proposals. Let me show you what this looks like. One more allegory. Imagine King Monty's... Now not King Markov, but King Markov's brother on the mainland is King Monty. He has a kingdom, but it's in a valley. It's a narrow valley. This is a north-south dimension. And there are mountains on the north end and the south end. Everybody just lives in this valley, is strung out along in a straight line. It's a continuous urban smear. And there are more people living in the bottom than on the edges. Same royal contract. You need to visit them in proportion to their population density. How can you do this? We can use the Hamiltonian Monte Carlo strategy to do it. How do you do it? He gets in the royal vehicle, starts driving randomly north or south at some random speed. The car obeys physics. It has to. If it starts going uphill, it slows down. You've got your foot at some constant pressure on the throttle. That's the idea.  Eventually you stop and then turn around. You go downhill. That same constant throttle. So you speed up when you're going downhill. You slow down when you're going uphill. And you drive for some sort of pre-specified duration. Then you stop. You get out.  Get back in the car. Repeat. This is the Hamiltonian Monte Carlo algorithm in silly story four. And if you do all this right, the stopping positions will be proportional to the population density at each point in the valley. In the book, I'll show you the simulation of this. This is what it looks like. We have time on the horizontal. On the vertical is the position. This is just a one-dimensional distribution. A long narrow valley with mountains on each end. So it's Gaussian.  Gaussian is not a valley. If you take the long government of Gaussian, it is. It's a parabola. So the minus long probability of a Gaussian is a valley. It's a parabolic valley. And that's what you end up with with a Gaussian. And so it starts in the middle. Randomly drives south. Starts to go uphill. He slows down in the car. He stops. And it turns around. He goes down. Then he stops. His timer goes off. Kisses some babies. Again, randomly drives uphill. Comes back down. Stops in the middle. Goes to the next point. Then the next point. It gets a lot of momentum after that third point. And drives uphill really far. But eventually turns around because it gets really steep. So he turns around. Comes back to Sample's point. Goes downhill. And he speeds up as he goes downhill. The width of the line segment here is the kinetic energy. It's his motion. It's momentum, actually, in this physics simulation. And over and over again. At each point, you randomize the momentum again. So that you're still cast dissipated. But because you're following the contour of the shape that you want to map out, you're making good proposals. You're not making silly proposals. But you need a lot more information to make this work. You need to know, well, the contour. So how does this work? Outside of the silly allegory, here's what's really going on.  We really are doing a physics simulation of a frictionless particle on a frictionless surface. The surface is the minus log posterior, which makes a bowl. So look at the bottom. If you've got a Gaussian posterior in one dimension on the left, if you take the logarithm of that and multiply it by minus one, you get a parabolic bowl. And we're going to roll our marble in this parabolic bowl and periodically stop it and record its position. And in the long run, you will get position samples, which are proportional to this shape, as if there was more probability in the bottom as there is on the left. You get more position recordings like that in exactly the right proportions. This is what's cool about it. You're not guessing and checking anymore. Let me show you what this looks like in a movie. So again, two-dimensional Gaussian hill. We're looking down on it from the top. But now we're going to run a Hamiltonian simulation on it. There's going to be a particle. We're going to flick it. This is a bowl now. So you want to see the dark part in the middle as the bottom of a parabolic bowl. Nobody has parabolic bowls at home, I know. But imagine you had a parabolic bowl. We're going to flick our simulation here and do the pass. Again, the green arrow is the transition. It's where you move to. The gray arrow is your initial momentum at the flick point. So we flick it and it rolls in the bowl, ends up somewhere inside the bowl. You're never going to get out of the bowl, the good part of the bowl, because of the physics that's governing it. When you start to go uphill, it makes you turn around. That's what stops you from getting into the silly parts of the parameter space that you don't want to be in. Now, to do this, though, you need to know the local curvature at each point. So the physics simulation has to calculate what's called the gradient, which is the derivative in each direction of the minus log posterior. So you have to do more math. This is a lot more computationally intensive. The computer is doing all the work, though, don't worry. But it's much more computationally intensive than the metropolis algorithm. You need more information. But you accept almost every proposal. You can get 100% acceptance rates routinely in Hamiltonian money problem. And so, yeah, speed it up. It had good times, right? So this is better living through physics. This works really, really well. In practice, since nearly all the proposals are accepted, you need many, many fewer steps to get a good picture of the posterior distribution than Hamiltonian money parlor. So each step, each proposal you consider takes a lot more computer cycles, but you need many, many fewer such steps, such proposals to get a good picture of the posterior distribution. So you get a lot of efficiency gain out of this. Hamiltonian money parlor for a two-dimensional Gaussian, so what? It just seems like a lot of bother. The reason we use it is because life isn't so benign. Instead, we get things like these high posterior correlations. It's actually the same problem I showed you before, where metropolis gets stuck. Now, using the code that I have in the chapter, I show you a simple R-strip to do this. It isn't that much code to run a physics simulation, actually, it's just a loop. And you get these tours. And so what you're seeing in these line segments is we start at the X in the upper left here, and it gets flicked sort of up and to the right, but then it curves back down and rolls down the valley, right? Rocking. Now, there's no friction, so it's never going to stop moving. It just keeps rolling around. And eventually we stop the clock, we record a position, and then we flick it again. So it rolls from the X to number one. You see the number one there, one point. And then from one, it gets flicked again and it rolls around a bit, curls back around to number two, and then we flick it again and it eventually gets over to number three and then to number four. So we're getting, not only do we have a high acceptance rate, but the sequential samples are very far from one another as well. What we call the sequential auto correlation is very low. So you're getting a lot of information about the whole shape very quickly with Hamiltonian Monte Carlo. That's what's great about it. And so then on the right, I'm showing you the same thing, but now for 50 trajectories and 50 samples, there's one open point, which is one rejection. You can sometimes reject things in Hamiltonian Monte Carlo, and we'll talk more about that later on and why that happens and what we can do about it. Okay, another movie. So Hamiltonian Monte Carlo does really well with the donut. That's why we like it. So we're back in the donut. Yeah? Let's run Hamiltonian Monte Carlo on the donut. So it's like it knows the donut's there. Of course it doesn't. All it knows is the curvature where it is. But it still obeys the donut. So now this donut is like a narrow valley.  And it just tours the whole thing. And so you don't get stuck in any particular area. So now imagine you have 27,000 dimensions in your model. This is great. This is really what you want, right? Something like this. So guessing and checking is dead to me. I am forever now going to run physics simulations. So there we go. Full speed. We can tour the donut. Get a really nice picture of it. At this point Hamiltonian Metropolis would still be stuck in a corner, right? So this is really, really important. Sort of innovation in the applied stats community. Okay. Richard, can I ask a side question on this real quick? Mm-hmm. How did you describe Osterian? You're trying to show the whole donut? No, of course you can't. You can't show anything more than three dimensions, actually. I mean, there's this joking in stats. It comes from someone named Jeff Hinton who's worked a lot on these. He says, if you want to picture something in, say, 14 dimensions, draw a three-dimensional shape, and then stare at it intensely and say 14. So nobody can do it. You don't. You can work with high-dimensional models the same way you've been working with the ones in the first part of the course. You look at the prediction scale. You don't interpret the parameters anyway, right?  Have I used that analogy already? It's the tide prediction engine.  But that's what you care about, is the prediction scale where it can make sense. There's some low-dimensional compression at the level of predictions that you can make sense of. And so that's true for these. Yeah, no one can visualize a hypersphere. It's really tough. This is a big problem. In the text, though, there's a short section where I try to explain this phenomenon called concentration of measure. And I've got a graph for you where I would talk about is if you consider the radial distance of your coordinate, from the mode, from the point of highest posterior probability, as the dimensionality of the distribution goes up, most of the samples get further and further in radial distance from the middle. So there's a graph where I show you this. That doesn't mean it's going to be intuitive, though. In my experience, this is all just like, what, radial distance? It's bizarre, but yeah. High-dimensional spaces are very confusing things. But I counsel you to just work with it the way you work with other models and not worry about plotting the whole thing. You can do Paris plots for interesting parts of it, stuff like that. Okay, let me try to summarize. All right, so why does HMC work so good? It doesn't get stuck. It follows the gradient. The gradient is the curvature, the local curvature, at each point. We've got to calculate that. There's extra stuff you have to put into this, including momentum variables now. So this little particle, it really is like, there's a real physics simulation going on here. It has mass. You have to choose that mass, and you have to give it momentum and all these other decisions. And it has energy. The system has energy, which is actually nice because this is how we check if it's working. Energy must be conserved. You may have heard this. Right? So if energy is not conserved, you can check that criterion and you know the simulation is broken. This is something you can do with Hamiltonian Monte Carlo that you can't do with metropolis. You cannot know very easily if metropolis is not working. When Hamiltonian breaks, it's dramatic and obvious. This is really nice. It gives you a sense of security. There'll be all kinds of warnings that come from these extra seemingly nuisance parameters. They're not nuisances at all. They help you feel secure that the thing is working. But you have to put in more. You need gradients, as I said, so we have to calculate that. What does that mean? It means taking derivatives of the log posterior at any given point. You don't need the whole posterior here, but at any particular point you can evaluate it because you have your model definition. But then you need derivatives of that. For any particular model you want to program, like the example in the text I gave you, I just did the derivatives by hand. There's a box where I show you that. Gaussian derivatives are super easy. The Gaussian distribution is always easy. In general, this is annoying, super annoying. This is why the tool we're going to use does this for you automatically. It is an aggressive user of the chain rule. Your computer can do all the derivatives for you. There is tuning that needs to be done, though. I'm going to tell you a little bit about that. We're not going to stay on this slide. In the book, if you're interested, there's a working R script that will run a Hamiltonian simulation. I used it to make those figures. It's not that complicated. If you want to tour through it so you really understand how dead simple this is. There's a particle. We know the curvature at this point so we can simulate its movement for a little while. That's all it's doing. This is not rocket science. The problem we often encounter with Hamiltonian Monte Carlo is they're stuck to pick. The two things that we need to pick to tune it are the step size. The step size determines the length of time we run each little segment of a simulation. If you're looking on the left, this is the two-dimensional Gaussian bowl. We start at the X and we cruise to one. We see those little dots I've plotted on here. Each of those is a step size. It's a time interval that the thing runs. You've got to choose that. This determines the coarseness of the simulation. Those of you who've done much numerical integration, you're used to this, right? You choose the stiffness of your numerical integrator when you do ODE's and things. This is like that. You want basically the biggest step size that you can get away with because it makes it more efficient, but if you make it too big, then you overshoot the curvature of the shape. You've got to choose that. Then you have to choose how many of those steps you're going to take in each so-called trajectory between the points where you stop. Those are the two things you need to choose. If you choose bad values for those, then you're going to have a bad time. In fact, you could have a maximally bad time with this algorithm. Here's an example. This only works maximally bad in the two-dimensional Gaussian. In general, it's not this bad, but something happens called the U-turn phenomenon where you can revisit where you started almost exactly every time if you pick exactly the right step size and number of what are called leapfrog steps. Each of those little steps is called a leapfrog step for reasons I cannot tell you. It's a good name. There's a little frog, I guess, jumping on the surface. The frog is counting one, two, three, four, and it always jumps a certain distance. Since it's a parabolic bowl, you get these parabolic orbits in the bowl.  Then you flick it again and then you're right back where you started. In the long run, this will still work because you're not coming back exactly where you were, but it's super inefficient again because you're getting a bunch of samples from the same local space. This is called the U-turn problem. How do you fix this? You've got to have good values for these tuning parameters, and that's annoying. Let me show you in the animated version what the U-turn problem looks like. Again, here's our Gaussian bowl. Since I'm evil, I've set the step size in a number of leapfrog steps so that it eats its own tail. This is Ouroboros, that snake that eats itself forever. Is that Mayan? Someone know. This is the Ouroboros of statistics, the Markov chain that eats its own tail. It gets stuck in this little region. It is wandering, but it's wandering very slowly. In the long run, this doesn't work. This doesn't ruin the chain. But we could do a lot better if we had better values for the step size. See, it's working in the long run, but it's not working as well as it did before. See, it's stuck in the regions for a little while. Since you don't know the distribution, though, it's hard to say before you start sampling what good step size and number of steps are. So how do you deal with this? Well, we're going to use an engine that deals with this for you. It uses a primitive form of artificial intelligence to find good step sizes, and then it runs the chain. So we're going to use a library called Stan. Stan is a C++ library. You will run it from within R, though. You just have to install Rstan, the package.  Yeah? And my package will interface with Stan for you in returning samples. And Stan does two things. First, it has something called a warm-up phase where it figures out good step size. It figures that out by trying a bunch of random simulations and finding a step size which maximizes the acceptance-raising proposals. Right? Good step sizes. And then it uses a really clever thing to determine the number of leapfrog steps called the no-U-turn algorithm. The no-U-turn sampler or nuts, it's called. And this is super clever, and I'm not going to explain the details of it, but I would be want to show it to you in animation, what it looks like. The key thing about the no-U-turn sampler is that it runs the physics simulation in both directions. What do I mean directions? I mean directions in time. So it imagines a simulation that loops back on itself and it runs that one backwards from your starting point and it goes forward at the same time so it can guess when the turning around is happening. And when it sees itself turning around, it stops and it takes a sample. Let me show you what this looks like. And it sounds like madness and maybe it is. Yes, here it goes. See, it goes in both directions. Yeah. It's a nice thing about Hamiltonian systems is their time invariant so you can run it in both directions. You can go backwards in time. You can figure out where the particle must have been if it's on the trajectory. So this is why it works. And so it guesses when it's turning around and it takes smart positions by doing that in both directions. I know it looks super weird but this means you don't have to pick the number of leapfrog steps. It figures out a good number of leapfrog steps for every individual trajectory that is custom and the location it's in. And this gives you good uncorrelated samples from the distribution. Now it starts to speed up. Yeah. So this means you don't have to pick that.  A smart artificial intelligence will do it for you. Okay. So I wanted to say why is this thing called Stan? Stan is not an acronym. It's the name of a person. Stanislaw Ulaw pictured here was one of the people who, one of the mathematicians working with Nick Metropolis and that whole group that was actually directed by John von Neumann that was part of the Manhattan Project. And Ulaan is really the father of Marco Tain Monte Carlo. He figured out a bunch of the basic things and built mechanical Marco Tain Monte Carlo simulators. You can find these online. These wooden things. You can demonstrate them. They're mechanical tools that will take random samples. And he also did a lot of really important work in biology which he's standing next to that molecule model there. And so Stan is named after Stanislaw Ulaw. And you can run Stan in anything you like. If you go to the Stan webpage you'll see this but we're going to run it in R. If you want to use Python I know some of you prefer Python. Python is an inherently better language than R. I'm sorry I got to say it. You can run Stan in Python too. It's fine. You can also run it just on the command line. When I run models on the servers the computing cluster in my department I just upload a C++ executable to the server and then put it on a million cores. Not a million. 80. And then let it go. So I don't... R is like a lot of overhead that you don't need necessarily, right? You just want samples. Okay. Here's my favorite picture I have to show you. This is Stanislaw Ulaw with his daughter Claire at Maniac to show you this on the same project. That was a computer, right? They couldn't even contain one cat, Jeff, right? I love that phone too. Phones used to look like that. Okay. There are problems with Hamilton and Monte Carlo still. If you've got a posterior distribution like this one that has multiple separated hills and this does happen this happens a lot in latent variable models like factor analytic models then it can be very hard to transition between the hills and your waiting time between hills can be long indeed. So you get stuck on the hill for a while and again in the long run you're guaranteed to explore this space but the long run now can be pretty long. So people are working on this in my experience you handle models like this by changing the geometry of the model. That sounds crazy. It's not cheating. You're not changing the answer but you collapse hills together by dealing with symmetries in the model. But you've got to realize there's a problem first before you do that. So in my department we get problems like this when we run item response theory models which have this kind of labeling exchangeability to them. Okay. Let me give you in the last five minutes here a little bit of advice on what this looks like when you run it but you're going to go home to the chapter and you're going to run all these examples yourself. You've got to get in there and do it and play around with it and walk through. I give you two examples at the end of this chapter of stuff that goes wrong and how to fix it because I want you to recognize a bad chain when you see one and sometimes they're bad. They're usually bad because the model is bad not because the algorithm is bad. So here's the sound of one-hand crapping. It's been a joke I've been waiting all turn to make. So when we run Quap this is a model you ran before for the terrain ruggedness example from Monday and Quap results. Now let's do this model with a Markov chain and I'll show you what it looks like. We're going to use Oolong which is in the rethinking package and it looks exactly the same as using Quap, same formula but we're going to typically I'm going to encourage you to make a little slim data list to pass in that only contains the variables of interest and then you need to decide how many chains to run and how many quarters to put them on and how long to run each chain down at the bottom. So here we're going to do four chains we're going to use put each one on a separate core I don't know how many cores your laptop has but probably four. Most laptops have two true cores and four virtual cores these days. If you don't have that many cores you can still put in four your computer won't explode. Nothing bad will happen. And when you run this what does Oolong do? First it translates this style model formula into raw stand code. I show you there's a whole box in the chapter where I show you what that looks like it doesn't look too different but there's a bunch of formal variable definitions that will be awkward if you haven't done C programming before. Stand code is basically C code and it's like a weird dialect of C++ and then stand builds a custom nut sampler to know you turn sampler and then the sampler runs feeds your data in and it runs and you get samples back and then it's just like samples it's just like you've been working with samples the whole term just samples and you analyze them the same way. So we run this what happens here you can see if you show it it'll report each chain how long warm up and sample took what is warm up warm up is when it's picking the step size it's not getting valid samples but it's doing little experiments to figure out how long the numerical integration link should be figures that out and then the sampling period is typically shorter and then the total so this goes pretty fast yeah this model did it in a second it's so fast so the total number of samples you're going to get back are the length of each chain minus the warm up and by default warm up is half and I encourage you not to change that you don't need more than a couple hundred samples from the posterior if for a Hamiltonian chain if it's efficient remember the central limit theorem you get a really good estimate of the posterior mean very rapidly from the central limit theorem okay and then preci works the same way but you get these extra diagnostic things on the end called in F and our hat I say more about these in the chapter in F is the number of effective samples what I want you to see here is that the number of effective samples is greater than the actual number of samples you took from the Markov chain this is not wrong how can that be true that can be true because Hamiltonian Monte Carlo takes dispersed samples so what does it mean effective number of samples that means it's a number of samples that you would get if there was no auto correlation between sequential samples you can do even better than the true number of samples in Hamiltonian Monte Carlo because your samples are really far apart they're exploring different regions these sequential samples have a completely different region of the space they're negatively auto correlated when your sequential samples are negatively auto correlated you're doing better than best and you can actually have the Markov chain that's better than totally random I know it's crazy but this is how good Stan is it's really amazing version one of Stan was not this good but they've gotten so good at tuning it but now it does this routinely for simple models so I just want to show you this because it's so cool okay so and then our hat is the Gellman Rubin convergence diagnostic it starts out above one and converges to one it's a ratio of variances it's an F statistic and so you want it to converge to one and that tells you that the chains aren't different the variance it's a ratio of the variance between chains and within chains right and you want if the chains are all the same then that should all collapse they should all be exchangeable chains they should all look the same the information in any one chain is the same as the information in the other chains okay so there's you can look at the I give you a little bit more detail in the chapter things like pairs plots which I don't tend to use a lot but you can look at the whole posterior this way just as before and then at the end of the chapter I give you two examples which I really want you to go home and run yourself and your computer and fight with and go through patiently in the text and read through it simple things that happen all the time ways to check your chains to make sure they work there are examples of bad chains I'm here because I don't want to keep you but there are cases where what you want to see is this I call it the hairy caterpillar ocular inspection test you want your chains look like hairy caterpillars that doesn't mean they're good but if they don't look like this they're probably not good it's one of those kind of things what you might get instead again you can walk through this on yourself in the chapter is something like this this is bad news yeah this patient is not going to make it so and I show you how to fix this there are lots of common things you can do to fix it so you can get the hairy caterpillar back alright now I give you two examples the thing I want to drive home where I let you go is that typically when you're having a problem getting your chain to work is because there's something wrong with your model definition this is my experience so Andrew Gelman calls this the folk theorem of statistical computing he learned this the hard way just like all of us who do applied stats you'll fight with the thing you'll try tuning parameters and all kinds of stuff and then you'll look at your model and you'll realize oh I don't have a prior for three of the parameters that's why it's not working so what you want to do first is check your model yeah it's probably not the Markov chain it's probably you that's nearly always alright and the problem is as we say in computer science before the keyboard alright so remember the folk theorem there's one more example with bad chains I have a new homework problem up for you in which you're going to model the Wine's 2012 data so this is a combined do your first Markov chains and work with some simple interactions at the same time I hope you find the data interesting I think they are and I want to remind you if you haven't updated the book PDF or the rethinking package lately you should do so there have been a couple minor bug fixes for chapters yet to come in the rethinking package and like every week I upload a new version of the book where there's stuff we're still going to get to that I'm editing as we go okay so the next week I promise you we'll do some more exciting stuff on Monday I want to show you the value of maximizing entropy it's really nice to have maximum entropy I know it sounds bad but actually it's the state you want to be in and I'm going to use that to teach you about generalized linear models so we can do much more powerful modeling of different kinds of variables okay thank you have a good weekend


Well, first of all, thank you very much, James, for inviting me, and this is my first trip to the annual meeting, and I hope it's not my last. I'm really impressed by it, thank you. The topic of the meeting is such that I have to wonder why, what part of it I fulfill. I think I'm pretty much the uncertainty, actually, as I look at it. But I did notice that James was suggesting that one of the things that we're continually doing is looking at new models, changes in models, and that's what I'd like to present today. Some of my co-authors are here, and they will answer any of the particularly curly questions, and one of the ones who isn't here, I will show some of his work, Sam Roy, who with his wife is expecting a child momentarily. So if we have questions, perhaps I'll cover those. First of all, the motivation, the roadmap, and what I'll be talking about. The Firm, which is the failure earth response model, is a model that we put together in response to concerns that many of the existing process models were not capturing some of what I was particularly interested in, and also it was not producing the kind of coupling between the geodynamics and the theory of geodynamics and the theory of the surface processes. So the two main principles of this are first one, which I promise not to read the slides in general, but these I think are things which will be centered to what we're doing. That is that the material detachment, no matter what the stresses applied are, still occurs as a function of the response to the stresses. So plastic failure or a discontinuous process occurs as the material is stressed and that is independent of whether the stresses happen to be applied by geomorphic processes or by geodynamic or tectonic processes. Secondly that the large displacements affect the strain field, the stress field, that is the strain affects stress and it also affects the strength and the strength field is going to be hugely important and this gives us a predictability in much of what we're going to talk about and a relationship in the topographic memory of the strain field and it's a very important component of what we'll be doing and I'll probably start out with that. So I'm really particularly interested in using the high frequency information that we have, that is the low frequency information, something I'll briefly go through, but we've been dealing with it for quite a few years now and it's really time to move on to the tougher problems I think and this is from the geodynamics, I think the rest of you are dealing with high frequency stuff all the time. So I will start out by looking first at the strength and strain relationship that is mostly number two of this, the failure earth response model firm and then we'll take a look at how that couples to geodynamics and how it's predictable and then from there we'll take a look at the introduction of an alternative formulation for surface dynamics that I hope we will develop in parallel with the existing one. Well first of all, let's clarify the kinds of frequencies that we're talking about, the low frequency behavior in geodynamics is something that we've been thrashing about for a long time and it's the relationships really between long term climates and converging or extensional or strake slip tectonics and they generally describe something like the, oh yeah, so that's not, if that pointer is offensive to anybody let me know, I won't pay any attention to you but you can let me know. This is a cross section across the Himalaya from India up in Tibet, it's a standard one, it's one of the early kinds of cross sections that we modeled and we modeled often and there's very few geodynamics modelers who haven't brought this up but all of this is, I think we're self-agreed to learn with it but it's also things that we can do relatively simply. What I'm more concerned with here and in fact in most of my research and most of what we're doing is looking at groups and back so maybe this isn't such a good idea with a pointer. It's mostly looking at the very high frequency spatial ridge valley relationships and the high frequency spatial relationships are usually associated with high frequency temporal problems too and I'm interested predominantly in this talk dealing not just with a gross Himalaya and the Tibetan plateau but also the small ridges and valleys and the rapid fluctuation of those ridges and valleys from time frames which exist at let's say in the hundreds of years down to the milliseconds of seismic displacements. So I guess one of the reasons is why bother and in fact as I look at this talk it's mostly about why bothering and what we're trying to produce. First of all the modern solutions these are solutions the kinds that I would prefer to work done by programs like child and so on do a fantastic job of representing the long term and sometimes short term fluxes and climate variations to the sediment that's being transferred but one of the obvious problems and one that is I find personally difficult to get over is that we look at the earth in terms of different physics depending upon whether we're looking down on the earth with surface dynamics or up at the earth with geodynamics and I'm trying to bring those two together. In addition the properties that we use in much of the surface dynamics are often difficult or impossible to measure things like erodibility and so on are important parameters intuitively they make sense but they're very difficult to measure but there are standard measurements of material properties that come from the engineering world and those are ones that I would like to move into. And third the kind of advances that we're seeing and many of them I think will be presented here allow us to begin to look at the full 3d stress field and rapidly fluctuating and slowly fluctuating and so the concept that we deal only with principal stresses and so on is something that we can now move away from and move on to something new. Another reason is because there are a number of parameters of a number of processes that occur in natural dynamics which are important processes I'll do it a couple times because I quite like it which are important processes we recognize their important processes such as seismic generated landslides and in addition there's a there fluctuations in poor pressure which we fluctuations in poor pressure fluctuations of poor pressure which also cause sliding and failure and they're very difficult to handle within the present scheme of our surface dynamics problems that doesn't mean that they can't be handled but the it becomes increasingly ad hoc sort of like a British automobile of the 1960s with more and more things hung on to a rather attractive initial model but by the time one gets it well developed it really no longer is the I feel is quite the correct way of handling many of these problems. So the I apologize I really have given a power point before in my life. Okay so let's go back first of all we're going to whip through the low frequency variations and in this case I'm taking a model of the Himalaya so it's about a 2000 kilometer model by a thousand kilometer model by about 100 vertically maybe a bit more than that and these are the ones that I suggest are relatively standard in this case of the in this particular region here we're indenting India or running India into Tibet from the left to the right and these are the representative strain fields above and I don't want to spend too much time on it except to say that much of what we're dealing with in geodynamics is of course is a very large field and we'll be looking at embedded melt models for the kinds of scales that we're interested for the surface dynamics and then feedback to the large models. In terms of the just advection link which is an important link then it's quite a straightforward picture that is where there's advection there's vertical displacement then we get topography constructed so if we start to one way in which we can couple the geodynamics with the with the surface models is by coupling it in with rheological changes as a function of rapid removal of material and one of the examples of this is the tectonic aneurysm the tectonic aneurysm model is developed really as if in response to these rapidly uplifted blocks which as they rapidly uplift they carry their temperature signal signature but also their strength signature they weaken the region there's a nonlinear feedback between the weakening and the erosional scheme can produce very large mountains and three examples which we've discussed are things like negaparbat where the erosional regime is a focused river power namtrabawar at the other side of the Himalayan chain which is also focused river power and another one that we've looked at more recently the St. Elias where they're focusing us through a vicious leisial erosion and that's a useful coupling and it's something that we've talked about for a while but one of the ones i'd like to look at now is the higher frequency one in which strength is coupled to strain that is to displacements and this is not something that is new certainly to those of us who deal with faulting and we recognize that a faulted region produces a gouge zone and weak zone shatter zone associated with it and the question that we're addressing is what does this how does this impact in the three-dimensional way on the evolution of the topography so we'll take an example which is from the channel river channel fault in New Zealand and the picture of the fault is looks like a this is this is actual data here by the way most of what i have won't be so let's take a quick look at it and the picture is a core zone with fault gouge very low that's a very low cohesion on the order of kilopascals and then going into the megapascals schist on either side with a zone a fault zone which i think everybody is familiar with but it's going to be important as we go on which is a decreasing intensity of fractures as one gets further and further from the zone if erotability is defined in most of the models in terms of an inverse function of cohesion it's actually of the one over the square of cohesion then a change in cohesion makes a big difference makes a huge difference in terms of what we're going to see and it also means that the orientation of these planes is hugely important so here's a generic model that's been published off and on for some time now and represents pretty clearly the way we see many of these faults although specific faults differ again a weak core and strong edges so the cohesion structure varies vastly in these systems from something which and they represent a more cool home space is down on the kilopascale level to something which is associated with the intact rock on the order of several orders of magnitude different and if we were to take a look at this we'd see that rivers in general have this kind of strength and the buttresses of negoparabhat in other places have a strength that exceeds that so the difference is huge particularly when we put it as a one over c squared term well what does it do to the to their standard model and this is the same standard model that we had before except in this case we're allowing the material the earth surface material to change its cohesion as a function of strain so this little video that we just played and that I'll play again is a is the result of a deformation field a three-dimensional deformation field which is imposed allowed to grow as a function of the strain localization so it's a strain softening those you could think of as fault zones those fault zones have the cohesion which is which is on the order of orders of magnitude weaker than the intact material and the erosional scheme is dominated the river scheme is dominated by that so the net effect of this is that the the strain field dominates the strength field and it the orientation 3d is something that's represented in the topographic field so at higher orders when these are fully coupled together with a flak model excuse me a deformation model and together with child we get something that in which we can begin to see the relative importance of the various components as it evolves and I'll play it again because it's quite nice one and in this case what we're saying is to sort of pick this apart a little bit more that is so the strain and strike field are coupled consequently this will be something in the memory of the topography which is associated with effectively mantle kinematics the in the first case there's the the elevation field that gets superimposed by developing cohesion field and in this cohesion field we see the strike strong bits the high peaks and the weak bits of the valleys and that then further if we put them both together we begin to see that elevation and cohesion aren't always the same but there is a relationship between the two and it also controls the velocity field so there's a complete feedback between the strain weakening which weakens both the erosional field and the material parameters the finite element solutions for the for the deformation field and clearly it is it's necessary as far as we're concerned to be able to represent this in a three-dimensional way so these patterns are very persistent that is the very strong regions can exist for a very long time the weak ones get eroded very very quickly which gives multiple multiple characteristic times to an origin and those characteristic times give characteristic rates which are very very different and an example of this is given by Sam Roy's earlier work in which he looks specifically at an uplifted block and the propagation of a in this case of a nick point like behavior as a function of weakening and the weakening occurs on these two fault zones one which is dipping at 30 degrees one which is dipping at 90 degrees the weakening has the same kind of shape as the fault zones that we were looking at before with a weak core and a more solid edge and if we run this it's not too surprising the the material outside of this by the way has no weakening it's on the order of intact rock and Albert again what Sam showed pretty clearly is that the rate of propagation is a function of certainly the weakening but also the orientation of the weakening so the fault zone that weakens very quickly and goes straight across as a strike slip fault equivalent whereas the other one is a dipping reverse fault equivalent in terms of the orientation of the strain plan this is important in that the the propagation of the rates is therefore coupled to the the strain the nature of the strain so the nature of the strain produce or the signal of the nature of the strain represents in the rates as well as the orientation Sam has since taken this with us and particularly helped from Greg has taken this I think in somewhat further and began to look at what happens when the dipping fault moves across into a system by which we have a the whole landscape develops as a function of the deforming region and what is quite clear is that the asymmetry associated with the dipping fault and erosion back along the fault plan produces a signal which again is characteristic that is it sits it will sit in the topographic memory and allow us to unpack these unpack from the topography something about the memory in the future so the intermediate conclusions of the strength strain relationships are that first about values a week and that the heterogeneity starts from the very beginning it is as soon as the origin begins to form heterogeneities there and from a modeling perspective almost all of probably all of us but certainly all of the things that I've done in my in my modeling life to assume homogeneity being with and it's a bad assumption how much of what we do it's we we will still continue to do it because it gets us started but it almost always gives us incorrect intuition from the beginning and I think that shows up particularly well the strength the anisotropy it's basic control in all origins and it reduces the complexity of the surface that is the increased complexity of the processes reduces the complexity of the remaining surface it's one of the nice things about systems which could be fractal otherwise and the 3d history is is a permanent memory memory in that topography but there's still a number of problems and opportunities that I'd like to go on about and one is the specifically the physics of the of the mechanisms which actually allow the displacement and bringing together the two sets of theories for surface behavior and geodynamic behavior so in this particular case what we will what I will try to do is look at the 3d strength field and this 3d strength stress field particularly and as examples take a look at the New Zealand Southern Alps this is a photo looking from the west to the east across the alpine fall and we'll talk about this region to be a little bit to begin with one of the problems of talking about new theory is that by necessity we criticize or implicitly criticize old theory which is not really my intention here so I want to make it clear that I'm responsible for some of this old theory so I can at least criticize myself on some of this and that is one of the points that James made he said that some of our results look like they look actually like some of the observations and therefore we think that they're right and in fact this is something which has got me a lot of free beer over the years because it does look very it is reminiscent of an origin in which we produce ridges and valleys and it's produced predominantly by a diffusional mechanism and diffusional advection coupled together it's been very useful it will continue to be useful but I think that we have to move beyond that or I have to move beyond that certainly case that I'm going to take a look at is the plate boundary of New Zealand and the area is in the Mount Cook region we take a look at two things one which is a glaciated valley and one which is a badly slumped valley right on the fault there so the point about this is that one of the one of the standard geomorphic surfaces that's examined are the iconic nick points and nick point behavior and nick points are generally described often described in terms of the relationship between the stresses which are generated from fluvial system and the response of the earth that is there's a threshold if it's overstepped and material can be removed in this particular case we're looking up a valley from the that I the first one I pointed out Franz Josef at Glacier and we're looking at the at a nick point which sits right about at the edge of the of the valley fluvial stresses are represented by that little piece of principal stress space that sits down in the same orientation as the river and this would be a fairly standard way of representing whether the stresses are going to be overcoming on in addition to that though there are also glacial stresses not too far away and the glacial stresses also operate on the same piece of rock but they operate with a different set of of conditions roge age shear stresses and so on in addition at the same place there are also the slope stresses those slope stresses not only act in the direction of the glacier in the direction of the valley but of course they also operate very strongly in the direction of the ridges which is right angles to the direction of the fluvial channel here and those stresses can be very large and they can be large enough to contribute to the tectonic behavior too so there's already a coupling between topography and tectonics through the stresses of the topography but in fact it gets a little more complex than that too so if we take a look at this we can see that that the the stresses at the bottom of the slope are very crudely represented by a reverse shear sense at the top of the slope they're represented by an extensional shear and there's also tensile failure possible further up as well in other words one slope has multiple different stress states it can't easily be represented in terms of a transport or failure in terms of a single value and then we put on tectonic stresses so again that same poor piece of rock is also affected by the tectonic stresses and then dynamic stresses those of earthquakes and that's really where i'm hoping to go today so if we then look at this through time that piece of rock is affected by all of those stresses plus the time dependent or climate dependent at glacial load on top all of these happening at a particular place so in order to construct firm to deal with the fact that these various stresses are responsible for the the failure therefore whether a threshold is overcome or not we reformulate the failure all of those stresses into a single tensor referred to as the total stress tensor and it has all of those components and it's this is mathematically not a difficult thing to do by the way it just means that we put into a single reference frame which isn't a river-based or glacier-based or slope-based but rather a geocentric based we then look at the earth failure in terms of the effective stress formulation allowing fluid pressure solutions this is something we normally do in geodynamics and engineering applications so it's not hard and we apply we allow a time-dependent strength material behavior and then we solve in 3D and in this case we use a constructed meshless method and that meshless method is something which has i think great hope for the future i'd like to talk to other groups about producing something here but this one is not a sophisticated one particularly and doesn't have any any transport in it so let's take a look first of all at the other system which is that great big landslide landslide system and federa has a poster on this that you can ask her questions about it later if we take a look at a model from this then that excuse me the model here is one which is a topographic model which then the stresses are solved for in this particular case the cohesion fields come from measured cohesions in the region and so we're looking up the valley across the alpine falls doesn't really matter too hugely at a three-dimensional model and the representation here is red is relatively strong blue is relatively weak and the strength stress relationship is such that we're going to look at the evolution of that strength stress relationship james is eyeing me so this is a strength stress relationship and its evolution as a function of fluid pressure if one takes a look at it you can see that there's a the weak region that is less than one migrates across the system as a function of changing the partial pressure fluid pressure we all know that that landslides are a function of partial pressure the local fluid pressure local fluid propagation we can solve for that now and in addition this other one is another example which we put on we take the same system this happens to be dried and put the wet one here but we look at the applied tectonic strains at the at the boundary and we look at the evolution of the failure system as a function of the application of those if you take a remember the first animation and this one you'll see that the blue marches at a very different way as that blue marches i don't know if you can see it in this slightly scruffy impact but the topography is changing significantly and it changes differently as a function of the kind of input that we that we would put on well there are ways in which we can i want to skip these because there's something else i want to get to particularly like we can because we now have the surface and the geodynamics in the same same reference frame we can solve continuously for both and therefore we no longer have an ad hoc replacement of one on the other they're actually fully integrated system but what i'm really interested in because i'm getting old and i'm interested in actually moving this thing forward fairly quickly and i'd like to talk to you about is the dynamics and that is the thing that i find most exciting about this is that we have the opportunity to look at some of the older paleo seismology in terms of a rather new technique and use all of the information available in the topography and that information is associated not just with a single scrap that may decay but rather with every piece of that of the topography that we can measure says something about its present history in terms of its stress state and its previous history in terms of acceleration and stress state so one of the things that i'm particularly keen on is and in fact is not particularly difficult is looking at the response spectra of a single of any particular part of the earth this gives two regions one in which is the response spectra of parts of new zealand which are calculated of course they're calculated for california and many other places and that response spectra is a function of the it's a function of the period and not only function of the period but also the nature of the earthquake so strike slips faults at different distances produce a different acceleration spectra than other faults and of course they have a different spatial relationship and the spatial relationship is also a good indicator of the spatial relationship of acceleration and the acceleration response spectra has a fantastic amount of information about the nature of tectonic geomorphic coupling at very high frequencies so there we go and we'll quickly look at one other example that we're doing i'll go through very quickly which is an ice sheet model applied to alaska in which we couple the presence and absence of the of an lgm load with a very high prism data set and we look in this particular case at the position of the high velocity ice streams and we couple it to our three-dimensional geodynamic model and these are some of the implications and i think that it's time to finish for questions thank you this could be quite disconcerting to me what you just said because if we have a block of homogeneous material and then you apply a a stress field to it propagate strain weakening let it interact with the surface processes then the steady state solution can has a lot of information in it that can get us back to what the stress field was and so on but in real earth materials with an inherited history hundreds of millions of years pre-existing faults and strains and so on it seems to me what you're saying is we're screwed because you comment on it. Rudy's asked me to comment on whether we're screwed or not with these and in fact my feeling is that we're less screwed with these than we are otherwise that is that these the weakness zones and so on are controlling variables and i think that we can demonstrate that in a number of ways so by not bringing them in we are we are making a very large mistake but beyond that i think that the the coupling through the firm approach allows us finally to look at the effect of the position of weak zones versus strong zones for instance in a river migration in valleys and ridges if we suggest that the the strength structure is only related to the one in the valley then and we play the same physics and the same extrapolation to the ridges we know the only thing we really know is that we have to have it wrong and so in this case this allows us to solve by a way which is in which we can go out and measure things like cohesion tensile strength and so on it allows us to solve and compare it directly to the to the observations so i appreciate the comment about the complexity which i think is real but i think this is a solution to it rather than the problem thanks you uh you noted that this approach can give us some solutions to things like poor fluid pressure landsliding that have traditionally given us trouble are there things in the surface processes that go the other way that you know you're modeling less than optimally by forcing into this other framework thank you yeah i think that there are a number of things i one of the real strengths of the of the current process directed regimes is that it allows us to look at integrated fluxes and basically it is an integrated flux if you use q over a certain amount of time and this is not particularly good at that but i don't as i look increasingly at this scheme i think that this does not lose any of those of the capabilities that is if we can define the flow field then we can calculate the stress field in addition the more we see about the the really exciting hydrodynamics where they produce a three-dimensional stress field the more we see that we can immediately put it on here so i don't except for the fact that we're early on on this i don't see it that there's any particular drawback other than my own clumsiness at bringing them together