13

This is not really 4D data. As Geoff said, it's 3D scalar data, i.e. you're visualizing a scalar function of three variables: $f(x,y,z)$. There are several ways to visualize this kind of data, and many tools that will help you. I'll show you a few styles of plots you can make. Contour plot showing one or more $f(x,y,z) = \text{(const.)}$ surfaces, ...


9

The traditional approach for scalar field-based data (temperature, velocity magnitude, pressure, density, etc.) plotted over two or three space dimensions uses color. It's important to note that choice of color scheme can distort your impressions of the data. For this reason, do not use a rainbow color scheme. (For why, see here, here, here, and here.) ...


7

So the way I went about formulating the problem was to essentially write the following equations: The state that will be estimated, which is defined as a column vector, is the following: $$w = [vec(A)^{T},\bar{N}^{T}]^{T}$$ where $\bar{N}$ is the unknown average $N$ vector and $vec(A)$ is the vectorization operator on the unknown matrix A. Based on the ...


6

This doesn't randomly sample points, but instead chooses representative points deterministically. scipy.stats.norm.ppf(np.linspace(0, 1, 1000+2)[1:-1])


5

The initial velocities are drawn from a Gaussian distribution with variance $$\sigma_i^2=\frac{k_{\textrm{B}}T}{m_i},$$ where $k_{\textrm{B}}$ denotes Boltzmann's constant, $T$ is the temperature and $m_i$ is the mass of the $i^{\textrm{th}}$ particle. Thus, the problem boils down to generate random numbers from a gaussian distribution using uniformly ...


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

Since you are assuming $\eta$ is normal what i would do is try to compute the expectation as fast as possible for every $\theta$. So I would compute the expectation using any numerical integration I might know. That is, if $n$ and $N$ are the Gaussian density and distribution respectively, $$ f( \theta ) = \int_{ -\infty }^\infty R(\theta, \eta ) n( \eta ) ...


4

If $X\sim \text{Normal}(\mu,\sigma)$, then the $Y\sim \text{Half-Normal}$ is obtainable via several approaches. Absolute Value: $\quad Y = |X|\quad $ (As pointed out by @DavidZ.) Truncation: $\quad Y = X_{(0,\infty)} \quad $ MATLAB does this nicely with truncate(). Conditioning (Logical Indexing): $\quad Y = (X|X>0)\quad $ Graphical examples below ...


4

I am familiar with the equation in a different application, it is a so called convection-diffusion equation. There is no problem to solve it for negative values of B. The scheme you mention (Chang-Cooper method) is known in my field as the exponential upwind scheme and it is valid for all nonzero values of B. Due to a lack of time I put here screen shots of ...


4

You could use a least-square solver with bounds. import numpy as np from scipy.optimize import lsq_linear W = np.array([ [1., 0., 0., -1.76690464, 0., -1.76690464, 0., -1.76690464, 0., -1.76690464], [0., -38.43501272, 1., 0., 0., -38.43501272, 0., -38.43501272, 0., -38.43501272], [0., -41.64051053, 0., -41....


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 ...


3

The assumptions you have are no more specific than what you need to assume to make something like Newton's method work. In fact, you don't even assume enough to make the problem unique: you only assume that $F(x)$ is non-decreasing, when of course you need it to be monotonically increasing to make finding an answer unique. As a consequence, there is really ...


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

Let's first rewrite this. From your formulas, you have that $$ A_{t+1} = A_t^\rho \exp(e_t) $$ where $A_t$ is just the previous value and, consequently, just a fixed parameter. So what you need to compute is then $$ E_t[\phi(A_{t+1},\eta_{t+1})|A_t] = \frac 1C \int_{-\infty}^\infty \int_{-\infty}^\infty \phi(A_t^\rho e^a,b) e^{-a^2/(2\...


3

Computing an approximation of $P_1=e^Q$ is typically done by some kind of polynomial expansion of the matrix exponential. One doesn't take the Taylor expansion (it converges too slowly) but for the sake of the argument, let me assume that you do and that you approximate $$ P_1 = e^Q \approx \sum_{k=0}^N \frac1{k!} Q^k = I+Q+\frac 1{2!} Q^2 + \frac 1{3!}Q^...


3

The Chinese Restaurant Process is a way of looking at a Dirichlet process. It is a distribution over distributions. There are various ways of thinking about either. One way of looking at it is that as you draw samples, for each new sample: there is a finite probability that the new sample is assigned to an existing cluster and otherwise it becomes the ...


3

In general, your type of question would be called a "multivariate goodness of fit test". If $F(x_1,\ldots,x_n)$ is the $n$-dimensional CDF for the theoretical distribution, and the random variables $(X_1,\ldots,X_n)$ are a sample from your ODE, then $$ Z_1 = F(X_1), \quad Z_2 = F(X_2\mid X_1), \quad\cdots\quad Z_n = F(X_n\mid X_1,\ldots,X_{n-1}) $$ are ...


2

A software package mentioned in Moler and van Loan's "twenty-five years after" portion of their paper (Sec. 13, see link in Christian Clason's comment) is Roger Sidje's Expokit, with an accompanying paper that explains the approach. For large sparse matrices a Taylor-like polynomial is applied to Krylov subspaces, obtaining the advantages of matrix-vector ...


2

If your PDF is bounded, you could try approximating its inverse with a high-degree polynomial interpolant. This is usually considered a bad thing, but that's just a myth. Some things to keep in mind: Instead of using an equispaced grid, interpolate at the Chebyshev nodes of the first kind, i.e. $x_i = \cos\left(\pi\frac{2i-1}{2N}\right)$, for $F(x)$ ...


2

To answer your edit... Basically is a sequence of 2^160 unique strings guaranteed to generate every possible combination of a sha1 hash? No, there is almost no chance. The probability is some double exponential like exp(-exp(100)). Say that by some miracle you have seen no collision among your first 2^160 - 1 unique strings. Then your last unique ...


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

I don't think you're missing anything, MCMC is used to sample points from a given distribution, known up to a normalization constant, and to evaluate expectations w.r.t. that distribution (not integrals over the volume). It's often used when the normalization constant of the distribution is not known, and doesn't require you to know it. On the other hand, ...


2

NumPy comes with a nifty random library with various distributions, including normal (Gaussian). From the Numpy documentation: mu, sigma = 0, 0.1 # mean and standard deviation s = np.random.normal(mu, sigma, 1000) which will give you 1000 normally distributed values with mean mu and standard deviation sigma.


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

If you can use C++11 you can use the built-in Poisson distribution #include <iostream> #include <random> static std::random_device rd; static std::mt19937 gen(rd()); int main() { std::poisson_distribution<int> pd(5); for (int i = 0; i < 10000; ++i) std::cout << pd(gen) << '\n'; } I mean, come on, it doesn't get ...


2

Posting a partial answer in hopes of inspiring a full one. For simplicity take $f(x)=x$. Reparameterize to $\xi\sim N(0,I)$. We are then interested in the quality of this approximation: $$ \nabla V(\theta)\approx \nabla \mathbb{E}V(\theta+\epsilon \xi) $$ Now, add an extra assumption that $V$ is $\beta$-smooth. Then we have enough regularity to interchange ...


1

In the loop over your instances of 100 random values for $k$, you could write the value for $k$ into a file, and call your finite element code; the finite element code could then open the file, read the value for $k$ from it and then use it for whatever computation you may want to do with it. As for your second question, visualizing uncertainty is difficult....


1

Mathematically, the probability density function ($\operatorname{PDF}$) for $Z$ is given by the integral: $$\operatorname{PDF}(Z) = \int \delta\left(Z - f_a(X,Y)\right) \operatorname{PDF}(X,Y)\operatorname{d}X \operatorname{d}Y.$$ If the transformation $$\left[\begin{array}{c} X \\ Y \end{array}\right] \rightarrow \left[\begin{array}{c} f_a(X,Y) \\ g_a(X,Y) \...


1

I don't know the answer to your question, but there is almost certainly a large amount of literature on it in the computational biology community. This is because your problem is in essence how genomes are reconstructed: genome sequencers break the many copies of the same DNA (which is a string consisting of 4 letters) into small pieces at essentially random ...


1

By taking differences of the cumulative distribution function $F$, you can find the probability density function $p$. What you're looking for are local maxima of $p$; from your second plot, it looks like there is one big maximum at $t = 7$ and another at $t \approx 1200$. You can spot these from the CDF directly as inflection points, where $F$ goes from ...


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