4

As I understand, your ultimate goal is to solve an inverse problem (i.e., infer some parameters from given data / observations). To this end, you want to apply Bayesian Inference, which relates the posterior (i.e., the probability distribution of the unknown parameters) to the likelihood (i.e., the probability model of observing some values given the ...


4

I'm in the middle of doing this for myself. I think the most appropriate analogue to the Gaussian would be the heat kernel in hyperbolic space. Fortunately, this has been figured out before: https://www.math.uni-bielefeld.de/~grigor/nog.pdf (also available in a Bulletin of the London Mathematical Society). If you use the standard decay ($e^{-dist^2/constant}...


4

You can do this with any balanced binary search tree data structure by additionally maintaining the total weight below each node. To randomly sample, compute a uniform random number between $0$ and the weight of the root node, and traverse down through the tree until you find the leaf whose range contains the random number. Unfortunately, while maintaining ...


3

Multidimensional scaling (MDS) is one algorithm that tries to find such a set of points. It's usually used for visualization, so it's usually used in 2 or 3D spaces, but the optimization procedure itself can be used in arbitrary dimensions. To get random samples, you can initialize with random vectors before optimizing, though it is not clear what ...


3

If I understand what information you have available, what you want is not possible: the information available to you is not enough to determine the entropy. It's not even enough to approximate the entropy. It sounds like you have a way to sample from the distribution $p(\cdot)$, and you have a way to compute the ratio $p(a_1)/p(a_2)$ for any pair of ...


2

Assuming C++, if you're willing to store duplicate items, a time efficient approach for the unweighted case is to store them unsorted in a std::vector or std::deque, v. You can then efficiently add new items by inserting them at the end of v. To draw a sample, simply pick a random index i into v and use sample = v[i]; v[i] = v.back(); v.pop_back(); For the ...


2

I'm a little late to reply but I found this review to be really informative. One particularly nice feature of Hamiltonian flows is that they preserve volume in phase space. If we take an infinitesimal box of size $(\delta q, \delta p)$ around a point $(q, p)$ in phase space, then propagate every point inside the box forward in time by $\delta t$ under the ...


2

Because the polynomial $ax^2+bx+c$ and $x^2+\frac{b}{a}x + \frac{c}{a}$ have the exact same roots, let's first simplify the problem by only considering the roots of the equation $$ x^2+b'x + c' = 0 $$ where you can think of $a'=\frac{b}{a}, c'=\frac{c}{a}$. This shows that really only two parameters matter. You could then, for example, plot the (real parts ...


2

For the second part of your question (estimation of entropy difference between distributions) you may be able to use the identity $$F = \langle E \rangle - T S,$$ where $\langle E \rangle$ is the average energy, $T$ is the temperature (it is proportional to $\theta$ in $p \propto \mathrm{e}^{\theta E}$), and $S$ is the entropy. For details, see: Jaynes, E. (...


2

If you want to use a single random number, you first need to find "a" solution (one vector w that meets the criterion), then rotate that vector about $\mathbf{v}$ by a randomly generated angle. We have some flexibility when finding that first vector $w_0$ - we just need to find "a" vector at right angles to $v$ and rotate about that axis by $\beta$. We can ...


2

I've seen in your comment that you want a uniform sampling. With the Eigen library, you can uniformly generate at random a unit quaternion: Eigen::Quaterniond q = Eigen::Quaterniond::UnitRandom(); and then convert it to a rotation (orthogonal) matrix: Eigen::MatrixXd M = q.toRotationMatrix();


2

The example you cited appears to be generating random Householder vectors and multiplying them out using backwards accumulation. Another simple thing to do would be to generate a random matrix $\mathbf A$, then compute its $\mathbf A=\mathbf Q \mathbf R$ decomposition and discard the $\mathbf R$ factor. The two LAPACK functions that you need are [geqrf] (to ...


2

I don't think $D$ is updated at all. It is just a temporary value for $D_{exp}$ times uniform random value in [0,1]. As a general remark, RRT means Rapidly-exploring Random Tree. Hence, the algorithm will involve branching and random directions with step sizes. On the left hand slide of the slide there is a column list of the random numbers used. Strike-...


2

Welcome to SE SciComp. First of all, I would suggest using Jupyter so that you have access to IPython and its nice timing magics (see %time and %timeit magic function). These magics take into account that you need to run the code a couple of times to get reliable measurements. If you need even more in-depth comparison, you also need to take into account ...


2

To generate a specified number of random points in a 3D ball, one possibility is to use rejection sampling to have an initial point distribution, then improve the sampling with Lloyd relaxation. Examples with 30, 300 and 10000 points Rejection sampling: Generate x,y,z coordinates using a uniform random number generator, ignore points that are outside the ...


2

What you’re seeing is a basic acceptance-rejection method. $\Delta S$ generated in that way will be uniform in the ball of radius $\Delta S_{\text{max}}$ centered at the origin. Added I imagine you're doing some kind of optimization by random search, in which case the goal is to find a set of spins that extremize some objective function. The way it looks is ...


2

I discuss the method you describe in more detail in this paper (Rackauckas and Nie 2017) as RSwM2. In that paper I am ever so slightly able to detect that it's sometimes doing something wrong, but since it only has issues with re-rejections it isn't that big of a deal. Those 3 methods (RSwM1, RSwM2, RSwM3) are now the basis of DiffEqNoiseProcess.jl and ...


1

One like this one? https://github.com/brandonckelly/CUDAHM What exactly are you looking for, which language, how plug&play and for which kind of problem? Depending on the case stuff like pyMC3 or tensorflow probability might already suit your needs.


1

The previous answer pretty sums up my understanding on this problem. I just want to add 2 solid references on this regard (Both are from an astrophysics context). The paper by Hogg et al provides a pretty hands-on approach while the the survey of Sharma is more of a survey of MCMC analysis usage in astrophysics. I am not from the astrophysics community, but ...


1

If $X$ is a $(n \times m)$-matrix whose entries are independently generated values from the standard normal distribution, then $X(X^{\top}X)^{-\frac12}$ is a uniformly generated random orthogonal matrix. Source. Here is an implementation with Eigen: #include<iostream> #include<random> #include<Eigen/Eigen> #include<ctime> using ...


1

If you want your sampling over the (a,b) interval to appear uniform when viewed on a log-log scale, linearly interpolate the logarithms: samples = 10.^linspace(log10(a),log10(b)). This is basically what samples=logspace(a,b) does anyway.


1

If I understood the first image correctly, you should do a linear fit with $\ln(N)$ being the independent variable and, presumably, $\ln(d)$ the dependent variable. Furthermore, it seems to me that in your code you use $c$ as $\log(d_{AVG})$, which is different from the first image, where $c$ appears to be a constant. (I'm not giving a full answer because ...


1

Find a vector $\vec a$ such that $\vec v \times \vec a \ne 0$. Define $\hat x = (\vec v\,\times \vec a)/||\vec v \times \vec a||$. By construction, $\hat x$ is a unit vector and is orthogonal to $\vec v$. Define $\hat y = \vec v \times \vec x$. The triple $(\hat v, \hat x, \hat y)$ form a right handed coordinate system. Form a second coordinate system $(\...


1

Well, for the visualization of the discrete and real/imaginary character of the roots, you can simply do a 3D plot of $$ \frac{-b\pm \sqrt{b^2-4ac}}{2a} $$ against $a$, $b$, and $c$. This will give you a surface linking parabolas/square roots/hyperbolas on the different ($aOc$), ($aOb$), ... planes where the surface disappears because the roots become ...


1

If random vector $X$ has variance $S$, then $LX$ has variance $LSL^\top$. So generate whatever random variables with mean 0 and identity covariance matrix, then transform it $LX+\mu$, where $\mu$ is your mean vector and $LL^\top$ equals to your covariance matrix. You can find $L$ by cholesky decomposition.


1

The algorithm you have chosen for your optimization problem sounds like a variant of the family of algorithms known as random optimization. I'd like to also add that since your objective and constraints seem to both be linear (please correct me if I'm wrong), your problem can be casted as a Linear Programming problem, which is a well studied problem with ...


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