24

You have to subclass the rv_continuous class in scipy.stats import scipy.stats as st class my_pdf(st.rv_continuous): def _pdf(self,x): return 3*x**2 # Normalized over its range, in this case [0,1] my_cv = my_pdf(a=0, b=1, name='my_pdf') now my_cv is a continuous random variable with the given PDF and range [0,1] Note that in this example ...


15

You should check out sympy.stats. It provides an interface to deal with random variables. The following example provides a random variable X defined on the unit interval with density 2x In [1]: from sympy.stats import * In [2]: x = Symbol('x') In [3]: X = ContinuousRV(x, 2*x, Interval(0, 1)) In [4]: P(X>.5) Out[4]: 0.750000000000000 In [5]: Var(X) # ...


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


10

As Wikipedia says, you can just compute a normally distributed random variable with a mean of zero and take the absolute value. The half-normal probability distribution function is identical to the normal PDF with $\mu = 0$ for $X > 0$, except for a factor of two for the normalization.


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

In principle one could preselect the number of samples to be drawn from each sub-distribution, then visit each sub distribution only once and draw than number of points. That is Find the random set $<n_1, n_2, \dots, n_k>$ such that $n = \sum_{i=1}^k n_i$ and respecting the weights. I believe that you do this by drawing a Poisson distribution a ...


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

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

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


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

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

Note: The original version of this question asked about a "weighted sum of normal distributions" to which the following answer might be useful. However, after a good bit of discussion on this answer, the answer by @Geoff, and on the question itself, it became clear the question was really on sampling a "mixture of normal distributions" to which this answer ...


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

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


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

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

Update: This answer is incorrect, stemming from confusion in terminology (see the comment chain below for details); I'm only leaving it up as a guidepost so that people do not repost this answer (besides Barron). Please do not vote it up or down. I'd just use properties of random variables to reduce it to a single normally distributed random variable. The ...


2

First of all, we can decorrelate the linear and quadratic terms by completing the square: $g + \sum A_i h_i + B_i h_i^2 \equiv g' + \sum B_i (h_i')^2$ The sum of squared nonstandard normal random variables appears to have no name or neat density, so I will compute it numerically, assuming that all the random variables are independent. Example Remember ...


2

Using the assumption of normality, we can get an expression for this in Maple, but it's not in closed form (it still has an integral sign in it). Here's what I did: with(Statistics): h := RandomVariable(Normal(mu[H], sigma[H])); g := RandomVariable(Normal(mu[G], sigma[G])); PDF(g + A * h + B * h^2, z) assuming A > 0, B > 0; Without the assumptions on ...


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

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

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

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


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