Databases
StackOverflow2013
Postsdbo
PO
 

Note that there are some explanatory texts on larger screens.

plurals
  1. PO
    primarykey
    data
    text
    <p>Okay, you appear to be pretty confused about several things. Let's start at the beginning: you mentioned a "multidimensional function", but then go on to discuss the usual one-variable Gaussian curve. This is <em>not</em> a multidimensional function: when you integrate it, you only integrate one variable (x). The distinction is important to make, because there <em>is</em> a monster called a "multivariate Gaussian distribution" which is a true multidimensional function and, if integrated, requires integrating over two or more variables (which uses the expensive Monte Carlo technique I mentioned before). But you seem to just be talking about the regular one-variable Gaussian, which is much easier to work with, integrate, and all that.</p> <p>The one-variable Gaussian distribution has two parameters, <code>sigma</code> and <code>mu</code>, and is a function of a single variable we'll denote <code>x</code>. You also appear to be carrying around a normalization parameter <code>n</code> (which is useful in a couple of applications). Normalization parameters are usually <em>not</em> included in calculations, since you can just tack them back on at the end (remember, integration is a linear operator: <code>int(n*f(x), x) = n*int(f(x), x)</code> ). But we can carry it around if you like; the notation I like for a normal distribution is then </p> <p><code>N(x | mu, sigma, n) := (n/(sigma*sqrt(2*pi))) * exp((-(x-mu)^2)/(2*sigma^2))</code></p> <p>(read that as "the normal distribution of <code>x</code> given <code>sigma</code>, <code>mu</code>, and <code>n</code> is given by...") So far, so good; this matches the function you've got. Notice that the only <em>true variable</em> here is <code>x</code>: the other three parameters are <em>fixed</em> for any particular Gaussian.</p> <p>Now for a mathematical fact: it is provably true that all Gaussian curves have the same shape, they're just shifted around a little bit. So we can work with <code>N(x|0,1,1)</code>, called the "standard normal distribution", and just translate our results back to the general Gaussian curve. So if you have the integral of <code>N(x|0,1,1)</code>, you can trivially calculate the integral of any Gaussian. This integral appears so frequently that it has a special name: the <em>error function</em> <code>erf</code>. Because of some old conventions, it's not <em>exactly</em> <code>erf</code>; there are a couple additive and multiplicative factors also being carried around. </p> <p>If <code>Phi(z) = integral(N(x|0,1,1), -inf, z)</code>; that is, <code>Phi(z)</code> is the integral of the standard normal distribution from minus infinity up to <code>z</code>, then it's true by the definition of the error function that</p> <p><code>Phi(z) = 0.5 + 0.5 * erf(z / sqrt(2))</code>.</p> <p>Likewise, if <code>Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z)</code>; that is, <code>Phi(z | mu, sigma, n)</code> is the integral of the normal distribution given parameters <code>mu</code>, <code>sigma</code>, and <code>n</code> from minus infinity up to <code>z</code>, then it's true by the definition of the error function that</p> <p><code>Phi(z | mu, sigma, n) = (n/2) * (1 + erf((x - mu) / (sigma * sqrt(2))))</code>.</p> <p>Take a look at <a href="http://en.wikipedia.org/wiki/Normal_distribution#Cumulative_distribution_function" rel="noreferrer">the Wikipedia article on the normal CDF</a> if you want more detail or a proof of this fact.</p> <p>Okay, that should be enough background explanation. Back to your (edited) post. You say "The erf(z) in scipy.special would require me to define exactly what t is initially". I have no idea what you mean by this; where does <code>t</code> (time?) enter into this at all? Hopefully the explanation above has demystified the error function a bit and it's clearer now as to why the error function is the right function for the job.</p> <p>Your Python code is OK, but I would prefer a closure over a lambda:</p> <pre><code>def make_gauss(N, sigma, mu): k = N / (sigma * math.sqrt(2*math.pi)) s = -1.0 / (2 * sigma * sigma) def f(x): return k * math.exp(s * (x - mu)*(x - mu)) return f </code></pre> <p>Using a closure enables precomputation of constants <code>k</code> and <code>s</code>, so the returned function will need to do less work each time it's called (which can be important if you're integrating it, which means it'll be called many times). Also, I have avoided any use of the exponentiation operator <code>**</code>, which is slower than just writing the squaring out, and hoisted the divide out of the inner loop and replaced it with a multiply. I haven't looked at all at their implementation in Python, but from my last time tuning an inner loop for pure speed using raw x87 assembly, I seem to remember that adds, subtracts, or multiplies take about 4 CPU cycles each, divides about 36, and exponentiation about 200. That was a couple years ago, so take those numbers with a grain of salt; still, it illustrates their relative complexity. As well, calculating <code>exp(x)</code> the brute-force way is a very bad idea; there are tricks you can take when writing a good implementation of <code>exp(x)</code> that make it significantly faster and more accurate than a general <code>a**b</code> style exponentiation.</p> <p>I've never used the numpy version of the constants pi and e; I've always stuck with the plain old math module's versions. I don't know why you might prefer either one.</p> <p>I'm not sure what you're going for with the <code>quad()</code> call. <code>quad(gen_gauss, -inf, inf, (10,2,0))</code> ought to integrate a renormalized Gaussian from minus infinity to plus infinity, and should always spit out 10 (your normalization factor), since the Gaussian integrates to 1 over the real line. Any answer far from 10 (I wouldn't expect <em>exactly</em> 10 since <code>quad()</code> is only an approximation, after all) means something is screwed up somewhere... hard to say what's screwed up without knowing the actual return value and possibly the inner workings of <code>quad()</code>.</p> <p>Hopefully that has demystified some of the confusion, and explained why the error function is the right answer to your problem, as well as how to do it all yourself if you're curious. If any of my explanation wasn't clear, I suggest taking a quick look at Wikipedia first; if you still have questions, don't hesitate to ask. </p>
    singulars
    1. This table or related slice is empty.
    1. POBest way to write a Python function that integrates a gaussian?
      singulars
      1. PTQuestion
    1. PTAnswer
    1. USali_m
    1. USkquinn
    plurals
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. POBest way to write a Python function that integrates a gaussian?
      singulars
      1. PTQuestion
    1. This table or related slice is empty.
    1. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    3. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    1. CONice response. BTW, in your make_gauss function, I assume that you meant to assign k and s in the body of make_gauss, not in the body of f.
      singulars
      1. PO
      1. USMr Fooz
    2. COThanks for putting the time to make such a complete response. When I said "multidimentional" I refer to passing quad a method that takes multiple arguments, one of which was to be the integrating variable x. I'm integrating finite widths around mu so erf won't do, but I'll use closure in the future.
      singulars
      1. PO
      1. USphysicsmichael
    3. CO@Mr Fooz: Fixed. Don't know how I missed that one.
      singulars
      1. PO
      1. USkquinn
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload