Riemann sum of a probability density

user_2625173 0
Viewed 1115

I am trying to find the probability of an event of a random variable past a specific value, i.e. pr(x>a), where a is some constant, typically much higher than the average of x, and x is not of any standard Gaussian distribution. So I wanted to fit some other probability density function, and take the integral of the pdf of x from a to inf. As this is a problem of modelling the spikes, I considered this an Extreme Value analysis problem, and found that the Weibull distribution might be appropriate.

Regarding extreme value distributions, the Weibull distribution has a very "not-easy-to-implement" integral, and I therefore figured I could just get the pdf from Scipy, and do a Riemann-sum. I also thought that I could as well simply evaluate the kernel density, get the pdf, and do the same with the Riemann sum, to approximate the integral.

I found a Q here on Stack which provided a neat method for doing Riemann sums in Python, and I adapted that code to fit my problem. But when I evaluate the integral I get weird numbers, indicating that something is either wrong with the KDE, or the Riemann sum-function.

Two scenarios, the first with the Weibull, in accordance with the Scipy documentation:

x = theData
x_grid = np.linspace(0,np.max(x),len(x))

p = ss.weibull_min.fit(x[x!=0], floc=0)
pd = ss.weibull_min.pdf(x_grid,p[0], p[1], p[2])

which looks like this:

enter image description here

and then also tried the KDE method as follows

pd = ss.gaussian_kde(x).pdf(x_grid)

which I subsequently run through the following function:

def riemannSum(a, b, n):
    dx = (b - a) / n
    s = 0.0
    x = a
    for i in range(n): 
        s += pd[x]
        x += dx
    return s * dx          
print(riemannSum(950.0, 1612.0, 10000))
print(riemannSum(0.0, 1612.0, 100000))

In the case of the Weibull, it gives me

>> 0.272502150549
>> 18.2860384829

and in the case of the KDE, I get

>> 0.448450460469
>> 18.2796021034

This is obviously wrong. Taking the integral of the entire thing should give me 1, and 18.2+ is quite far off.

Am I wrong in my assumptions of what I can do with these density functions? Or have I made some mistake in the Riemann sum function

weibull integral probability-density python
3 Answers

I know there is an accepted answer that worked for you but I stumbled across this while looking to see how to do a Riemann sum of a probability density and others may too so I will give this a go.

Basically, I think you had (what is now) an older version of numpy that allowed floating point indexing and your pd variable pointed to an array of values drawn from the pdf corresponding to the values at xgrid. Nowadays you will get an error in numpy when trying to use a floating point index but since you didn't you were accessing the value of the pdf at the grid values corresponding to that index. What you needed to do was calculate the pdf with the new values you wanted to use in your Riemann sum.

I edited the code from the question to create a method that works for calculating the integral of the pdf.

def riemannSum(a, b, n):
     dx = (b-a)/n
     s = 0.0
     x = 0
     pd = weibull_min.pdf(np.linspace(a, b, n), p[0], p[1], p[2])
     for i in range(n):
         s += pd[x]
         x += 1
     return s*dx
user_26007840

Below Riemann implementation can also be used (it uses Java instead of Python) sorry.

import static java.lang.Math.exp;
import static java.lang.Math.pow;

import java.util.Optional;
import java.util.function.BiFunction;
import java.util.function.BinaryOperator;
import java.util.function.Function;
import java.util.stream.IntStream;

public class WeibullPDF
{
    public interface Riemann extends BiFunction<Function<Double, Double>, Integer, 
    BinaryOperator<Double>>     { }

    public static void main(String args[])
    {
        int N=100000;
        Riemann s = (f, n) -> (a, b) -> 
        IntStream.range(0, n).
        .mapToDouble(i->f.apply(a+i*((b-a)/n))*((b-a)/n)).sum();
        double k=1.5;
        Optional<Double> weibull = 
        Optional.of(s.apply(x->k*pow(x,k-1)*exp(-pow(x,k)),N).apply(0.0,1612.0));
        weibull.ifPresent(System.out::println); //prints 0.9993617886716168
    }
}
user_100404150
Related