9

You can generate a noise sequence with whatever noise spectrum you want (including $1/f$, also known as pink noise) by generating the noise coefficients in spectral space. The magnitudes of the coefficients should be chosen to give the desired spectrum and the phases should be chosen randomly. You then simply perform an inverse Fourier transform to give the ...


7

I can understand that this property is useful in some applications where the derivative is difficult or computationally infeasible to obtain or does not exist. However, I would not expect such problems to be very relevant in application. If an analytical solution for the derivative is not known, it's very costly and error prone. Calculating the Jacobian is $...


6

The Wikipedia article on this topic is rather bad (from an educational point of view) and somewhat misleading. Do not use it. You'll find a better summary of the method here. In particular, I find it confusing to call $a_k$ a step sequence. To get an intuitive understanding of the method, let's look at Newton's method for the (deterministic) function $M(...


5

Response surface models (a kind of surrogate model) are often used in situations like this. The idea is to sample values of the parameters $a$ and $b$, compute $R(a,b)$ at each point, and then build a regression model (typically quadratic or even higher order) of the function. You can then optimize over the the fitted model.


5

Since $\phi$ is a scalar between $0$ and $1$, the easiest method for finding a root is bisection. If you cannot calculate the expectation of the nonlinear function $$f_\phi(\theta) = \left(\phi r_z +(1-\phi)(r_k+\theta)\right)^{1-\gamma}(r_k+\theta-r_z)$$ in terms of $\phi$ analytically, you can use quadrature to approximate. For the first variant, this ...


4

Programs like Visit and Paraview can do "volume rendering", which is what you show in your figure. You just need to export the data you have in a format that either of these programs can read.


4

I am not an expert in specifically stochastic differential equations, but I would assume that my answer will still be of some value. Computation of the derivative can be challenging, as you mentioned in your question. However, this would be even more pronounced in a multidimensional case, as one would have to calculate Jacobian matrices ($n^2$ entries). So, ...


4

I'd need a TimestepInterval of less than a second - this increases the algorithms computational burden by more than a factor of 10, and is wasteful is food areas with lower RoS. Nonetheless, your existing, simple model is mathematically defensible, easy to explain, and easy to debug. Therefore, I'd recommend that you work on improving the performance of ...


3

You're interested in the solution $u(x,t)$ of the heat equation subject to some boundary conditions and an initial condition $u(x,0)=\delta(x-x_{0})$. You can approximate this at times $t_{1}$, $t_{2}$, $\ldots$, $t_{n}$ by simulating the diffusion of a bunch of particles from $x_{0}$ and constructing histograms of the particle distributions at times $t_{1}$...


3

DifferentialEquations.jl in Julia can do it if you can write it in mass-matrix form. You won't find it mentioned in the tutorial, but you can provide a mass matrix as part of the SDEProblem. Some of the stiff solvers can handle the problem (I see it's not well-documented yet which ones, but it's the symplectic and implicit Euler forms). I will caution that ...


3

There is a huge literature on "stochastic programming", but you're probably interested in what is called "chance constrained programming", in which the constraint coefficients are random variables, and you want to find a solution $x$ such that each constraint is individually satisfied with probability $\eta$ ($\eta$ might 0.95, 0.99, etc.) These problems ...


3

The most straightforward way to solve your SDE is with an Euler-Maruyama scheme. This is a simple and effective method for additive noise, i.e., the diffusion/noise term is not a function of the state, as appears to be the case for your example. Here is some Matlab code to solve your system: n = 2; % Order of system t0 = 0; % Initial time dt = ...


3

There is a nice volume rendering toolbox for MATLAB: http://www.mathworks.com/matlabcentral/fileexchange/22940-vol3d-v2 I think you could tweak it for your purposes.


3

There is no reason to believe that two random fields with the same arithmetic mean would yield solutions that have anything to do with each other. In fact, for the case you consider, one might imagine that maybe the harmonic mean is actually a better indicator, but even that is unclear -- it could also be the geometric mean. Apart from this, you have to ...


3

Achieving bit-for-bit reproducibility of results for a code that runs on a variable number of processors in a distributed computing cluster is extremely difficult if not impossible- there are just too many random factors that can change the result. Even on a shared memory multiprocessor, obtaining bit-for-bit reproducibility can be extremely difficult and ...


2

$$dx_{t}=2x_{t} dt +x_{t} dw_{t}\\x_{0}=1$$ % EM Euler-Maruyama method on linear SDE % % SDE is dX = lambda*X dt + mu*X dW, X(0) = Xzero, % where lambda = 2, mu = 1 and Xzero = 1. % % Discretized Brownian path over [0,1] has dt = 2^(-8). % Euler-Maruyama uses timestep R*dt. state=randn(100) lambda = 2; mu = 1; Xzero = 1; % problem ...


2

The cumulative distribution function is the integral (antiderivative) of the probability distribution function. In other words, the PDF is the derivative of the CDF. You can therefore compute the PDF by computing the derivative of your data, for example by forming a difference quotient to approximate the derivative from a finite set of points.


2

(1.8) is a simple reformulation of the deterministic LP (1.2) It's still a deterministic LP. This may be somewhat confusing since the authors are going back forth between a deterministic LP formulation (assuming that $d$ is known) and a stochastic LP formulation (in which we're minimizing the expected value over random values of $d$.) See (1.9) for the ...


2

Active learning (aka experimental design) strategies are suitable for this. One of them being the response surface modeling suggested by Brian Borchers. The idea is to choose the next point to evaluate to learn about the optimal M as fast as possible. Here are some old and new papers to start from: Harold J. Kushner. A new method of locating the maximum of ...


2

How is the KDE used for calculating the new residual value here? Assume that you have already processed $n$ samples $x_1,x_2,\ldots,x_n \in \mathbb{R}^3$. Then the residual images $r_j : \mathbb{R}^2 \rightarrow \mathbb{R}$, $j=1,2,\ldots$, take the form $$r_j(x) = m_j(x) - \sum_{i=1}^n e_i K(x-proj_j(x_i))$$ where $m_j$ is the $j$-th projection image, ...


2

You are applying Richardson extrapolation to increasingly-accurate independent sample paths of an OU process. How would this work even for independent identically normally distributed random variables $y_1,y_2,y_3,y_4\sim\mathrm{N}(\mu,\sigma^2)$? These have already "converged" to the true distribution $\mathrm{N}(\mu,\sigma^2)$, and you might expect the ...


2

You cannot take deterministic methods and slap a kind of Milstein correction on them to get higher order. Instead, you should be looking at methods which are specifically designed for integrating SDEs at high order if you want good accuracy. If pure efficiency is what matters and the noise is additive/diagonal, then the adaptive Rossler methods mentioned by @...


2

First of all, it is impossible to intertwine a multi-step Runge–Kutta method and the Milstein–Itō methods for a multitude of reasons that go beyond the scope of this question¹. So the best you can possibly do in is: Make a deterministic Runge–Kutta step, ignoring the noise term. Apply noise (in a Milstein–Itō fashion). Go to 1. This has the problem that ...


2

The setup that is used in DifferentialEquations.jl and QuantumOptics.jl is what's known as time-adaptive jumping. It's nice because it allows for jump events to do things like change the number of DEs, and the jumps are computed exactly. However, it does have the limitation that jumps are computed exactly, so if you have a high jump rate then this slows down....


2

If $H = H(q, p)$ is the Hamiltonian of your mechanical system, $T$ is the temperature and $\gamma$ is the friction, then Langevin dynamics (LD), \begin{equation*} \left\{ \begin{aligned} \mathrm{d} q &= \nabla_p H(q, p) \, \mathrm{d} t, \\ \mathrm{d} p &= -\nabla_q H(q, p) \, \mathrm{d} t - \gamma \, p \, \mathrm{d} t + \sqrt{2 \...


2

Take a classic paper like this one from Davie and Gaines on solving the stochastic heat equation. By equation (2) they say We consider finite-difference approximations to (1). The simplest such approximation is the explicit scheme Equation 2 is then really simple: it's just the 3 point [1 -2 1] stencil of the Laplacian for the spatial part with an Euler-...


2

It turns out that stable summation of numbers is a topic that is still being researched today -- but you can get the basics by looking up "Kahan's summation algorithm". That said, there really is only a stability issue if you have numbers of widely varying size. In that case, you need to do the summation in a particular order -- intuitively from ...


1

In a more general sense, you could formulate your problem as follows. Suppose you have a number of resources, i.e. places where you can find CPU clicks (I use clicks but depending on your granularity this could be hours, days, months). Number these resources $i=1, 2, \ldots$. All these resources have a price to pay. As you indicate in your post, these ...


1

Not sure if I understood correctly, but I think that if your parameters µ and $\Delta t$ are constant during all the simulations and for all nodes, you can generate a list of possible $p_j$ values in the beginning of the simulations. Because you have a limited possiblity of $b_j$ values (you have maximum 8 neighbors in a 2D grid if you look at first-...


1

Here's a very simple way to create pink noise in C#, which just sums lots of waves (spaced logarithmically apart) together! The functions outputs a double array with values from -1 to 1. This represents the lowest and highest points in the waveform respectively. The 'quality' parameter represents the number of waves produced to make the sound. I find 5000 ...


Only top voted, non community-wiki answers of a minimum length are eligible