概率密度的黎曼和
Riemann sum of a probability density
我试图找出随机变量事件超过特定值的概率,即 pr(x>a),其中 a 是某个常数,通常远高于 x 的平均值,而 x 是不是任何标准的高斯分布。所以我想拟合一些其他的概率密度函数,并将 x 的 pdf 从 a 积分到 inf。由于这是一个尖峰建模问题,我认为这是一个极值分析问题,并发现威布尔分布可能是合适的。
关于极值分布,威布尔分布有一个非常 "not-easy-to-implement" 的积分,因此我想我可以从 Scipy 得到 pdf,然后做一个黎曼求和。我还认为我也可以简单地评估核密度,得到 pdf,并对黎曼和做同样的事情,来近似积分。
我在 Stack 上找到了一个 Q,它提供了一种在 Python 中进行黎曼求和的简洁方法,我修改了该代码以解决我的问题。但是当我计算积分时,我得到了奇怪的数字,表明 KDE 或黎曼求和函数有问题。
根据 Scipy 文档,根据 Weibull 的第一种情况:
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])
看起来像这样:
然后也试了下KDE的方法
pd = ss.gaussian_kde(x).pdf(x_grid)
我随后 运行 通过以下函数:
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))
就 Weibull 而言,它给了我
>> 0.272502150549
>> 18.2860384829
对于 KDE,我得到
>> 0.448450460469
>> 18.2796021034
这显然是错误的。对整个事情进行积分应该给我1,而18.2+还差得很远。
我对这些密度函数可以做什么的假设有误吗?或者我在黎曼和函数中犯了一些错误
the Weibull distribution has a very "not-easy-to-implement" integral
咦?!
Weibull distribution 具有非常明确的 CDF,因此实现积分几乎是一行(好吧,为清楚起见,将其设为两行)
def WeibullCDF(x, lmbd, k):
q = pow(x/lmbd, k)
return 1.0 - exp(-q)
当然还有ss.weibull_min.cdf(x_grid,p[0], p[1], p[2])
如果你想从标准库中选择
我知道有一个公认的答案对你有用,但我在寻找如何计算概率密度的黎曼和时偶然发现了这个答案,其他答案也可能如此,所以我会试一试。
基本上,我认为你有(现在是什么)旧版本的 numpy 允许浮点索引,你的 pd
变量指向从 pdf 中提取的值数组,对应于 xgrid 的值.如今,当您尝试使用浮点索引时,您会在 numpy 中遇到错误,但由于您没有这样做,您正在访问与该索引对应的网格值处的 pdf 值。您需要做的是用您想在黎曼和中使用的新值计算 pdf。
我编辑了问题中的代码以创建一种用于计算 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
下面的黎曼实现也可以使用(它使用Java而不是Python)抱歉。
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
}
}
我试图找出随机变量事件超过特定值的概率,即 pr(x>a),其中 a 是某个常数,通常远高于 x 的平均值,而 x 是不是任何标准的高斯分布。所以我想拟合一些其他的概率密度函数,并将 x 的 pdf 从 a 积分到 inf。由于这是一个尖峰建模问题,我认为这是一个极值分析问题,并发现威布尔分布可能是合适的。
关于极值分布,威布尔分布有一个非常 "not-easy-to-implement" 的积分,因此我想我可以从 Scipy 得到 pdf,然后做一个黎曼求和。我还认为我也可以简单地评估核密度,得到 pdf,并对黎曼和做同样的事情,来近似积分。
我在 Stack 上找到了一个 Q,它提供了一种在 Python 中进行黎曼求和的简洁方法,我修改了该代码以解决我的问题。但是当我计算积分时,我得到了奇怪的数字,表明 KDE 或黎曼求和函数有问题。
根据 Scipy 文档,根据 Weibull 的第一种情况:
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])
看起来像这样:
然后也试了下KDE的方法
pd = ss.gaussian_kde(x).pdf(x_grid)
我随后 运行 通过以下函数:
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))
就 Weibull 而言,它给了我
>> 0.272502150549
>> 18.2860384829
对于 KDE,我得到
>> 0.448450460469
>> 18.2796021034
这显然是错误的。对整个事情进行积分应该给我1,而18.2+还差得很远。
我对这些密度函数可以做什么的假设有误吗?或者我在黎曼和函数中犯了一些错误
the Weibull distribution has a very "not-easy-to-implement" integral
咦?!
Weibull distribution 具有非常明确的 CDF,因此实现积分几乎是一行(好吧,为清楚起见,将其设为两行)
def WeibullCDF(x, lmbd, k):
q = pow(x/lmbd, k)
return 1.0 - exp(-q)
当然还有ss.weibull_min.cdf(x_grid,p[0], p[1], p[2])
如果你想从标准库中选择
我知道有一个公认的答案对你有用,但我在寻找如何计算概率密度的黎曼和时偶然发现了这个答案,其他答案也可能如此,所以我会试一试。
基本上,我认为你有(现在是什么)旧版本的 numpy 允许浮点索引,你的 pd
变量指向从 pdf 中提取的值数组,对应于 xgrid 的值.如今,当您尝试使用浮点索引时,您会在 numpy 中遇到错误,但由于您没有这样做,您正在访问与该索引对应的网格值处的 pdf 值。您需要做的是用您想在黎曼和中使用的新值计算 pdf。
我编辑了问题中的代码以创建一种用于计算 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
下面的黎曼实现也可以使用(它使用Java而不是Python)抱歉。
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
}
}