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) # ...
15
There is a straightforward solution with only two passes through the data:
First compute
$$K := \max_i\; a_i,$$
which tells you that, if there are $n$ terms, then
$$\sum_i e^{a_i} \le n e^K.$$
Since you presumably don't have $n$ anywhere near as large as even $10^{20}$, you should have no worry about overflowing in the computation of
$$\tau := \sum_i e^{...
10
To keep precision while you add doubles together you need to use Kahan Summation, this is the software equivalent to having a carry register.
This is fine for most values, but if you are getting overflow then you are hitting the limits of IEEE 754 double-precision which would be about $e^{709.783}$. At this point you need a new representation. You can ...
9
In addition to all that Bill Barth has already said above, let me mention that people often report the fastest of several runs. The rationale is that the actual run time is the ideal run time plus any number of slow downs resulting from other processes running, OS delays, network delays, etc. Since these are all noise we are not interested in, using the ...
8
$\mathbf{A}$ is an $(n+1) \times (n+1)$ matrix. It can be obtained as follows:
$\textbf{A} = \left[ \matrix{1 & 1 & 1 & \cdots & 1 \cr
x_0 & x_1 & x_2 & \cdots & x_{n} \cr
\vdots & \vdots & \vdots & \cdots & \vdots \cr
x_0^n & x_1^n & x_2^n & \cdots & ...
6
Histograms are not useful for high dimensional data. The curse of dimensionality affects one quite fast. As in your case if the grid is of size 7**6, you have on average one point in one bin. Kernel density estimator are better suited as long as you keep the kernel bandwidth large enough. In my experience the top hat kernel as k-nearest neighbor yields ...
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])
6
It's a question of how you choose $\lambda$ and $\Gamma$ (of course). Think for a moment about what happens if you choose $\Gamma=I$ and make $\lambda$ large: in that case you say that it is more important to you to minimize the regularization term $\|Ix\|^2=\|x\|^2$ than to minimize the misfit $\|Ax-b\|^2$. Obviously, making the term $\|x\|^2$ small means ...
5
Knuth (TAOCP Volume 2, 3rd ed., pg. 232) suggests using the formula $M_k = M_{k-1} + (x_k - M_{k-1})/k$ to calculate the mean, where $x$ is a vector of your samples.
See also: a stackoverflow question, an article on accurately computing running variance, and a paper on computing covariances and arbitrary-order statistical moments
5
You haven't specified the distribution of $x(t)$. I'll assume that you want to use a complex normal distribution, since that choice makes it reasonably easy to solve the problem and because this assumption is quite common in signal processing. I'll also discretize the problem so that you're generating a vector $X$ of $N$ entries with a specified complex ...
5
It's actually not all that complicated to calculate the Tracy-Widom CDF just from its definition: see On The Numerical Evaluation Of Fredholm Determinants by Folkmar Bornemann. The Wikipedia page gives the definition as
$$ F_2(s) = \det(I - A_s), $$
where $A_s$ is the integral operator on $[s,\infty)$ with the kernel
$$
\DeclareMathOperator{\Ai}{Ai}
K(x,y) =...
5
Would a decomposition of the form $A = XX^T$ suffice? This would be enough, e.g., if the end goal is sampling from the Gaussian distribution with this given covariance.
If so, you can use the following formula, which is quite similar to your approximation:
$$X = D^{1/2} + \frac{\sqrt{u^T D^{-1} u+1}-1}{u^T D^{-1} u} u u^T D^{-1/2}$$
This follows from ...
4
One solution would be to keep your rolling window of data sorted using a self-balancing binary tree or a hash table. On every update, you pay a constant insertion cost, a $O(n)$ deletion cost, and a $O(n)$ traversal cost to update the quintiles.
4
There are many reasons that a coarser simulation would give different results that a finer grained simulation. A few examples:
are boundary layers being resolved differently?
am I resolving new features (vortices / blockages to flow)
Thinking about a grid result as purely a convolution of a fine-grid result with a Gaussian will work very well in viscous ...
4
Geoffrey already answered the question, but I'd like to add another perspective to it. One of the things you will have to do one way or another is to debug code. It may not necessarily be that you have to debug anything that has to do with the random numbers themselves, but, say, a bug in the function evaluation that depends on the randomly selected sample.
...
4
There is no reason to use hardware random number generation for anything other than full cryptography. For everything else, including computational physics, pseudorandom generators are fine. I would suggest using the Random123 library of Salmon et al.: it's fast, trivially parallelizable, and strong (in particular stronger than Mersenne Twister).
It is ...
4
My experience in these types of problems is that almost any perturbation will do. I suggest you do not try to perturb $u$ or $p$, because generating a divergence-free $u$ with a perturbed component is more complicated than necessary. You can easily perturb $b$ with lots of different things including a small random field or a randomly-distributed set of ...
4
You can easily do this with for loops. Just use n to define your loop limits. Here's a fully functional solution for MATLAB (just define n and x first):
%pre-allocate A
A = zeros(n+1);
%first row:
for j=0:n
A(1,j+1)=sum(x.^j);
end
%rows 2 through n
for i=1:n
A(i+1,1:n)=A(i,2:n+1); %copy from previous row
A(i+1,n+1)=sum(x.^(n+i)); %compute last ...
3
From the Computer Science perspective I do not think that make sense to make an general statistical model for memory access time (latency) and memory bandwidth.
It does make sense to create an statistical model for a algorithm. That is because each algorithm has a specific memory access pattern, the memory access patterns are relevant to the cache hierarchy,...
3
It is imperative to compare many runs with refined meshes until you detect convergence. A single solution with no mesh refinement study shouldn't give you much confidence in your results.
Comparing runs with different fluid properties tells you something different. If you think that a set of runs with different viscosities is relevant to your ultimate ...
3
I'd calculate the mean $\overline x$ and the covariance matrix $C$ of the joint data set, and then do a K/S test on the univariate quantity $V(x):=(x-\overline x)^TC^{-1}(x-\overline x)$ evaluated on the parts. If the K/S test give a significant difference between the parts, there is one. If it gives no significant difference, the test is to be regarded as ...
3
The quantity you are measuring currently is something akin to "prominence" which is a better formed topographical quantity than numerical one. The wiki page in that link describes all sorts of strange cases that arise, even without having any image/surface borders to worry about.
If you want to stick with this alignment metric I would suggest using a ...
3
You don't need to sort your data first, if you have access to additional temporary storage, in which case, you should use a selection algorithm.
3
Testing whether or not the mean is correct, or even if the histogram of your generated random variants "looks" like a certain distribution is not sufficient. Stick with much more rigorous test suites such as TestU01 or Diehard.
Also, you really only have TWO random numbers in each row, because of the constraint that they sum to 1.
This requires more ...
3
You may use Kahan summation algorithm [1]
The idea is to reschedule the sum operations in such a way precision loss is limited. The code is very simple (reproduced from [1] below).
If this does not suffice, you may use multiprecision representations, such as quad doubles [2]. They are supported by several languages / compilers (including GNU c). Finally, if ...
3
Mathematica has the TW distributions:
http://reference.wolfram.com/language/ref/TracyWidomDistribution.html
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 ...
3
(Converted to an answer from my comments and expanded.)
Basically you need to compute the fourth moments $E[x_i x_j y_k y_l]$ for all $i,j,k,l$, given the second moments. These fourth moments are not universally determined, but they depend on the distribution of your variables, even in the one-dimensional case (that's exactly the reason why kurtosis is a ...
2
Did you have a look at dieharder?
Dieharder is a random number generator (rng) testing suite. It is intended to test generators, not files of possibly random numbers as the latter is a fallacious view of what it means to be random. Is the number 7 random? If it is generated by a random process, it might be. If it is made up to serve the purpose of some ...
Only top voted, non community-wiki answers of a minimum length are eligible
