Global sensitivity analysis in stochastic systems
Le Maître, Olivier
Stochastic models are used in many scientific fields, including mechanics, physics, life sciences, queues and social-network studies, chemistry. Stochastic modeling is necessary when deterministic ones cannot capture features of the dynamics, for instance, to represent effects of unresolved small-scale fluctuations, or when systems are subjected to important inherent noise. Often, stochastic models are not completely known and involve some calibrated parameters that should be considered as uncertain. In this case, it is critical to assess the impact of the uncertain model parameters on the stochastic model predictions. This is usually achieved by performing a sensitivity analysis (SA) which characterizes changes in a model output when the uncertain parameters are varied. In the case of a stochastic model, one classically applies the SA to statistical moments of the prediction, estimating, for instance, the derivatives with respect to the uncertain parameters of the output mean and variance. In this presentation, we introduce new approaches of SA in a stochastic system based on variance decomposition methods (ANOVA, Sobol). Compared to previous methods, our SA methods are global, with respect to both the parameters and stochasticity, and decompose the variance into stochastic, parametric and mixed contributions. We consider first the case of uncertain Stochastic Differential Equations (SDE), that is systems with external noisy forcing and uncertain parameters. A polynomial chaos (PC) analysis with stochastic expansion coefficients is proposed to approximate the SDE solution. We first use a Galerkin formalism to determine the expansion coefficients, leading to a hierarchy of SDEs. Under the mild assumption that the noise and uncertain parameters are independent, the Galerkin formalism naturally separates parametric uncertainty and stochastic forcing dependencies, enabling an orthogonal decomposition of the variance, and consequently identify contributions arising from the uncertainty in parameters, the stochastic forcing, and a coupled term. Non-intrusive approaches are subsequently considered for application to more complex systems hardly amenable to Galerkin projection. We also discuss parallel implementations and application to derived quantity of interest, in particular, a novel sampling strategy for non-smooth quantities of interest but smooth SDE solution. Numerical examples are provided to illustrate the output of the SA and the computational complexity of the method. Second, we consider the case of stochastic simulators governed by a set of reaction channels with stochastic dynamics. Reformulating the system dynamics in terms of independent standardized Poisson processes permits the identification of individual realizations of each reaction channel dynamic and a quantitative characterization of the inherent stochasticity sources. By judiciously exploiting the inherent stochasticity of the system, we can then compute the global sensitivities associated with individual reaction channels, as well as the importance of channel interactions. This approach is subsequently extended to account for the effects of uncertain parameters and we propose dedicated algorithms to perform the Sobols decomposition of the variance into contributions from an arbitrary subset of uncertain parameters and stochastic reaction channels. The algorithms are illustrated in simplified systems, including the birth-death, Schlgl, and Michaelis-Menten models. The sensitivity analysis output is also contrasted with a local derivative-based sensitivity analysis method.
 Thank you. I think it's perfect time for an app, so I'm very grateful to see you here. I'd like to thank the organizer for the invitation. I think it's a very great event and I was very pleased to be here this week. So my talk will discuss how to perform some global sensitivity analysis in a stochastic system and perhaps to motivate a little bit this question and the methods I'm going to introduce. I must say that for many years we have been working on uncertainty quantification but essentially parametric uncertainty. So models having some parameters that are not well known and we developed techniques to deal with the situation, so parametric uncertainty. But in many situations the model itself is not deterministic. So it's not that we have parameters that we don't really know which value to give but the dynamic, the evolution of the system in time is truly stochastic. So we have some inherent stochasticity. So many physical systems have these situations. It can be because at very small scales usually you have noise, typically thermal noise that can appear and so the evolution is really inherently random. So here I'm thinking of molecular dynamic chemical systems for very small size reactors for instance. Reactors could be a biological sense. There are also situations where there are some forcing acting on the system we want to simulate and there's no way we can predict what will be this forcing. So this could be for instance a case in finance, it could be the case in many situations in energy, think for instance of wind turbine, it is submitted to turbulent winds, storms or whatever and you cannot predict a priori what condition the machine will have to stand on support. So in engineering stochastic modeling becomes very important in that case. And there are obviously the situation of all models that in order to maintain tractability, computational tractability, we have to reduce, we have to solve at say coarse scales. While we know that at smaller scales there are some phenomena that we are completely neglecting. Sometimes we can average or buy some homogenization techniques upscale these effects but there are situations where we don't know how to perform such upscaling and one possibility in this case is to add some random noise that mimics the effect of the small scale process that we cannot model. That's the case for instance in climate modeling where they add some random noise to force the flow at small scale to mimic the effect of turbulence mixing. So stochastic modeling and here I'm just thinking about physical models but you can think of all models related to social science for instance. So stochastic modeling is very common. I can be very complex. It's a very nice way to deal with things that we cannot prescribe or that we don't know. But the cost of doing so is that usually you have to prescribe some parameters typically. Say I don't know what will be the wind so I decide to model it with stochastic process but now I have to prescribe some parameters say time correlations, intensity of the wind gusts or the six line. So we have now a system that has some inherent stochasticity because this is a modeling choice that we made and some parameter and this parameter usually cannot be measured. They need to be learned or to be inferred from information. If we have some experimental data for instance but they are not certain. So now we have two sources of uncertainty. One that is inherent to the model and the second one that comes from lack of knowledge about the parameter. And this is really the question we wanted to address with my colleague Omar Knew to see how we could propagate and perform some sensitivity analysis in such systems. So I will have two parts in this presentation. In the first part I will focus on the case of stochastic ordinary differential equations. So typically scalar ones, so equations or ordinary differential equations driven by some vener noise. And I will show how we can consider this equation plus some parametric uncertainty on how to solve this equation by means of polynomial chaos expansion. And the second part I will go to the stochastic simulators and treating this type of system and we will have to introduce new questions that we believe are quite interesting. So I am considering a very simple SDE. So X is a state. Here everything is scalar but you can extend everything to a higher dimension. So vector case for instance. So you have the evolution of your solution X. You have a drift term with a C that may depend on X. It could also depend on T. And you have diffusion terms here with the vener process. So very classical. So what you can do is for instance to simulate this equation to generate samples or trajectories of this system by considering for instance the other schemes. The simplest one you can devise and what you have to do is to draw some IID normal random variables with variance depends on the time step. So you can integrate a trajectory and repeat multiple times going independently realization of this increment. And when you have a set of trajectories you can build some statistic, mean standard deviation, correlation, whatever is needed. And if you consider a... Oh, I thought it was a spider. So if you have a functional of your trajectory, it could be just the value at a certain time or it could be an exit time or whatever it is. So I call it G and now you can estimate the expectation of this quantity of interest X. The same thing is done by simply by sample average, generating M replicas of a trajectory of the system. So the situation I want to tackle now is exactly the same except that now the drift and the diffusion coefficient may have some non-parameter. That could represent special variability. That could be just the magnitude of... if it's just a random variable, it's just the magnitude of the noise that is entering the system. So I did not queue these random parameters and I'm assuming that they have a known probability. So how do you know this probability? Well, perhaps they are coming from some Bayesian inference based on some experimental observation and that would be your posterior of your Bayesian. And the main assumption I'm doing here is that the knowledge of the parameter doesn't say anything about the forcing process as it is. Two are independent, which is quite, I think, natural and I actually don't see situations where this would not be the case. So it's not a strong assumption, however, it is crucial for what I'm going to develop later. Okay, I'm just here making some assumption about square-integrability. So the process has to be a second-order process. If I condition on the random walks, so I still have the random parameter here, so this conditional variance is finite and that is the same if I condition on a certain value of the parameter, the variance is also finite. This equality here are 0.1. Okay, and our goal is now really to say I want to investigate how the variability in X depends on the different source of the variabilities, which are the linear process and the parameter and I want to separate this different effect and being able to draw some conclusion. So the typical approach that people are following consists in conditioning on a certain value of the parameter, then you average with respect to the process and you do that for some moment of X. So typically you're trying to look to the dependence of the mean of the solution of your process as a function of the parameter. So you average with the process and you keep the parameter out and you can do that for a second moment. So you have now the dependence of the variance with the parameter. And now you have functions here that only depends on the parameter and you can apply the classical sensitivity analysis or I mean the parametric one that I've been developing and I will give you a flavor of the idea that are behind. But here is the idea and one possibility is to go through the functional approximation of this two first moments, but you could go to higher moments on a polynomial basis for instance. Okay, and once you have this representation, you can sample as much as you want this approximation and study some statistics that can characterize the dependencies. Okay, now that's the classical approach. So in my view, you are losing some information because the first step you do is always averaging with respect to the venerable process. So there is a loss of information and somehow you're not too penalized because you can do these analyzes for higher and higher moments. But how many moments you need and it becomes very difficult to interpret the outcome of the sensitivity analysis. So in that among different parameters, you can say one is very important for the moment of order five. What does that mean? It's not clear. So what we want to do is slightly different. It's an orthogonal decomposition of the variance. So X has a certain variance and I want to split these variance into different contributions, contributions that is due to the inner noise, so the venerable process, contributions that is due entirely on the parameter and the contribution to the variance of the interaction. And the way we do that is by the soboil of the decomposition of X where the process that I see here as a function of the vener and the parameter is decomposed as a mean value, a function that depends only on the vener, a function that depends only on the parameter and a mixture. And we ask this composition to be, this function is a decomposition to be mutually orthogonal and this makes the decomposition unique. In the case where W and Q are independent, that's where my assumptions become crucial. And here you have the definition of this function. So for instance, the function that depends only on the vener is the conditional expectation of X given W minus C, full expectation. So because this function is orthogonal, you have immediately the decomposition of the variance. You can normalize the partial variance and you get what we call the sensitivity indices. So they are all between 0 and 1. Their sum here is equal to 1. This one characterizes the impact of the vener process on the variance. This one of the parameter and here the mixed contribution. I will illustrate this idea of characterization on an example later. So here is the interpretation of this coefficient and what you can show is that the variance of the mean with respect to W of X for given value of the parameter, so this is really the variance with respect to the parameter of the mean of the solution. This is actually the variance due to W only. Sorry, to Q only. So clear interpretation of what we have classically. However, if we average with respect to the parameter, the dependence of the variance to the parameter, we see here that we have the sum of two contributions. So the contribution due to W to the vener process and also the mixture. And actually there is no way that given this term here, you can split into these two contributions. That's why we believe that the approach we are proposing is richer. We have a finer analysis in particular. If you discover that SW is really large and dominant compared to the other term, there is very few chance that you can learn your parameter by observation. This is kind of important conclusion that you can do. So how to compute this coefficient? What we propose to do is then because we see the process as a function of W on Q, it is to us quite natural to have a series representation. So I'm considering basis of a complete orthogonal set of function psi alpha, which are functional in Q. Actually, I introduce here the measure of the parameter PQ. So there are an orthogonal family. And we are expanding the solution X in this form. So the coefficient in the expansion, the X alpha are again stochastic processes. But they don't depend on the parameter. So I have this expansion. Now I can plan this expansion in my original equation. I will get a residual. If I truncate my expansion, I get a residual. So I'm requiring this residual to be orthogonal to the representation basis. So to every of the psi alpha I have in my expansion. So I'm forming here the garb of projection of my problem. And I end up with a system for the mode or the coefficient of the expansion. And you will notice that this equation, so if you gather all the coefficient X beta in a vector, you end up with an equation that exactly the same structure as previously, except you were a scalar equation, it's a vector equation. But otherwise you can still recognize here a drift term and here a diffusion coefficient term. And what's quite special is that W is still unique. So if you started with W scalar, you still have the same stochastic process that acts on all the modes, force all the modes. Okay. So I'm not going to detail how this garb of machine is implemented or how it is solved. But basically instead of having a single ODE solved by a pair of schemes, now you have a system of ODE, a coupled system. It's not very complicated. Now based on the expansion coefficient that I have in the series here, if I take for, say, convention that the first of this function is a constant. So X0 will be the mean of the solution, well, mean with respect to the parameter. And now I can have the expression of the different function in the Sobolevding expansion and the expression of the different partial variance, which basically if you want to compute, say, the partial variance due to W, only you have to compute the variance of X0. Okay. So basically for each trajectory, you average with respect to the parameter and you look to the variance of the result. But this is made quite easily thanks to the polynomial chaos expansion. Well, maybe I didn't mention, but typical basis we use here are polynomial chaos. So it's polynomial functions. Okay. So let's take the first example. So linear additive system. And I have two uncertain parameter. One is Q1. It appears in the drift term. So it's basically the location of the attracting point. And for the diffusion, I have another parameter Q2 that appears here. If I make nu equals zero, I have an additive noise problem. And if nu is different than zero, we have a multiplicative noise problem. Okay. So I will start with nu equals zero. And I will consider also that Q1 and Q2 are independent. So in the case where we have an additive model. So nu equals zero. You can show that you need just a degree one polynomial expansion to exactly capture the solution. Okay. And here are trajectories of the mean mode, X0. So you, it really looks like typical trajectories for the system or you process. Here is trajectories for the mode that is linear in Q1. So due to the drift term, and you see that there is absolutely no variability. Because basically the system has a completely deterministic evolution. If you fix Q1, plus some noise that comes and had, but they don't mix. Okay. So if you fix Q1, so here you have a range and we respect to W. And well, this is the average we suspect to W. That's all you have. And this is the mixed term, sorry, the term related to the magnitude of the forcing. So this one obviously has an impact on the solution. Okay. If we have a multiplicative noise term, now it's a bit different because the two parametric sources can interact. So now you have variability on the linear term in Q1. And obviously as before, some stochasticity in other mode related to the linear term in Q2. You can even go to higher order, but very quickly, due to the structure of the system, you see that there is no variability in this mode. So you could even, if you were smart enough to analyze the structure of the galerking system, you could drastically reduce the size of the basis. There's absolutely no need to simulate these trajectories. There is no variability at all. Okay. So these are just densities of the coefficients in my expansion to show that some of them are really Gaussian, as you would expect for Gaussian forcing, but some of them are not Gaussians because this is due to the interaction between the modes. Here are also some protection of coordinates in planes to show how there exists some strong dependence between the coefficients. So how the galerking system is really coupling the dynamic between the different modes. Here's perhaps the simplest plot you could imagine to illustrate the idea. What we have here is a single realization of the Brunyan. We have fixed the noise and we are varying the parameter. Okay. So you see here's different trajectories for different values of the parameter, but you observe that they are very similar. I mean, they are shifted. They can be enlarged, but they are quite self-similar. Okay. And this is really this similarity that allows to have a polynomial expansion that converges very quickly. Here's the opposite. We have fixed the parameter to a certain value and we are running different values. And you see that the stochasticity is clearly very much different. Okay. But now we can decompose. So here is a set of trajectories and you can decompose into the function that depends only on the parameter. So essentially we have removed all the noise and we just have variability with the parameter. We can do the opposite. That is, we look or we average over the parameter and we have here variability that is due to WLE. And here's a mixture. Just for illustration. Okay. This one is not that much important. And here are the, not the sensitivity coefficient I have introduced before, but the partial variance directly. So it's not normalized because the variance varies a long time. So it will be difficult to report everything in terms of normalize. So if you have absolutely no variability or no uncertainty on the parameter, so they have variance that are both zero, then all the variance of the system is due to the noise. Okay. That's only W as a source of variability. Now if you have variability also in the drift term, but not in the magnitude of the diffusion, the total variance in red here is essentially due to the variance due to the noise and you have a contribution of the parameter, but you have no interaction. That's what we have seen before where the old trajectory was exactly the same. So this is what is happening here. On the contrary, if you have here a determinist, well, a perfectly known drift term and some uncertainty on the magnitude of the noise, what you have now is a contribution of the, of the winner, of course, but also of the mixture. Okay. And now if you have uncertainty in all the parameters, so drifts and diffusion, now you have all the contributions that are non-zero. So with this way, we can really distinguish what is the inner run variability, what is due to the parameter and have a very fine understanding of what's going on. Okay. Here's another more complicated system. It's, I have just changed the drift function. So the system has two attracting branch. Okay. And we are starting from a point somewhere in the middle and depending on the realization of the Brunyan, we are led to go to a first attracting branch or to another one. Of course, there is a non-zero, but very low probability that when you have selected a branch, you could also branch to the other one, but this is a very rare event. Okay. And for this system, I'm considering some uncertainty which are on the initial location. So clearly if I'm starting from the very high point, it's very likely that I'll end up on this branch and not on the second one. So I have uncertainty on the starting point and I have also uncertainty on the magnitude of the stochastic force. Okay. So here is just an illustration of what is happening depending on my starting point. So if I start from the left, I always go on the same branch or essentially always go to the same branch. If I stop somewhere on the unstable point, I may select with equal probability one branch or the other. And as a contrary, if I start on the other side, I tend to always end up on the same branch. Okay. Here is the effect of the noise magnitude. So I'm starting always at the same point. If I have a very low stochastic forcing, basically I always go into the same branch. I mean it's like deterministic attractor while I stay on the deterministic attractor. And if I increase the noise level, I am able to change this picture. Okay. But of course we would like to do this in a more quantitative way. Okay. Although we understand perfectly what's happening in the system, we want to provide numbers. So as you've seen, even for a given starting point, depending on the brilliant or for a given brilliant, depending on my starting point or intensity of the noise, I can either finish on a branch or another, meaning that for a large time, my solution can be discontinuous. I have a branching process. I have a bifurcation. So the polynomial basis I was using before are not going to work because or I would need an extremely high polynomial degree to capture the discontinuity. So we have here to use some multi-weblet basis with some activity in the parameter space. So this is just to illustrate here. You have the mesh of the parametric space. And we have here, I think, initial condition and noise magnitude on the other direction. And what is represented here is the size of the element over which we have a polynomial approximation. Okay, and here is the result. So you have the evolution of the total variance. What you can conclude is that initially it's essentially the variance that is due to the parameter that dominates. Okay, essentially it depends where you start before the noise can have time to act. So as time increases, you have more and more contribution of the noise, but it's not very much important. Okay, and later on there is the mixed contribution that emerged, which in fact represents a case of events where without noise you would have uncertainty. Only the blue one with the noise is just fluctuation around the branch, so to say. And here the fact that the noise can make you switch from branch to another, which also depends on the parameter. Okay, this is exactly the same, but for higher noise level, and you see that now the magnitude of the different contributions tend to be more balanced as before, because you have more noise in the system. It's quite what you would expect. Okay, as a final note, here it was made with a gallerking projection. This can be difficult if you have a complicated nonlinear drift function, perform the actual projection of the system, but you can also proceed in a nonintrusive way by selecting some value in the parameter space, computing the trajectories for these specific values, and then construct by say regression for instance, or quadrature, the coefficient of the expansion. And you can do that, and perhaps this is a possibility to better illustrate the interest of going through this polynomial chaos approximation. You see that in this case, the coefficients here, or to advance, or to compute the coefficient for a given realization of the vener, you would have to solve NQ trajectories for NQ different values of the parameter. The complexity here is from one trajectory, you go to NQ trajectory, but if your approximation here is sufficiently converged, then you can query your approximation for any value of the parameter, and so you can resample as much as you want. And if most of the variance is carried by the parameter, as opposed to the vener, this can be extremely effective. But there remains a lot of work to really analyze and come up with algorithms that somehow balance the parametric error, sampling error, and do something very effectively. Okay, one final comment here. Here I'm working directly on the representation of the analysis of the whole trajectory, and there is a good reason. For instance, you may want to apply this technique for the case of quantities of interest like exit time, but exit time are not smooth functions, not even with respect to the parameter. You could imagine that for a given realization of the noise trajectories that goes and exit, you change the value of the parameter, you may be missing the exit boundary and continue and exist much later. So there is really an issue with the regularity. So working directly on the expansion of this function, it will be extremely difficult because you cannot construct such approximation for G directly. However, you can use this trick to approximate the trajectory and resample, and by resampling this trajectory, you can very easily generate many, many samples of the exit time. Second part I want to discuss is the case of stochastic simulators. So Raoul discussed this, and it was a work with actually Alvaro Morais, who was a former PhD student of Raoul. So it's a different class of system, systems that are governed by a master equation. If your background is chemistry, that's how you know it. So basically what it says, it specifies the advancement in time of the probability of your system to be in some state x as a certain time t, knowing that it was in a state x0 at time t0. And on the right hand side, you have basically two contributions, the contributions that makes you leaving the state x at time t, and the contributions that makes you going to that state. So you're receiving some probability mass from elsewhere, and you are dumping probability mass by going elsewhere. So this is really a balanced equation. So in our case, x will be on the lattice. It will represent typically a number of molecules in a certain system, a molecule of different species. And this balance equation here has a sum over k reaction channel. So there are k ways for the system to evolve. So in the chemical system, you have k reaction. The function Hj here are the propensity function that basically fix the time rate at which a certain reaction j occurs. And they are generally dependent on the state. And the mu g here are the state change vectors. So each time you have a reaction of time j that occurs, you have to update the state of your system, and that's the quantity by which the state is changing for a single event. So it's a Markov process, and it's used in many, many different applications. But we have to solve this type of equation, while there are plenty of techniques to solve this equation, but what we consider is the Gillespie algorithm, which is an exact approach for this system. What you do is that you really keep in time, or keep track in time of the evolution of the system by jumping from a reaction to the next. So if you have k reaction with current propensity function Hj, you know that in the next time interval t plus t plus dt, this is the probability of the next reaction to work. And it's exponentially distributed. So the more reaction you have, the more likely it is that it occurred in a small period of time. So the Gillespie algorithm is quite easy to understand, so you start from your initial state, x0, then you draw the time to the next reaction, given that you are in a certain state, you can evaluate the rate of reaction here. You decide of the time to the next reaction, so you will advance your system to this time. And you pick at random the reaction that fires first with relative probabilities that are given by the propensity function. So the ones that have large propensity function are more likely to have occurred first than the ones that have low probability function. And then you update your system, and you continue till you reach your final time. So it's really mimicking all reaction in the system. So our goal here is to say, what if I have some uncertainty or some parameter in this propensity function, how can I deal and perform the same type of analysis that I had before using these type of simulators? But, well, actually we had absolutely no clue on how to do that. So the first question we wanted to address is, even for this simple problem, so I can generate multiple trajectories, I can estimate the variance of the population size at some time t, but I'd like to be able to answer and say, if my question is among this k reaction, which one are responsible for most of the variability? So I have no parametric variability, but already I'd like to address this question, which are the reactions that are causing the most of the variance. And surprisingly, we didn't find any answer to this question, which seems quite important, but no. What people are doing is generally they average, they look to the average trajectory and eventually they look to the sensitivity of this average trajectory with respect to the coefficient here, but it's really direct. So how can we do that? So here's a typical example of a set of trajectories. You have here two reactions, how to distinguish the effect of the two is really unclear. So let me just generalize a little bit more the symbol of the decomposition. So I have a functional that depends on n independent random functional, and here I'm just writing the orthogonal decomposition that I had before, but it is generalized. So I have the random input in my functional, here's just an ensemble notation for the symbol of the decomposition. And here it's the definition of the orthogonality, and here's a variance where you sum about all possible, well, all functions in your decomposition, except for the first one, which is a mean. So how can you define this decomposition? Well, here is the definition of this partial variance. So if you consider a set of the direction specified by the vector u, then the variability due to this vector is simply the variance of the expectation of f given this stochastic input minus the sum of all the variance of the functions that contain this set of indices. And you have here the other definition. And from that you can define two important indices, what I call the first order indices, that is the part of the variance that is due to only a certain subset of input, and only this subset. And then I have the total variance, oops, sorry for that, here, so for subset of variable or input in u, u is the index of the input that I want to characterize. So the total variance is not only the variance that is due to this input, but also this input and their interaction with any other input. And I think that in this situation where I have independent input, I have a nice way, a Monte Carlo approach that was proposed by Subaul, to estimate this partial variance and then this sensitivity indices, so first order and second order. The computational effort you have to do is basically compute correlations between your stochastic output when you keep part of your inputs exactly the same, so you generate exactly the same replica for part of the inputs and you redraw the inputs that you want to analyze. So basically for each coefficient you want to compute, you have an additional trajectory to solve, so you double, but the cost is eventually scaling linearly with the number of indices you want to compute. So if you just want to separate or if you have k-channel in your system, you will have to increase as compared to say estimate classically by Monte Carlo moment or an expectation, say now you have to estimate end time or as many input more trajectory, but it's still linear. Okay, but how to apply this to the stochastic simulator? Here in this Gillespie algorithm it's very unclear what are the stochastic inputs on the channel, but more exactly when I draw a random number, I can hardly say okay, this is a number drawn for certain reaction channel. Actually, if I change the indexation of the channel, I may completely change the dynamic. So how can I decide of things or inputs that are staying the same for some channel and others that are varying for the channel? And this was really smart idea by Alvaro Morais was to say okay, what if you go to the Kurt's representation where you represent here the evolution of your system in terms of independent standard Poisson processes because now each channel in this representation has its own Poisson processes. So you can repeat them or you can resample them and you can decide to keep the same sum of them and we compute others. And then you will be able to compute some correlations. And of course, it doesn't mean that by conditioning on a certain realization of Poisson process for a certain channel, you don't change the reaction time of this channel because everything is linked. Okay, so here is typically you have to compute some conditional expectations. So you may fix the value of some of the Poisson processes associated to a subset of channels and you can resample the other one. And then you can compute the correlation to get the sensitivity in this. So if you already have this representation, Kurt's representation, you already have all you need to answer the question of what channels are important in terms of variables. In fact, this is the key position in some of Poisson processes, variable processes. It's done under which assumptions on the coefficients Ae? I don't think there is any particular... But this is... We can make it for any... Any type of systems. Any type of systems? Yeah. Okay, without any assumption on the coefficients Ae? I mean... No, no, no. Those on the intensities of the germs, in fact. No. Now in this problem, you don't have any restriction on our assumption. Okay. It's also known as the time change representation. Okay, it's Tg, so it's a scale... It's Tg here, so the Poisson here are normal, but the argument here depends on the... I see, it's linked to the integral of Ae. Okay, thank you. Okay, so let's see what the conclusion of such analysis is. So here I'm considering the simplest system, so the death birth process. So I have individuals that are with birth rate B. It's independent of the state because it's spontaneous generation. It's not couples that are making children. No, it's just... So you have a birth that occurs once in a while at a rate given by B, and then people are also dying, but at a rate D that is linear in the size of the population, so it's D times X. The more individual you have, the more frequent is to observe a death, or people dying. Right? Okay, so if you apply the Gillespie or the Kurtz algorithm, you can generate trajectories. There is an invariant measure for T going to infinity, which is actually a Poisson distribution that you have on the right. Okay? And now what we would like to say in this variance that we may have, say, at T equal 8, is it B that is more important? Is it my coefficient D, or is it the birth or the death channel? Actually, I should say. Okay. Do you have ideas? I had no idea. Okay, so here are the results. So at very early time, all the variability is due to the birth process, because we start with a zero population, so if there is no individual, no one is going to die, and there will be no variability on the number of dead people, right? Okay. So initially, it's all the birth process that carries the, sorry, here, the variability, and as time is increasing, variability due to the death channel is building up. The sickness of this line here represents the distance between the first order sensitivity indices and the total order sensitivity indices. So the sickness of the line represents the interaction between the two channels. Okay? And actually, if you run the computation over a very long period of time, okay, here it's a log scale, you see that eventually the contribution of the death channel is going to die, okay? Because if you don't have variability in the birth, well, nobody remains. Okay? So over a long time, you need to have variability in the death to have some, in the birth to have variability in the death, okay? And eventually, everything will be due to the birth channel and interaction with death. Okay? Now, a more complicated system, this is a Schlogel system. It's just four reactions, well, actually two reverse reactions. We have three species, B1, S, B2, but B1 and B2 are considered in large excess and are considered constant. So we just modeled the evolution of the number of individual species S that we call X. Here as the propensity function, so you have two body reactions, three body reactions, et cetera. And here are some trajectories. For this system that represents the infection of a cell by a virus, you have basically two asymptotic states. One is the virus has replicated a lot, okay? Or the virus is not replicating and remain at a very low level. So infection is low or infection is high. And you have here the density, not the histogram, more exactly at T equal 8 of the solution. Okay? So for this system, when we carry the analysis we had before, so we want to characterize the effect of the four channels that we have, this is what we obtain as a variance is essentially due to the first and second reaction, while the second and third are much less important. We see that we have strong interactions. But the main conclusion is that first and fourth reaction are the dominant one in terms of inducing variability. And how can this be explained? It can be explained by the fact that it's essentially these two channels and their variability that plays on the selection of containing a low level of viruses or giving birth to lots of replicas. So selecting the upper branch or the lower branch is essentially dictated by channel one and four. While the two other channels induces some noise, but they are way more, way less likely to induce a complete shift of branch. That's what you have here in this plot. So here it's on each of these small plots you have a draw of N1 and N4, so the process of channel one and channel four, and we repeat for other realization of channel two and three. And you see that most of the time you stay on the same branch, it's quite unlikely that you shift. While on the other case, which is completely the opposite, you see that you tend to be less selective in your branch. Okay, that's another system, but I have a taxi to take, and I don't think you would be very much interested into the details of these results, but it's just too stress that you can have a quite fine understanding of how the different channels are combining their stochasticity. It's quite obvious that now we are, and perhaps I even have this, you can add some parameters in the propensity function and you consider them as additional source, independent source of uncertainty. So this is a case for instance where you add in the propensity function an uncertain parameter. This is a situation where you draw different question processes. This is a picture you had before, but it's always the same parameter. While here, it's the same question processes, but you change the value of the parameter. Okay, and again, the trajectories are more similar. Okay, and now you are able to separate and say, okay, my sensitivity indices have a contribution due to the parameters, that's a blue one here, due to the channel, stochasticity, and one that is due to the mixture. And obviously in that case, it would be quite easy to, I think, to learn the parameter because our uncertainty on the parameter is a dominant source of variability, so we should be able to learn quite easily these parameters. The opposite, if we had here a yellow contribution that would be far up, why I would say that you don't really care about the actual value of the parameter because it doesn't give you much variance, and it will be very difficult to learn them. Okay, I think maybe I would prefer to answer questions rather than commenting these slides. Any questions? Like when you are discussing the sensitivity analysis for the stochastic ordinary differential equation, you have like capital Q1 and capital Q2, and you assume that they're independent and you've done the sensitivity analysis in that way. If you just have an idea, if you don't have this independence between this Q1 and Q2, how it's changing, it makes it easier, the sensitivity analysis. Or to consider them independent. If I consider them dependent. Dependent? Yes. How it changes these sensitivity analysis? Well, first you could still continue to perform a separation between inner and stochasticity and parametric one, but it would have to be globally the global contribution of the parameters. Now, if you want to have a finer analysis and distinguish between the contribution of Q1 and the contribution of Q2, then if they are independent, you can apply the symbol of being the composition. If they are not independent, you cannot do that anymore. So there are some, recently they were proposed some techniques to perform this kind of variance decomposition in the case of dependent variable, and perhaps you could reuse them. But the interpretation of the sensitivity coefficient that you get at the end or the characterization that you get is not very well understood at the moment. Typically, you have to account between correlation between parameters, and you can have sensitivity coefficients that can even become negative, which is really strange, right? So thank you Olivier for the presentation. So to come back to the sensitivity analysis, so here you have presented something related to path 5 sensitivity, where you conditioned one motion and so on, but you could also perform a sensitivity analysis regarding the PDE underlying the stochastic differential equation, so that is we know by the distribution of linear PDE. So you could perform also the PCE on these PDE. So could you comment a bit about the choice of doing one way compared to the other? So what is your motivation? There are some, well first I wanted to work by Swise, which is a good motivation. One issue with trying to work directly on the densities or on the master equation is rationality. First, I mean, systems that we want eventually to solve are high dimensional, so solving this master equation in high dimensions is very complicated. That's one point. Second point is that the approximation techniques we are using, so polynomial chaos, honestly they are not that good at enforcing positivity. So expanding density can be problematic. That's another motivation. And basically for us it was making more sense to say, okay, these simulators are existing, stochastic simulators or your favorite tools. Let's try to reuse this machinery as much as possible. And as you've seen, this is actually what we are doing. Not with the gallerking, but it's, gallerking projection is my data, but you could do it non-intrusively and then you just need to have a solver, which I think is a very strong point for users. I'm not very much aware of people working on the direct solution of the master equation. But honestly, positivity is a big issue. Okay, so I think it's time for you to leave. Thank you. Thank you very much.
 Thank you.  I think it's perfect time for an app, so I'm very grateful to see you here.  I'd like to thank the organizer for the invitation.  I think it's a very great event and I was very pleased to be here this week.  So my talk will discuss how to perform some global sensitivity analysis in a stochastic  system and perhaps to motivate a little bit this question and the methods I'm going to  introduce.  I must say that for many years we have been working on uncertainty quantification but  essentially parametric uncertainty.  So models having some parameters that are not well known and we developed techniques  to deal with the situation, so parametric uncertainty.  But in many situations the model itself is not deterministic.  So it's not that we have parameters that we don't really know which value to give but  the dynamic, the evolution of the system in time is truly stochastic.  So we have some inherent stochasticity.  So many physical systems have these situations.  It can be because at very small scales usually you have noise, typically thermal noise that  can appear and so the evolution is really inherently random.  So here I'm thinking of molecular dynamic chemical systems for very small size reactors  for instance.  Reactors could be a biological sense.  There are also situations where there are some forcing acting on the system we want  to simulate and there's no way we can predict what will be this forcing.  So this could be for instance a case in finance, it could be the case in many situations in  energy, think for instance of wind turbine, it is submitted to turbulent winds, storms  or whatever and you cannot predict a priori what condition the machine will have to stand  on support.  So in engineering stochastic modeling becomes very important in that case.  And there are obviously the situation of all models that in order to maintain tractability,  computational tractability, we have to reduce, we have to solve at say coarse scales.  While we know that at smaller scales there are some phenomena that we are completely  neglecting.  Sometimes we can average or buy some homogenization techniques upscale these effects but there  are situations where we don't know how to perform such upscaling and one possibility  in this case is to add some random noise that mimics the effect of the small scale process  that we cannot model.  That's the case for instance in climate modeling where they add some random noise to force  the flow at small scale to mimic the effect of turbulence mixing.  So stochastic modeling and here I'm just thinking about physical models but you can  think of all models related to social science for instance.  So stochastic modeling is very common.  I can be very complex.  It's a very nice way to deal with things that we cannot prescribe or that we don't know.  But the cost of doing so is that usually you have to prescribe some parameters typically.  Say I don't know what will be the wind so I decide to model it with stochastic process  but now I have to prescribe some parameters say time correlations, intensity of the wind  gusts or the six line.  So we have now a system that has some inherent stochasticity because this is a modeling choice  that we made and some parameter and this parameter usually cannot be measured.  They need to be learned or to be inferred from information.  If we have some experimental data for instance but they are not certain.  So now we have two sources of uncertainty.  One that is inherent to the model and the second one that comes from lack of knowledge  about the parameter.  And this is really the question we wanted to address with my colleague Omar Knew to see  how we could propagate and perform some sensitivity analysis in such systems.  So I will have two parts in this presentation.  In the first part I will focus on the case of stochastic ordinary differential equations.  So typically scalar ones, so equations or ordinary differential equations driven by some  vener noise.  And I will show how we can consider this equation plus some parametric uncertainty on how to  solve this equation by means of polynomial chaos expansion.  And the second part I will go to the stochastic simulators and treating this type of system  and we will have to introduce new questions that we believe are quite interesting.  So I am considering a very simple SDE.  So X is a state.  Here everything is scalar but you can extend everything to a higher dimension.  So vector case for instance.  So you have the evolution of your solution X.  You have a drift term with a C that may depend on X.  It could also depend on T. And you have diffusion terms here with the vener process.  So very classical.  So what you can do is for instance to simulate this equation to generate samples or trajectories  of this system by considering for instance the other schemes.  The simplest one you can devise and what you have to do is to draw some IID normal random  variables with variance depends on the time step.  So you can integrate a trajectory and repeat multiple times going independently realization  of this increment.  And when you have a set of trajectories you can build some statistic, mean standard deviation,  correlation, whatever is needed.  And if you consider a...  Oh, I thought it was a spider.  So if you have a functional of your trajectory, it could be just the value at a certain time  or it could be an exit time or whatever it is.  So I call it G and now you can estimate the expectation of this quantity of interest X.  The same thing is done by simply by sample average, generating M replicas of a trajectory  of the system.  So the situation I want to tackle now is exactly the same except that now the drift and the  diffusion coefficient may have some non-parameter.  That could represent special variability.  That could be just the magnitude of... if it's just a random variable, it's just the magnitude  of the noise that is entering the system.  So I did not queue these random parameters and I'm assuming that they have a known probability.  So how do you know this probability?  Well, perhaps they are coming from some Bayesian inference based on some experimental observation  and that would be your posterior of your Bayesian.  And the main assumption I'm doing here is that the knowledge of the parameter doesn't  say anything about the forcing process as it is.  Two are independent, which is quite, I think, natural and I actually don't see situations  where this would not be the case.  So it's not a strong assumption, however, it is crucial for what I'm going to develop  later.  Okay, I'm just here making some assumption about square-integrability.  So the process has to be a second-order process.  If I condition on the random walks, so I still have the random parameter here, so this conditional  variance is finite and that is the same if I condition on a certain value of the parameter,  the variance is also finite.  This equality here are 0.1.  Okay, and our goal is now really to say I want to investigate how the variability in  X depends on the different source of the variabilities, which are the linear process and the parameter  and I want to separate this different effect and being able to draw some conclusion.  So the typical approach that people are following consists in conditioning on a certain value  of the parameter, then you average with respect to the process and you do that for some moment  of X.  So typically you're trying to look to the dependence of the mean of the solution of  your process as a function of the parameter.  So you average with the process and you keep the parameter out and you can do that for  a second moment.  So you have now the dependence of the variance with the parameter.  And now you have functions here that only depends on the parameter and you can apply  the classical sensitivity analysis or I mean the parametric one that I've been developing  and I will give you a flavor of the idea that are behind.  But here is the idea and one possibility is to go through the functional approximation  of this two first moments, but you could go to higher moments on a polynomial basis  for instance.  Okay, and once you have this representation, you can sample as much as you want this approximation  and study some statistics that can characterize the dependencies.  Okay, now that's the classical approach.  So in my view, you are losing some information because the first step you do is always averaging  with respect to the venerable process.  So there is a loss of information and somehow you're not too penalized because you can do  these analyzes for higher and higher moments.  But how many moments you need and it becomes very difficult to interpret the outcome of  the sensitivity analysis.  So in that among different parameters, you can say one is very important for the moment  of order five.  What does that mean?  It's not clear.  So what we want to do is slightly different.  It's an orthogonal decomposition of the variance.  So X has a certain variance and I want to split these variance into different contributions,  contributions that is due to the inner noise, so the venerable process, contributions that  is due entirely on the parameter and the contribution to the variance of the interaction.  And the way we do that is by the soboil of the decomposition of X where the process that  I see here as a function of the vener and the parameter is decomposed as a mean value,  a function that depends only on the vener, a function that depends only on the parameter  and a mixture.  And we ask this composition to be, this function is a decomposition to be mutually orthogonal  and this makes the decomposition unique.  In the case where W and Q are independent, that's where my assumptions become crucial.  And here you have the definition of this function.  So for instance, the function that depends only on the vener is the conditional expectation  of X given W minus C, full expectation.  So because this function is orthogonal, you have immediately the decomposition of the  variance.  You can normalize the partial variance and you get what we call the sensitivity indices.  So they are all between 0 and 1.  Their sum here is equal to 1.  This one characterizes the impact of the vener process on the variance.  This one of the parameter and here the mixed contribution.  I will illustrate this idea of characterization on an example later.  So here is the interpretation of this coefficient and what you can show is that the variance  of the mean with respect to W of X for given value of the parameter, so this is really  the variance with respect to the parameter of the mean of the solution.  This is actually the variance due to W only.  Sorry, to Q only.  So clear interpretation of what we have classically.  However, if we average with respect to the parameter, the dependence of the variance  to the parameter, we see here that we have the sum of two contributions.  So the contribution due to W to the vener process and also the mixture.  And actually there is no way that given this term here, you can split into these two contributions.  That's why we believe that the approach we are proposing is richer.  We have a finer analysis in particular.  If you discover that SW is really large and dominant compared to the other term, there  is very few chance that you can learn your parameter by observation.  This is kind of important conclusion that you can do.  So how to compute this coefficient?  What we propose to do is then because we see the process as a function of W on Q, it is  to us quite natural to have a series representation.  So I'm considering basis of a complete orthogonal set of function psi alpha, which are functional  in Q.  Actually, I introduce here the measure of the parameter PQ.  So there are an orthogonal family.  And we are expanding the solution X in this form.  So the coefficient in the expansion, the X alpha are again stochastic processes.  But they don't depend on the parameter.  So I have this expansion.  Now I can plan this expansion in my original equation.  I will get a residual.  If I truncate my expansion, I get a residual.  So I'm requiring this residual to be orthogonal to the representation basis.  So to every of the psi alpha I have in my expansion.  So I'm forming here the garb of projection of my problem.  And I end up with a system for the mode or the coefficient of the expansion.  And you will notice that this equation, so if you gather all the coefficient X beta in  a vector, you end up with an equation that exactly the same structure as previously,  except you were a scalar equation, it's a vector equation.  But otherwise you can still recognize here a drift term and here a diffusion coefficient  term.  And what's quite special is that W is still unique.  So if you started with W scalar, you still have the same stochastic process that acts  on all the modes, force all the modes.  Okay.  So I'm not going to detail how this garb of machine is implemented or how it is solved.  But basically instead of having a single ODE solved by a pair of schemes, now you have  a system of ODE, a coupled system.  It's not very complicated.  Now based on the expansion coefficient that I have in the series here, if I take for,  say, convention that the first of this function is a constant.  So X0 will be the mean of the solution, well, mean with respect to the parameter.  And now I can have the expression of the different function in the Sobolevding expansion  and the expression of the different partial variance, which basically if you want to compute,  say, the partial variance due to W, only you have to compute the variance of X0.  Okay.  So basically for each trajectory, you average with respect to the parameter and you look  to the variance of the result.  But this is made quite easily thanks to the polynomial chaos expansion.  Well, maybe I didn't mention, but typical basis we use here are polynomial chaos.  So it's polynomial functions.  Okay.  So let's take the first example.  So linear additive system.  And I have two uncertain parameter.  One is Q1.  It appears in the drift term.  So it's basically the location of the attracting point.  And for the diffusion, I have another parameter Q2 that appears here.  If I make nu equals zero, I have an additive noise problem.  And if nu is different than zero, we have a multiplicative noise problem.  Okay.  So I will start with nu equals zero.  And I will consider also that Q1 and Q2 are independent.  So in the case where we have an additive model.  So nu equals zero.  You can show that you need just a degree one polynomial expansion to exactly capture the  solution.  Okay.  And here are trajectories of the mean mode, X0.  So you, it really looks like typical trajectories for the system or you process.  Here is trajectories for the mode that is linear in Q1.  So due to the drift term, and you see that there is absolutely no variability.  Because basically the system has a completely deterministic evolution.  If you fix Q1, plus some noise that comes and had, but they don't mix.  Okay.  So if you fix Q1, so here you have a range and we respect to W. And well, this is the  average we suspect to W.  That's all you have.  And this is the mixed term, sorry, the term related to the magnitude of the forcing.  So this one obviously has an impact on the solution.  Okay.  If we have a multiplicative noise term, now it's a bit different because the two parametric  sources can interact.  So now you have variability on the linear term in Q1.  And obviously as before, some stochasticity in other mode related to the linear term in  Q2.  You can even go to higher order, but very quickly, due to the structure of the system,  you see that there is no variability in this mode.  So you could even, if you were smart enough to analyze the structure of the galerking  system, you could drastically reduce the size of the basis.  There's absolutely no need to simulate these trajectories.  There is no variability at all.  Okay.  So these are just densities of the coefficients in my expansion to show that some of them  are really Gaussian, as you would expect for Gaussian forcing, but some of them are not  Gaussians because this is due to the interaction between the modes.  Here are also some protection of coordinates in planes to show how there exists some strong  dependence between the coefficients.  So how the galerking system is really coupling the dynamic between the different modes.  Here's perhaps the simplest plot you could imagine to illustrate the idea.  What we have here is a single realization of the Brunyan.  We have fixed the noise and we are varying the parameter.  Okay.  So you see here's different trajectories for different values of the parameter, but you  observe that they are very similar.  I mean, they are shifted.  They can be enlarged, but they are quite self-similar.  Okay.  And this is really this similarity that allows to have a polynomial expansion that converges  very quickly.  Here's the opposite.  We have fixed the parameter to a certain value and we are running different values.  And you see that the stochasticity is clearly very much different.  Okay.  But now we can decompose.  So here is a set of trajectories and you can decompose into the function that depends  only on the parameter.  So essentially we have removed all the noise and we just have variability with the parameter.  We can do the opposite.  That is, we look or we average over the parameter and we have here variability that is due to  WLE.  And here's a mixture.  Just for illustration.  Okay.  This one is not that much important.  And here are the, not the sensitivity coefficient I have introduced before, but the partial  variance directly.  So it's not normalized because the variance varies a long time.  So it will be difficult to report everything in terms of normalize.  So if you have absolutely no variability or no uncertainty on the parameter, so they  have variance that are both zero, then all the variance of the system is due to the noise.  Okay.  That's only W as a source of variability.  Now if you have variability also in the drift term, but not in the magnitude of the diffusion,  the total variance in red here is essentially due to the variance due to the noise and you  have a contribution of the parameter, but you have no interaction.  That's what we have seen before where the old trajectory was exactly the same.  So this is what is happening here.  On the contrary, if you have here a determinist, well, a perfectly known drift term and some  uncertainty on the magnitude of the noise, what you have now is a contribution of the , of the winner, of course, but also of the mixture.  Okay.  And now if you have uncertainty in all the parameters, so drifts and diffusion, now you  have all the contributions that are non-zero.  So with this way, we can really distinguish what is the inner run variability, what is  due to the parameter and have a very fine understanding of what's going on.  Okay.  Here's another more complicated system.  It's, I have just changed the drift function.  So the system has two attracting branch.  Okay.  And we are starting from a point somewhere in the middle and depending on the realization  of the Brunyan, we are led to go to a first attracting branch or to another one.  Of course, there is a non-zero, but very low probability that when you have selected a  branch, you could also branch to the other one, but this is a very rare event.  Okay.  And for this system, I'm considering some uncertainty which are on the initial location.  So clearly if I'm starting from the very high point, it's very likely that I'll end up  on this branch and not on the second one.  So I have uncertainty on the starting point and I have also uncertainty on the magnitude  of the stochastic force.  Okay.  So here is just an illustration of what is happening depending on my starting point.  So if I start from the left, I always go on the same branch or essentially always go  to the same branch.  If I stop somewhere on the unstable point, I may select with equal probability one branch  or the other.  And as a contrary, if I start on the other side, I tend to always end up on the same  branch.  Okay.  Here is the effect of the noise magnitude.  So I'm starting always at the same point.  If I have a very low stochastic forcing, basically I always go into the same branch.  I mean it's like deterministic attractor while I stay on the deterministic attractor.  And if I increase the noise level, I am able to change this picture.  Okay.  But of course we would like to do this in a more quantitative way.  Okay.  Although we understand perfectly what's happening in the system, we want to provide numbers.  So as you've seen, even for a given starting point, depending on the brilliant or for a  given brilliant, depending on my starting point or intensity of the noise, I can either  finish on a branch or another, meaning that for a large time, my solution can be discontinuous.  I have a branching process.  I have a bifurcation.  So the polynomial basis I was using before are not going to work because or I would need  an extremely high polynomial degree to capture the discontinuity.  So we have here to use some multi-weblet basis with some activity in the parameter space.  So this is just to illustrate here.  You have the mesh of the parametric space.  And we have here, I think, initial condition and noise magnitude on the other direction.  And what is represented here is the size of the element over which we have a polynomial  approximation.  Okay, and here is the result.  So you have the evolution of the total variance.  What you can conclude is that initially it's essentially the variance that is due to the  parameter that dominates.  Okay, essentially it depends where you start before the noise can have time to act.  So as time increases, you have more and more contribution of the noise, but it's not very  much important.  Okay, and later on there is the mixed contribution that emerged, which in fact represents a case  of events where without noise you would have uncertainty.  Only the blue one with the noise is just fluctuation around the branch, so to say.  And here the fact that the noise can make you switch from branch to another, which also  depends on the parameter.  Okay, this is exactly the same, but for higher noise level, and you see that now the magnitude  of the different contributions tend to be more balanced as before, because you have more  noise in the system.  It's quite what you would expect.  Okay, as a final note, here it was made with a gallerking projection.  This can be difficult if you have a complicated nonlinear drift function, perform the actual  projection of the system, but you can also proceed in a nonintrusive way by selecting  some value in the parameter space, computing the trajectories for these specific values,  and then construct by say regression for instance, or quadrature, the coefficient of the expansion.  And you can do that, and perhaps this is a possibility to better illustrate the interest  of going through this polynomial chaos approximation.  You see that in this case, the coefficients here, or to advance, or to compute the coefficient  for a given realization of the vener, you would have to solve NQ trajectories for NQ different  values of the parameter.  The complexity here is from one trajectory, you go to NQ trajectory, but if your approximation  here is sufficiently converged, then you can query your approximation for any value  of the parameter, and so you can resample as much as you want.  And if most of the variance is carried by the parameter, as opposed to the vener, this  can be extremely effective.  But there remains a lot of work to really analyze and come up with algorithms that somehow  balance the parametric error, sampling error, and do something very effectively.  Okay, one final comment here.  Here I'm working directly on the representation of the analysis of the whole trajectory, and  there is a good reason.  For instance, you may want to apply this technique for the case of quantities of interest like  exit time, but exit time are not smooth functions, not even with respect to the parameter.  You could imagine that for a given realization of the noise trajectories that goes and exit,  you change the value of the parameter, you may be missing the exit boundary and continue  and exist much later.  So there is really an issue with the regularity.  So working directly on the expansion of this function, it will be extremely difficult because  you cannot construct such approximation for G directly.  However, you can use this trick to approximate the trajectory and resample, and by resampling  this trajectory, you can very easily generate many, many samples of the exit time.  Second part I want to discuss is the case of stochastic simulators.  So Raoul discussed this, and it was a work with actually Alvaro Morais, who was a former  PhD student of Raoul.  So it's a different class of system, systems that are governed by a master equation.  If your background is chemistry, that's how you know it.  So basically what it says, it specifies the advancement in time of the probability of  your system to be in some state x as a certain time t, knowing that it was in a state x0  at time t0.  And on the right hand side, you have basically two contributions, the contributions that  makes you leaving the state x at time t, and the contributions that makes you going to  that state.  So you're receiving some probability mass from elsewhere, and you are dumping probability  mass by going elsewhere.  So this is really a balanced equation.  So in our case, x will be on the lattice.  It will represent typically a number of molecules in a certain system, a molecule of different  species.  And this balance equation here has a sum over k reaction channel.  So there are k ways for the system to evolve.  So in the chemical system, you have k reaction.  The function Hj here are the propensity function that basically fix the time rate at which  a certain reaction j occurs.  And they are generally dependent on the state.  And the mu g here are the state change vectors.  So each time you have a reaction of time j that occurs, you have to update the state  of your system, and that's the quantity by which the state is changing for a single event.  So it's a Markov process, and it's used in many, many different applications.  But we have to solve this type of equation, while there are plenty of techniques to solve  this equation, but what we consider is the Gillespie algorithm, which is an exact approach  for this system.  What you do is that you really keep in time, or keep track in time of the evolution of  the system by jumping from a reaction to the next.  So if you have k reaction with current propensity function Hj, you know that in the next time  interval t plus t plus dt, this is the probability of the next reaction to work.  And it's exponentially distributed.  So the more reaction you have, the more likely it is that it occurred in a small period of  time.  So the Gillespie algorithm is quite easy to understand, so you start from your initial  state, x0, then you draw the time to the next reaction, given that you are in a certain  state, you can evaluate the rate of reaction here.  You decide of the time to the next reaction, so you will advance your system to this time.  And you pick at random the reaction that fires first with relative probabilities that are  given by the propensity function.  So the ones that have large propensity function are more likely to have occurred first than  the ones that have low probability function.  And then you update your system, and you continue till you reach your final time.  So it's really mimicking all reaction in the system.  So our goal here is to say, what if I have some uncertainty or some parameter in this  propensity function, how can I deal and perform the same type of analysis that I had before  using these type of simulators?  But, well, actually we had absolutely no clue on how to do that.  So the first question we wanted to address is, even for this simple problem, so I can  generate multiple trajectories, I can estimate the variance of the population size at some  time t, but I'd like to be able to answer and say, if my question is among this k reaction,  which one are responsible for most of the variability?  So I have no parametric variability, but already I'd like to address this question, which are  the reactions that are causing the most of the variance.  And surprisingly, we didn't find any answer to this question, which seems quite important,  but no.  What people are doing is generally they average, they look to the average trajectory and eventually  they look to the sensitivity of this average trajectory with respect to the coefficient  here, but it's really direct.  So how can we do that?  So here's a typical example of a set of trajectories.  You have here two reactions, how to distinguish the effect of the two is really unclear.  So let me just generalize a little bit more the symbol of the decomposition.  So I have a functional that depends on n independent random functional, and here I'm just writing  the orthogonal decomposition that I had before, but it is generalized.  So I have the random input in my functional, here's just an ensemble notation for the symbol  of the decomposition.  And here it's the definition of the orthogonality, and here's a variance where you sum about  all possible, well, all functions in your decomposition, except for the first one, which  is a mean.  So how can you define this decomposition?  Well, here is the definition of this partial variance.  So if you consider a set of the direction specified by the vector u, then the variability  due to this vector is simply the variance of the expectation of f given this stochastic  input minus the sum of all the variance of the functions that contain this set of indices.  And you have here the other definition.  And from that you can define two important indices, what I call the first order indices,  that is the part of the variance that is due to only a certain subset of input, and only  this subset.  And then I have the total variance, oops, sorry for that, here, so for subset of variable  or input in u, u is the index of the input that I want to characterize.  So the total variance is not only the variance that is due to this input, but also this input  and their interaction with any other input.  And I think that in this situation where I have independent input, I have a nice way,  a Monte Carlo approach that was proposed by Subaul, to estimate this partial variance  and then this sensitivity indices, so first order and second order.  The computational effort you have to do is basically compute correlations between your  stochastic output when you keep part of your inputs exactly the same, so you generate exactly  the same replica for part of the inputs and you redraw the inputs that you want to analyze.  So basically for each coefficient you want to compute, you have an additional trajectory  to solve, so you double, but the cost is eventually scaling linearly with the number of indices  you want to compute.  So if you just want to separate or if you have k-channel in your system, you will have  to increase as compared to say estimate classically by Monte Carlo moment or an expectation, say  now you have to estimate end time or as many input more trajectory, but it's still linear.  Okay, but how to apply this to the stochastic simulator?  Here in this Gillespie algorithm it's very unclear what are the stochastic inputs on  the channel, but more exactly when I draw a random number, I can hardly say okay, this  is a number drawn for certain reaction channel.  Actually, if I change the indexation of the channel, I may completely change the dynamic.  So how can I decide of things or inputs that are staying the same for some channel and  others that are varying for the channel?  And this was really smart idea by Alvaro Morais was to say okay, what if you go to the Kurt's  representation where you represent here the evolution of your system in terms of independent  standard Poisson processes because now each channel in this representation has its own  Poisson processes.  So you can repeat them or you can resample them and you can decide to keep the same sum  of them and we compute others.  And then you will be able to compute some correlations.  And of course, it doesn't mean that by conditioning on a certain realization of Poisson process  for a certain channel, you don't change the reaction time of this channel because everything  is linked.  Okay, so here is typically you have to compute some conditional expectations.  So you may fix the value of some of the Poisson processes associated to a subset of channels  and you can resample the other one.  And then you can compute the correlation to get the sensitivity in this.  So if you already have this representation, Kurt's representation, you already have all  you need to answer the question of what channels are important in terms of variables.  In fact, this is the key position in some of Poisson processes, variable processes.  It's done under which assumptions on the coefficients Ae?  I don't think there is any particular...  But this is...  We can make it for any...  Any type of systems.  Any type of systems?  Yeah.  Okay, without any assumption on the coefficients Ae?  I mean...  No, no, no.  Those on the intensities of the germs, in fact.  No.  Now in this problem, you don't have any restriction on our assumption.  Okay.  It's also known as the time change representation.  Okay, it's Tg, so it's a scale...  It's Tg here, so the Poisson here are normal, but the argument here depends on the...  I see, it's linked to the integral of Ae.  Okay, thank you.  Okay, so let's see what the conclusion of such analysis is.  So here I'm considering the simplest system, so the death birth process.  So I have individuals that are with birth rate B. It's independent of the state because  it's spontaneous generation.  It's not couples that are making children.  No, it's just...  So you have a birth that occurs once in a while at a rate given by B, and then people  are also dying, but at a rate D that is linear in the size of the population, so it's D times  X.  The more individual you have, the more frequent is to observe a death, or people dying.  Right?  Okay, so if you apply the Gillespie or the Kurtz algorithm, you can generate trajectories.  There is an invariant measure for T going to infinity, which is actually a Poisson distribution  that you have on the right.  Okay?  And now what we would like to say in this variance that we may have, say, at T equal  8, is it B that is more important?  Is it my coefficient D, or is it the birth or the death channel?  Actually, I should say.  Okay.  Do you have ideas?  I had no idea.  Okay, so here are the results.  So at very early time, all the variability is due to the birth process, because we start  with a zero population, so if there is no individual, no one is going to die, and there  will be no variability on the number of dead people, right?  Okay.  So initially, it's all the birth process that carries the, sorry, here, the variability,  and as time is increasing, variability due to the death channel is building up.  The sickness of this line here represents the distance between the first order sensitivity  indices and the total order sensitivity indices.  So the sickness of the line represents the interaction between the two channels.  Okay?  And actually, if you run the computation over a very long period of time, okay, here it's  a log scale, you see that eventually the contribution of the death channel is going  to die, okay?  Because if you don't have variability in the birth, well, nobody remains.  Okay?  So over a long time, you need to have variability in the death to have some, in the birth to  have variability in the death, okay?  And eventually, everything will be due to the birth channel and interaction with death.  Okay?  Now, a more complicated system, this is a Schlogel system.  It's just four reactions, well, actually two reverse reactions.  We have three species, B1, S, B2, but B1 and B2 are considered in large excess and are  considered constant.  So we just modeled the evolution of the number of individual species S that we call X.  Here as the propensity function, so you have two body reactions, three body reactions,  et cetera.  And here are some trajectories.  For this system that represents the infection of a cell by a virus, you have basically two  asymptotic states.  One is the virus has replicated a lot, okay?  Or the virus is not replicating and remain at a very low level.  So infection is low or infection is high.  And you have here the density, not the histogram, more exactly at T equal 8 of the solution.  Okay?  So for this system, when we carry the analysis we had before, so we want to characterize the  effect of the four channels that we have, this is what we obtain as a variance is essentially  due to the first and second reaction, while the second and third are much less important.  We see that we have strong interactions.  But the main conclusion is that first and fourth reaction are the dominant one in terms  of inducing variability.  And how can this be explained?  It can be explained by the fact that it's essentially these two channels and their variability  that plays on the selection of containing a low level of viruses or giving birth to  lots of replicas.  So selecting the upper branch or the lower branch is essentially dictated by channel  one and four.  While the two other channels induces some noise, but they are way more, way less likely  to induce a complete shift of branch.  That's what you have here in this plot.  So here it's on each of these small plots you have a draw of N1 and N4, so the process  of channel one and channel four, and we repeat for other realization of channel two and three.  And you see that most of the time you stay on the same branch, it's quite unlikely that  you shift.  While on the other case, which is completely the opposite, you see that you tend to be  less selective in your branch.  Okay, that's another system, but I have a taxi to take, and I don't think you would be  very much interested into the details of these results, but it's just too stress that you  can have a quite fine understanding of how the different channels are combining their  stochasticity.  It's quite obvious that now we are, and perhaps I even have this, you can add some parameters  in the propensity function and you consider them as additional source, independent source  of uncertainty.  So this is a case for instance where you add in the propensity function an uncertain parameter.  This is a situation where you draw different question processes.  This is a picture you had before, but it's always the same parameter.  While here, it's the same question processes, but you change the value of the parameter.  Okay, and again, the trajectories are more similar.  Okay, and now you are able to separate and say, okay, my sensitivity indices have a contribution  due to the parameters, that's a blue one here, due to the channel, stochasticity, and  one that is due to the mixture.  And obviously in that case, it would be quite easy to, I think, to learn the parameter because  our uncertainty on the parameter is a dominant source of variability, so we should be able  to learn quite easily these parameters.  The opposite, if we had here a yellow contribution that would be far up, why I would say that  you don't really care about the actual value of the parameter because it doesn't give you  much variance, and it will be very difficult to learn them.  Okay, I think maybe I would prefer to answer questions rather than commenting these slides.  Any questions?  Like when you are discussing the sensitivity analysis for the stochastic ordinary differential  equation, you have like capital Q1 and capital Q2, and you assume that they're independent  and you've done the sensitivity analysis in that way.  If you just have an idea, if you don't have this independence between this Q1 and Q2,  how it's changing, it makes it easier, the sensitivity analysis.  Or to consider them independent.  If I consider them dependent.  Dependent?  Yes.  How it changes these sensitivity analysis?  Well, first you could still continue to perform a separation between inner and stochasticity  and parametric one, but it would have to be globally the global contribution of the parameters.  Now, if you want to have a finer analysis and distinguish between the contribution of  Q1 and the contribution of Q2, then if they are independent, you can apply the symbol  of being the composition.  If they are not independent, you cannot do that anymore.  So there are some, recently they were proposed some techniques to perform this kind of variance  decomposition in the case of dependent variable, and perhaps you could reuse them.  But the interpretation of the sensitivity coefficient that you get at the end or the  characterization that you get is not very well understood at the moment.  Typically, you have to account between correlation between parameters, and you can have sensitivity  coefficients that can even become negative, which is really strange, right?  So thank you Olivier for the presentation.  So to come back to the sensitivity analysis, so here you have presented something related  to path 5 sensitivity, where you conditioned one motion and so on, but you could also perform  a sensitivity analysis regarding the PDE underlying the stochastic differential equation, so that  is we know by the distribution of linear PDE.  So you could perform also the PCE on these PDE.  So could you comment a bit about the choice of doing one way compared to the other?  So what is your motivation?  There are some, well first I wanted to work by Swise, which is a good motivation.  One issue with trying to work directly on the densities or on the master equation is  rationality.  First, I mean, systems that we want eventually to solve are high dimensional, so solving  this master equation in high dimensions is very complicated.  That's one point.  Second point is that the approximation techniques we are using, so polynomial chaos, honestly  they are not that good at enforcing positivity.  So expanding density can be problematic.  That's another motivation.  And basically for us it was making more sense to say, okay, these simulators are existing,  stochastic simulators or your favorite tools.  Let's try to reuse this machinery as much as possible.  And as you've seen, this is actually what we are doing.  Not with the gallerking, but it's, gallerking projection is my data, but you could do it  non-intrusively and then you just need to have a solver, which I think is a very strong  point for users.  I'm not very much aware of people working on the direct solution of the master equation.  But honestly, positivity is a big issue.  Okay, so I think it's time for you to leave.  Thank you.  Thank you very much.
2017-07-21_LeMaitre_0_47.png 2017-07-21_LeMaitre_48_95.png 2017-07-21_LeMaitre_96_155.png 2017-07-21_LeMaitre_156_203.png 2017-07-21_LeMaitre_204_275.png 2017-07-21_LeMaitre_276_395.png 2017-07-21_LeMaitre_396_467.png 2017-07-21_LeMaitre_468_1691.png 2017-07-21_LeMaitre_1692_7403.png 2017-07-21_LeMaitre_7404_8423.png 2017-07-21_LeMaitre_8424_8519.png 2017-07-21_LeMaitre_8520_11147.png 2017-07-21_LeMaitre_11148_14207.png 2017-07-21_LeMaitre_14208_17123.png 2017-07-21_LeMaitre_17124_17951.png 2017-07-21_LeMaitre_17952_19727.png 2017-07-21_LeMaitre_19728_22991.png 2017-07-21_LeMaitre_22992_24623.png 2017-07-21_LeMaitre_24624_25139.png 2017-07-21_LeMaitre_25140_25187.png 2017-07-21_LeMaitre_25188_25247.png 2017-07-21_LeMaitre_25248_25511.png 2017-07-21_LeMaitre_25512_25631.png 2017-07-21_LeMaitre_25632_26543.png 2017-07-21_LeMaitre_26544_29135.png 2017-07-21_LeMaitre_29136_29519.png 2017-07-21_LeMaitre_29520_30227.png 2017-07-21_LeMaitre_30228_30695.png 2017-07-21_LeMaitre_30696_31055.png 2017-07-21_LeMaitre_31056_32279.png 2017-07-21_LeMaitre_32280_32915.png 2017-07-21_LeMaitre_32916_33623.png 2017-07-21_LeMaitre_33624_36323.png 2017-07-21_LeMaitre_36324_36419.png 2017-07-21_LeMaitre_36420_36875.png 2017-07-21_LeMaitre_36876_38075.png 2017-07-21_LeMaitre_38076_41159.png 2017-07-21_LeMaitre_41160_41543.png 2017-07-21_LeMaitre_41544_41723.png 2017-07-21_LeMaitre_41724_41855.png 2017-07-21_LeMaitre_41856_42083.png 2017-07-21_LeMaitre_42084_42227.png 2017-07-21_LeMaitre_42228_44051.png 2017-07-21_LeMaitre_44052_44231.png 2017-07-21_LeMaitre_44232_44279.png 2017-07-21_LeMaitre_44280_45983.png 2017-07-21_LeMaitre_45984_46487.png 2017-07-21_LeMaitre_46488_49607.png 2017-07-21_LeMaitre_49608_50075.png 2017-07-21_LeMaitre_50076_54227.png 2017-07-21_LeMaitre_54228_56975.png 2017-07-21_LeMaitre_56976_57587.png 2017-07-21_LeMaitre_57588_59795.png 2017-07-21_LeMaitre_59796_59963.png 2017-07-21_LeMaitre_59964_60023.png 2017-07-21_LeMaitre_60024_60863.png 2017-07-21_LeMaitre_60864_60911.png 2017-07-21_LeMaitre_60912_61727.png 2017-07-21_LeMaitre_61728_62243.png 2017-07-21_LeMaitre_62244_63803.png 2017-07-21_LeMaitre_63804_64787.png 2017-07-21_LeMaitre_64788_64835.png 2017-07-21_LeMaitre_64836_65075.png 2017-07-21_LeMaitre_65076_67163.png 2017-07-21_LeMaitre_67164_67439.png 2017-07-21_LeMaitre_67440_69443.png 2017-07-21_LeMaitre_69444_71267.png 2017-07-21_LeMaitre_71268_74051.png 2017-07-21_LeMaitre_74052_74219.png 2017-07-21_LeMaitre_74220_74711.png 2017-07-21_LeMaitre_74712_75683.png 2017-07-21_LeMaitre_75684_76955.png 2017-07-21_LeMaitre_76956_77519.png 2017-07-21_LeMaitre_77520_77879.png 2017-07-21_LeMaitre_77880_78503.png 2017-07-21_LeMaitre_78504_78719.png 2017-07-21_LeMaitre_78720_84575.png 2017-07-21_LeMaitre_84576_84647.png 2017-07-21_LeMaitre_84648_84695.png 2017-07-21_LeMaitre_84696_84811.png
Thank you. I think it's perfect time for an app, so I'm very grateful to see you here. I'd like to thank the organizer for the invitation. I think it's a very great event and I was very pleased to be here this week. So my talk will discuss how to perform some global sensitivity analysis in a stochastic system and perhaps to motivate a little bit this question and the methods I'm going to introduce. I must say that for many years we have been working on uncertainty quantification but essentially parametric uncertainty. So models having some parameters that are not well known and we developed techniques to deal with the situation, so parametric uncertainty. But in many situations the model itself is not deterministic. So it's not that we have parameters that we don't really know which value to give but the dynamic, the evolution of the system in time is truly stochastic. So we have some inherent stochasticity. So many physical systems have these situations. It can be because at very small scales usually you have noise, typically thermal noise that can appear and so the evolution is really inherently random. So here I'm thinking of molecular dynamic chemical systems for very small size reactors for instance. Reactors could be a biological sense. There are also situations where there are some forcing acting on the system we want to simulate and there's no way we can predict what will be this forcing. So this could be for instance a case in finance, it could be the case in many situations in energy, think for instance of wind turbine, it is submitted to turbulent winds, storms or whatever and you cannot predict a priori what condition the machine will have to stand on support. So in engineering stochastic modeling becomes very important in that case. And there are obviously the situation of all models that in order to maintain tractability, computational tractability, we have to reduce, we have to solve at say coarse scales. While we know that at smaller scales there are some phenomena that we are completely neglecting. Sometimes we can average or buy some homogenization techniques upscale these effects but there are situations where we don't know how to perform such upscaling and one possibility in this case is to add some random noise that mimics the effect of the small scale process that we cannot model. That's the case for instance in climate modeling where they add some random noise to force the flow at small scale to mimic the effect of turbulence mixing. So stochastic modeling and here I'm just thinking about physical models but you can think of all models related to social science for instance. So stochastic modeling is very common. I can be very complex. It's a very nice way to deal with things that we cannot prescribe or that we don't know. But the cost of doing so is that usually you have to prescribe some parameters typically. Say I don't know what will be the wind so I decide to model it with stochastic process but now I have to prescribe some parameters say time correlations, intensity of the wind gusts or the six line. So we have now a system that has some inherent stochasticity because this is a modeling choice that we made and some parameter and this parameter usually cannot be measured. They need to be learned or to be inferred from information. If we have some experimental data for instance but they are not certain. So now we have two sources of uncertainty. One that is inherent to the model and the second one that comes from lack of knowledge about the parameter. And this is really the question we wanted to address with my colleague Omar Knew to see how we could propagate and perform some sensitivity analysis in such systems. So I will have two parts in this presentation. In the first part I will focus on the case of stochastic ordinary differential equations. So typically scalar ones, so equations or ordinary differential equations driven by some vener noise. And I will show how we can consider this equation plus some parametric uncertainty on how to solve this equation by means of polynomial chaos expansion. And the second part I will go to the stochastic simulators and treating this type of system and we will have to introduce new questions that we believe are quite interesting. So I am considering a very simple SDE. So X is a state. Here everything is scalar but you can extend everything to a higher dimension. So vector case for instance. So you have the evolution of your solution X. You have a drift term with a C that may depend on X. It could also depend on T. And you have diffusion terms here with the vener process. So very classical. So what you can do is for instance to simulate this equation to generate samples or trajectories of this system by considering for instance the other schemes. The simplest one you can devise and what you have to do is to draw some IID normal random variables with variance depends on the time step. So you can integrate a trajectory and repeat multiple times going independently realization of this increment. And when you have a set of trajectories you can build some statistic, mean standard deviation, correlation, whatever is needed. And if you consider a... Oh, I thought it was a spider. So if you have a functional of your trajectory, it could be just the value at a certain time or it could be an exit time or whatever it is. So I call it G and now you can estimate the expectation of this quantity of interest X. The same thing is done by simply by sample average, generating M replicas of a trajectory of the system. So the situation I want to tackle now is exactly the same except that now the drift and the diffusion coefficient may have some non-parameter. That could represent special variability. That could be just the magnitude of... if it's just a random variable, it's just the magnitude of the noise that is entering the system. So I did not queue these random parameters and I'm assuming that they have a known probability. So how do you know this probability? Well, perhaps they are coming from some Bayesian inference based on some experimental observation and that would be your posterior of your Bayesian. And the main assumption I'm doing here is that the knowledge of the parameter doesn't say anything about the forcing process as it is. Two are independent, which is quite, I think, natural and I actually don't see situations where this would not be the case. So it's not a strong assumption, however, it is crucial for what I'm going to develop later. Okay, I'm just here making some assumption about square-integrability. So the process has to be a second-order process. If I condition on the random walks, so I still have the random parameter here, so this conditional variance is finite and that is the same if I condition on a certain value of the parameter, the variance is also finite. This equality here are 0.1. Okay, and our goal is now really to say I want to investigate how the variability in X depends on the different source of the variabilities, which are the linear process and the parameter and I want to separate this different effect and being able to draw some conclusion. So the typical approach that people are following consists in conditioning on a certain value of the parameter, then you average with respect to the process and you do that for some moment of X. So typically you're trying to look to the dependence of the mean of the solution of your process as a function of the parameter. So you average with the process and you keep the parameter out and you can do that for a second moment. So you have now the dependence of the variance with the parameter. And now you have functions here that only depends on the parameter and you can apply the classical sensitivity analysis or I mean the parametric one that I've been developing and I will give you a flavor of the idea that are behind. But here is the idea and one possibility is to go through the functional approximation of this two first moments, but you could go to higher moments on a polynomial basis for instance. Okay, and once you have this representation, you can sample as much as you want this approximation and study some statistics that can characterize the dependencies. Okay, now that's the classical approach. So in my view, you are losing some information because the first step you do is always averaging with respect to the venerable process. So there is a loss of information and somehow you're not too penalized because you can do these analyzes for higher and higher moments. But how many moments you need and it becomes very difficult to interpret the outcome of the sensitivity analysis. So in that among different parameters, you can say one is very important for the moment of order five. What does that mean? It's not clear. So what we want to do is slightly different. It's an orthogonal decomposition of the variance. So X has a certain variance and I want to split these variance into different contributions, contributions that is due to the inner noise, so the venerable process, contributions that is due entirely on the parameter and the contribution to the variance of the interaction. And the way we do that is by the soboil of the decomposition of X where the process that I see here as a function of the vener and the parameter is decomposed as a mean value, a function that depends only on the vener, a function that depends only on the parameter and a mixture. And we ask this composition to be, this function is a decomposition to be mutually orthogonal and this makes the decomposition unique. In the case where W and Q are independent, that's where my assumptions become crucial. And here you have the definition of this function. So for instance, the function that depends only on the vener is the conditional expectation of X given W minus C, full expectation. So because this function is orthogonal, you have immediately the decomposition of the variance. You can normalize the partial variance and you get what we call the sensitivity indices. So they are all between 0 and 1. Their sum here is equal to 1. This one characterizes the impact of the vener process on the variance. This one of the parameter and here the mixed contribution. I will illustrate this idea of characterization on an example later. So here is the interpretation of this coefficient and what you can show is that the variance of the mean with respect to W of X for given value of the parameter, so this is really the variance with respect to the parameter of the mean of the solution. This is actually the variance due to W only. Sorry, to Q only. So clear interpretation of what we have classically. However, if we average with respect to the parameter, the dependence of the variance to the parameter, we see here that we have the sum of two contributions. So the contribution due to W to the vener process and also the mixture. And actually there is no way that given this term here, you can split into these two contributions. That's why we believe that the approach we are proposing is richer. We have a finer analysis in particular. If you discover that SW is really large and dominant compared to the other term, there is very few chance that you can learn your parameter by observation. This is kind of important conclusion that you can do. So how to compute this coefficient? What we propose to do is then because we see the process as a function of W on Q, it is to us quite natural to have a series representation. So I'm considering basis of a complete orthogonal set of function psi alpha, which are functional in Q. Actually, I introduce here the measure of the parameter PQ. So there are an orthogonal family. And we are expanding the solution X in this form. So the coefficient in the expansion, the X alpha are again stochastic processes. But they don't depend on the parameter. So I have this expansion. Now I can plan this expansion in my original equation. I will get a residual. If I truncate my expansion, I get a residual. So I'm requiring this residual to be orthogonal to the representation basis. So to every of the psi alpha I have in my expansion. So I'm forming here the garb of projection of my problem. And I end up with a system for the mode or the coefficient of the expansion. And you will notice that this equation, so if you gather all the coefficient X beta in a vector, you end up with an equation that exactly the same structure as previously, except you were a scalar equation, it's a vector equation. But otherwise you can still recognize here a drift term and here a diffusion coefficient term. And what's quite special is that W is still unique. So if you started with W scalar, you still have the same stochastic process that acts on all the modes, force all the modes. Okay. So I'm not going to detail how this garb of machine is implemented or how it is solved. But basically instead of having a single ODE solved by a pair of schemes, now you have a system of ODE, a coupled system. It's not very complicated. Now based on the expansion coefficient that I have in the series here, if I take for, say, convention that the first of this function is a constant. So X0 will be the mean of the solution, well, mean with respect to the parameter. And now I can have the expression of the different function in the Sobolevding expansion and the expression of the different partial variance, which basically if you want to compute, say, the partial variance due to W, only you have to compute the variance of X0. Okay. So basically for each trajectory, you average with respect to the parameter and you look to the variance of the result. But this is made quite easily thanks to the polynomial chaos expansion. Well, maybe I didn't mention, but typical basis we use here are polynomial chaos. So it's polynomial functions. Okay. So let's take the first example. So linear additive system. And I have two uncertain parameter. One is Q1. It appears in the drift term. So it's basically the location of the attracting point. And for the diffusion, I have another parameter Q2 that appears here. If I make nu equals zero, I have an additive noise problem. And if nu is different than zero, we have a multiplicative noise problem. Okay. So I will start with nu equals zero. And I will consider also that Q1 and Q2 are independent. So in the case where we have an additive model. So nu equals zero. You can show that you need just a degree one polynomial expansion to exactly capture the solution. Okay. And here are trajectories of the mean mode, X0. So you, it really looks like typical trajectories for the system or you process. Here is trajectories for the mode that is linear in Q1. So due to the drift term, and you see that there is absolutely no variability. Because basically the system has a completely deterministic evolution. If you fix Q1, plus some noise that comes and had, but they don't mix. Okay. So if you fix Q1, so here you have a range and we respect to W. And well, this is the average we suspect to W. That's all you have. And this is the mixed term, sorry, the term related to the magnitude of the forcing. So this one obviously has an impact on the solution. Okay. If we have a multiplicative noise term, now it's a bit different because the two parametric sources can interact. So now you have variability on the linear term in Q1. And obviously as before, some stochasticity in other mode related to the linear term in Q2. You can even go to higher order, but very quickly, due to the structure of the system, you see that there is no variability in this mode. So you could even, if you were smart enough to analyze the structure of the galerking system, you could drastically reduce the size of the basis. There's absolutely no need to simulate these trajectories. There is no variability at all. Okay. So these are just densities of the coefficients in my expansion to show that some of them are really Gaussian, as you would expect for Gaussian forcing, but some of them are not Gaussians because this is due to the interaction between the modes. Here are also some protection of coordinates in planes to show how there exists some strong dependence between the coefficients. So how the galerking system is really coupling the dynamic between the different modes. Here's perhaps the simplest plot you could imagine to illustrate the idea. What we have here is a single realization of the Brunyan. We have fixed the noise and we are varying the parameter. Okay. So you see here's different trajectories for different values of the parameter, but you observe that they are very similar. I mean, they are shifted. They can be enlarged, but they are quite self-similar. Okay. And this is really this similarity that allows to have a polynomial expansion that converges very quickly. Here's the opposite. We have fixed the parameter to a certain value and we are running different values. And you see that the stochasticity is clearly very much different. Okay. But now we can decompose. So here is a set of trajectories and you can decompose into the function that depends only on the parameter. So essentially we have removed all the noise and we just have variability with the parameter. We can do the opposite. That is, we look or we average over the parameter and we have here variability that is due to WLE. And here's a mixture. Just for illustration. Okay. This one is not that much important. And here are the, not the sensitivity coefficient I have introduced before, but the partial variance directly. So it's not normalized because the variance varies a long time. So it will be difficult to report everything in terms of normalize. So if you have absolutely no variability or no uncertainty on the parameter, so they have variance that are both zero, then all the variance of the system is due to the noise. Okay. That's only W as a source of variability. Now if you have variability also in the drift term, but not in the magnitude of the diffusion, the total variance in red here is essentially due to the variance due to the noise and you have a contribution of the parameter, but you have no interaction. That's what we have seen before where the old trajectory was exactly the same. So this is what is happening here. On the contrary, if you have here a determinist, well, a perfectly known drift term and some uncertainty on the magnitude of the noise, what you have now is a contribution of the, of the winner, of course, but also of the mixture. Okay. And now if you have uncertainty in all the parameters, so drifts and diffusion, now you have all the contributions that are non-zero. So with this way, we can really distinguish what is the inner run variability, what is due to the parameter and have a very fine understanding of what's going on. Okay. Here's another more complicated system. It's, I have just changed the drift function. So the system has two attracting branch. Okay. And we are starting from a point somewhere in the middle and depending on the realization of the Brunyan, we are led to go to a first attracting branch or to another one. Of course, there is a non-zero, but very low probability that when you have selected a branch, you could also branch to the other one, but this is a very rare event. Okay. And for this system, I'm considering some uncertainty which are on the initial location. So clearly if I'm starting from the very high point, it's very likely that I'll end up on this branch and not on the second one. So I have uncertainty on the starting point and I have also uncertainty on the magnitude of the stochastic force. Okay. So here is just an illustration of what is happening depending on my starting point. So if I start from the left, I always go on the same branch or essentially always go to the same branch. If I stop somewhere on the unstable point, I may select with equal probability one branch or the other. And as a contrary, if I start on the other side, I tend to always end up on the same branch. Okay. Here is the effect of the noise magnitude. So I'm starting always at the same point. If I have a very low stochastic forcing, basically I always go into the same branch. I mean it's like deterministic attractor while I stay on the deterministic attractor. And if I increase the noise level, I am able to change this picture. Okay. But of course we would like to do this in a more quantitative way. Okay. Although we understand perfectly what's happening in the system, we want to provide numbers. So as you've seen, even for a given starting point, depending on the brilliant or for a given brilliant, depending on my starting point or intensity of the noise, I can either finish on a branch or another, meaning that for a large time, my solution can be discontinuous. I have a branching process. I have a bifurcation. So the polynomial basis I was using before are not going to work because or I would need an extremely high polynomial degree to capture the discontinuity. So we have here to use some multi-weblet basis with some activity in the parameter space. So this is just to illustrate here. You have the mesh of the parametric space. And we have here, I think, initial condition and noise magnitude on the other direction. And what is represented here is the size of the element over which we have a polynomial approximation. Okay, and here is the result. So you have the evolution of the total variance. What you can conclude is that initially it's essentially the variance that is due to the parameter that dominates. Okay, essentially it depends where you start before the noise can have time to act. So as time increases, you have more and more contribution of the noise, but it's not very much important. Okay, and later on there is the mixed contribution that emerged, which in fact represents a case of events where without noise you would have uncertainty. Only the blue one with the noise is just fluctuation around the branch, so to say. And here the fact that the noise can make you switch from branch to another, which also depends on the parameter. Okay, this is exactly the same, but for higher noise level, and you see that now the magnitude of the different contributions tend to be more balanced as before, because you have more noise in the system. It's quite what you would expect. Okay, as a final note, here it was made with a gallerking projection. This can be difficult if you have a complicated nonlinear drift function, perform the actual projection of the system, but you can also proceed in a nonintrusive way by selecting some value in the parameter space, computing the trajectories for these specific values, and then construct by say regression for instance, or quadrature, the coefficient of the expansion. And you can do that, and perhaps this is a possibility to better illustrate the interest of going through this polynomial chaos approximation. You see that in this case, the coefficients here, or to advance, or to compute the coefficient for a given realization of the vener, you would have to solve NQ trajectories for NQ different values of the parameter. The complexity here is from one trajectory, you go to NQ trajectory, but if your approximation here is sufficiently converged, then you can query your approximation for any value of the parameter, and so you can resample as much as you want. And if most of the variance is carried by the parameter, as opposed to the vener, this can be extremely effective. But there remains a lot of work to really analyze and come up with algorithms that somehow balance the parametric error, sampling error, and do something very effectively. Okay, one final comment here. Here I'm working directly on the representation of the analysis of the whole trajectory, and there is a good reason. For instance, you may want to apply this technique for the case of quantities of interest like exit time, but exit time are not smooth functions, not even with respect to the parameter. You could imagine that for a given realization of the noise trajectories that goes and exit, you change the value of the parameter, you may be missing the exit boundary and continue and exist much later. So there is really an issue with the regularity. So working directly on the expansion of this function, it will be extremely difficult because you cannot construct such approximation for G directly. However, you can use this trick to approximate the trajectory and resample, and by resampling this trajectory, you can very easily generate many, many samples of the exit time. Second part I want to discuss is the case of stochastic simulators. So Raoul discussed this, and it was a work with actually Alvaro Morais, who was a former PhD student of Raoul. So it's a different class of system, systems that are governed by a master equation. If your background is chemistry, that's how you know it. So basically what it says, it specifies the advancement in time of the probability of your system to be in some state x as a certain time t, knowing that it was in a state x0 at time t0. And on the right hand side, you have basically two contributions, the contributions that makes you leaving the state x at time t, and the contributions that makes you going to that state. So you're receiving some probability mass from elsewhere, and you are dumping probability mass by going elsewhere. So this is really a balanced equation. So in our case, x will be on the lattice. It will represent typically a number of molecules in a certain system, a molecule of different species. And this balance equation here has a sum over k reaction channel. So there are k ways for the system to evolve. So in the chemical system, you have k reaction. The function Hj here are the propensity function that basically fix the time rate at which a certain reaction j occurs. And they are generally dependent on the state. And the mu g here are the state change vectors. So each time you have a reaction of time j that occurs, you have to update the state of your system, and that's the quantity by which the state is changing for a single event. So it's a Markov process, and it's used in many, many different applications. But we have to solve this type of equation, while there are plenty of techniques to solve this equation, but what we consider is the Gillespie algorithm, which is an exact approach for this system. What you do is that you really keep in time, or keep track in time of the evolution of the system by jumping from a reaction to the next. So if you have k reaction with current propensity function Hj, you know that in the next time interval t plus t plus dt, this is the probability of the next reaction to work. And it's exponentially distributed. So the more reaction you have, the more likely it is that it occurred in a small period of time. So the Gillespie algorithm is quite easy to understand, so you start from your initial state, x0, then you draw the time to the next reaction, given that you are in a certain state, you can evaluate the rate of reaction here. You decide of the time to the next reaction, so you will advance your system to this time. And you pick at random the reaction that fires first with relative probabilities that are given by the propensity function. So the ones that have large propensity function are more likely to have occurred first than the ones that have low probability function. And then you update your system, and you continue till you reach your final time. So it's really mimicking all reaction in the system. So our goal here is to say, what if I have some uncertainty or some parameter in this propensity function, how can I deal and perform the same type of analysis that I had before using these type of simulators? But, well, actually we had absolutely no clue on how to do that. So the first question we wanted to address is, even for this simple problem, so I can generate multiple trajectories, I can estimate the variance of the population size at some time t, but I'd like to be able to answer and say, if my question is among this k reaction, which one are responsible for most of the variability? So I have no parametric variability, but already I'd like to address this question, which are the reactions that are causing the most of the variance. And surprisingly, we didn't find any answer to this question, which seems quite important, but no. What people are doing is generally they average, they look to the average trajectory and eventually they look to the sensitivity of this average trajectory with respect to the coefficient here, but it's really direct. So how can we do that? So here's a typical example of a set of trajectories. You have here two reactions, how to distinguish the effect of the two is really unclear. So let me just generalize a little bit more the symbol of the decomposition. So I have a functional that depends on n independent random functional, and here I'm just writing the orthogonal decomposition that I had before, but it is generalized. So I have the random input in my functional, here's just an ensemble notation for the symbol of the decomposition. And here it's the definition of the orthogonality, and here's a variance where you sum about all possible, well, all functions in your decomposition, except for the first one, which is a mean. So how can you define this decomposition? Well, here is the definition of this partial variance. So if you consider a set of the direction specified by the vector u, then the variability due to this vector is simply the variance of the expectation of f given this stochastic input minus the sum of all the variance of the functions that contain this set of indices. And you have here the other definition. And from that you can define two important indices, what I call the first order indices, that is the part of the variance that is due to only a certain subset of input, and only this subset. And then I have the total variance, oops, sorry for that, here, so for subset of variable or input in u, u is the index of the input that I want to characterize. So the total variance is not only the variance that is due to this input, but also this input and their interaction with any other input. And I think that in this situation where I have independent input, I have a nice way, a Monte Carlo approach that was proposed by Subaul, to estimate this partial variance and then this sensitivity indices, so first order and second order. The computational effort you have to do is basically compute correlations between your stochastic output when you keep part of your inputs exactly the same, so you generate exactly the same replica for part of the inputs and you redraw the inputs that you want to analyze. So basically for each coefficient you want to compute, you have an additional trajectory to solve, so you double, but the cost is eventually scaling linearly with the number of indices you want to compute. So if you just want to separate or if you have k-channel in your system, you will have to increase as compared to say estimate classically by Monte Carlo moment or an expectation, say now you have to estimate end time or as many input more trajectory, but it's still linear. Okay, but how to apply this to the stochastic simulator? Here in this Gillespie algorithm it's very unclear what are the stochastic inputs on the channel, but more exactly when I draw a random number, I can hardly say okay, this is a number drawn for certain reaction channel. Actually, if I change the indexation of the channel, I may completely change the dynamic. So how can I decide of things or inputs that are staying the same for some channel and others that are varying for the channel? And this was really smart idea by Alvaro Morais was to say okay, what if you go to the Kurt's representation where you represent here the evolution of your system in terms of independent standard Poisson processes because now each channel in this representation has its own Poisson processes. So you can repeat them or you can resample them and you can decide to keep the same sum of them and we compute others. And then you will be able to compute some correlations. And of course, it doesn't mean that by conditioning on a certain realization of Poisson process for a certain channel, you don't change the reaction time of this channel because everything is linked. Okay, so here is typically you have to compute some conditional expectations. So you may fix the value of some of the Poisson processes associated to a subset of channels and you can resample the other one. And then you can compute the correlation to get the sensitivity in this. So if you already have this representation, Kurt's representation, you already have all you need to answer the question of what channels are important in terms of variables. In fact, this is the key position in some of Poisson processes, variable processes. It's done under which assumptions on the coefficients Ae? I don't think there is any particular... But this is... We can make it for any... Any type of systems. Any type of systems? Yeah. Okay, without any assumption on the coefficients Ae? I mean... No, no, no. Those on the intensities of the germs, in fact. No. Now in this problem, you don't have any restriction on our assumption. Okay. It's also known as the time change representation. Okay, it's Tg, so it's a scale... It's Tg here, so the Poisson here are normal, but the argument here depends on the... I see, it's linked to the integral of Ae. Okay, thank you. Okay, so let's see what the conclusion of such analysis is. So here I'm considering the simplest system, so the death birth process. So I have individuals that are with birth rate B. It's independent of the state because it's spontaneous generation. It's not couples that are making children. No, it's just... So you have a birth that occurs once in a while at a rate given by B, and then people are also dying, but at a rate D that is linear in the size of the population, so it's D times X. The more individual you have, the more frequent is to observe a death, or people dying. Right? Okay, so if you apply the Gillespie or the Kurtz algorithm, you can generate trajectories. There is an invariant measure for T going to infinity, which is actually a Poisson distribution that you have on the right. Okay? And now what we would like to say in this variance that we may have, say, at T equal 8, is it B that is more important? Is it my coefficient D, or is it the birth or the death channel? Actually, I should say. Okay. Do you have ideas? I had no idea. Okay, so here are the results. So at very early time, all the variability is due to the birth process, because we start with a zero population, so if there is no individual, no one is going to die, and there will be no variability on the number of dead people, right? Okay. So initially, it's all the birth process that carries the, sorry, here, the variability, and as time is increasing, variability due to the death channel is building up. The sickness of this line here represents the distance between the first order sensitivity indices and the total order sensitivity indices. So the sickness of the line represents the interaction between the two channels. Okay? And actually, if you run the computation over a very long period of time, okay, here it's a log scale, you see that eventually the contribution of the death channel is going to die, okay? Because if you don't have variability in the birth, well, nobody remains. Okay? So over a long time, you need to have variability in the death to have some, in the birth to have variability in the death, okay? And eventually, everything will be due to the birth channel and interaction with death. Okay? Now, a more complicated system, this is a Schlogel system. It's just four reactions, well, actually two reverse reactions. We have three species, B1, S, B2, but B1 and B2 are considered in large excess and are considered constant. So we just modeled the evolution of the number of individual species S that we call X. Here as the propensity function, so you have two body reactions, three body reactions, et cetera. And here are some trajectories. For this system that represents the infection of a cell by a virus, you have basically two asymptotic states. One is the virus has replicated a lot, okay? Or the virus is not replicating and remain at a very low level. So infection is low or infection is high. And you have here the density, not the histogram, more exactly at T equal 8 of the solution. Okay? So for this system, when we carry the analysis we had before, so we want to characterize the effect of the four channels that we have, this is what we obtain as a variance is essentially due to the first and second reaction, while the second and third are much less important. We see that we have strong interactions. But the main conclusion is that first and fourth reaction are the dominant one in terms of inducing variability. And how can this be explained? It can be explained by the fact that it's essentially these two channels and their variability that plays on the selection of containing a low level of viruses or giving birth to lots of replicas. So selecting the upper branch or the lower branch is essentially dictated by channel one and four. While the two other channels induces some noise, but they are way more, way less likely to induce a complete shift of branch. That's what you have here in this plot. So here it's on each of these small plots you have a draw of N1 and N4, so the process of channel one and channel four, and we repeat for other realization of channel two and three. And you see that most of the time you stay on the same branch, it's quite unlikely that you shift. While on the other case, which is completely the opposite, you see that you tend to be less selective in your branch. Okay, that's another system, but I have a taxi to take, and I don't think you would be very much interested into the details of these results, but it's just too stress that you can have a quite fine understanding of how the different channels are combining their stochasticity. It's quite obvious that now we are, and perhaps I even have this, you can add some parameters in the propensity function and you consider them as additional source, independent source of uncertainty. So this is a case for instance where you add in the propensity function an uncertain parameter. This is a situation where you draw different question processes. This is a picture you had before, but it's always the same parameter. While here, it's the same question processes, but you change the value of the parameter. Okay, and again, the trajectories are more similar. Okay, and now you are able to separate and say, okay, my sensitivity indices have a contribution due to the parameters, that's a blue one here, due to the channel, stochasticity, and one that is due to the mixture. And obviously in that case, it would be quite easy to, I think, to learn the parameter because our uncertainty on the parameter is a dominant source of variability, so we should be able to learn quite easily these parameters. The opposite, if we had here a yellow contribution that would be far up, why I would say that you don't really care about the actual value of the parameter because it doesn't give you much variance, and it will be very difficult to learn them. Okay, I think maybe I would prefer to answer questions rather than commenting these slides. Any questions? Like when you are discussing the sensitivity analysis for the stochastic ordinary differential equation, you have like capital Q1 and capital Q2, and you assume that they're independent and you've done the sensitivity analysis in that way. If you just have an idea, if you don't have this independence between this Q1 and Q2, how it's changing, it makes it easier, the sensitivity analysis. Or to consider them independent. If I consider them dependent. Dependent? Yes. How it changes these sensitivity analysis? Well, first you could still continue to perform a separation between inner and stochasticity and parametric one, but it would have to be globally the global contribution of the parameters. Now, if you want to have a finer analysis and distinguish between the contribution of Q1 and the contribution of Q2, then if they are independent, you can apply the symbol of being the composition. If they are not independent, you cannot do that anymore. So there are some, recently they were proposed some techniques to perform this kind of variance decomposition in the case of dependent variable, and perhaps you could reuse them. But the interpretation of the sensitivity coefficient that you get at the end or the characterization that you get is not very well understood at the moment. Typically, you have to account between correlation between parameters, and you can have sensitivity coefficients that can even become negative, which is really strange, right? So thank you Olivier for the presentation. So to come back to the sensitivity analysis, so here you have presented something related to path 5 sensitivity, where you conditioned one motion and so on, but you could also perform a sensitivity analysis regarding the PDE underlying the stochastic differential equation, so that is we know by the distribution of linear PDE. So you could perform also the PCE on these PDE. So could you comment a bit about the choice of doing one way compared to the other? So what is your motivation? There are some, well first I wanted to work by Swise, which is a good motivation. One issue with trying to work directly on the densities or on the master equation is rationality. First, I mean, systems that we want eventually to solve are high dimensional, so solving this master equation in high dimensions is very complicated. That's one point. Second point is that the approximation techniques we are using, so polynomial chaos, honestly they are not that good at enforcing positivity. So expanding density can be problematic. That's another motivation. And basically for us it was making more sense to say, okay, these simulators are existing, stochastic simulators or your favorite tools. Let's try to reuse this machinery as much as possible. And as you've seen, this is actually what we are doing. Not with the gallerking, but it's, gallerking projection is my data, but you could do it non-intrusively and then you just need to have a solver, which I think is a very strong point for users. I'm not very much aware of people working on the direct solution of the master equation. But honestly, positivity is a big issue. Okay, so I think it's time for you to leave. Thank you. Thank you very much.
perform some global sensitivity analysis in a uncertainty. So models having some parameters that are not well known and we developed parametric uncertainty. But in many situations molecular dynamic chemical systems for very small maintain tractability, computational inherent stochasticity because this is a modeling parameter usually cannot be measured. They need sources of uncertainty. One that is inherent to the model and the second one that comes from lack presentation. In the first part I will focus on the stochastic simulators and treating this type of system and we will have to introduce new what you can do is for instance to simulate this this system by considering for instance the other estimate the expectation of this quantity of forcing process as it is. Two are independent, to do is slightly different. It's an orthogonal variance and I want to split these variance into different contributions, contributions that is decomposition of the variance. You can normalize characterizes the impact of the vener process on the mixed contribution. I will illustrate this and what you can show is that the variance of the of what we have classically. However, if we compute this coefficient? What we propose to do system of ODE, a coupled system. It's not very complicated. Now based on the expansion the solution, well, mean with respect to the you have to compute the variance of X0. Okay. So we use here are polynomial chaos. So it's expansion to exactly capture the solution. Okay. the system has a completely deterministic illustrate the idea. What we have here is a decompose. So here is a set of trajectories and removed all the noise and we just have so they have variance that are both zero, then you have uncertainty in all the parameters, so like to do this in a more quantitative way. Okay. this is just to illustrate here. You have the and more contribution of the noise, but it's not gallerking projection. This can be difficult if perform the actual projection of the system, but coefficients here, or to advance, or to compute complexity here is from one trajectory, you go to for the case of quantities of interest like exit time, but exit time are not smooth functions, not approximate the trajectory and resample, and by Second part I want to discuss is the case of governed by a master equation. If your background Markov process, and it's used in many, many consider is the Gillespie algorithm, which is an exact approach for this system. What you do is output when you keep part of your inputs exactly each coefficient you want to compute, you have an drawn for certain reaction channel. Actually, if the Poisson processes associated to a subset of of what channels are important in terms of particular... But this is... We can make it for our assumption. Okay. It's also known as the time it's D times X. The more individual you have, the and stochasticity and parametric one, but it moment. Typically, you have to account between correlation between parameters, and you can have