概述
前言:
本章我们将利用算法工具Matlab,对常用的几个韦伯分布进行计算分析。
本章实例包括显示面板的数据分析、轮胎数据分析和纺线的数据分析。
文后附录,转分享韦伯分布的Python 和 C++的代码。
例一:一个纺线张力的韦伯分布模型和MT应用
1.1 图形的生成
这个例子,分析纺线的张力测量。
- 已知经验数据:参数shap β =2,参数scale η=0.5的案例
【Matlab的源码如下:】
beta =2;eta =0.5;
time = 0:0.1:5;
fx =(beta/eta)*(time/eta).^(beta-1).*exp(-(time/eta).^beta);
ft=1-exp(-(time./eta).^beta);
subplot(211);plot(time,fx,'linewidth',2.5);xlabel('Time');ylabel('PDF');
subplot(212);plot(time,ft,'linewidth',2.5);xlabel('Time');ylabel('CDF');
【如果改变参数:下一段代码,展示了连续变换参数的结果】
clc;clear;close all;
beta = [0.1 0.5 1 2]; eta = [1 2];time = 0:0.1:10;
for i = 1:length(beta)
for j = 1:length(eta)
fx =(beta(i)/eta(j))*(time/eta(j)).^(beta(i)-1).*exp(-(time/eta(j)).^beta(i));
ft=1-exp(-(time./eta(j)).^beta(i));
subplot(211);plot(time,fx,'linewidth',2.5);xlabel('Time');ylabel('PDF');
subplot(212);plot(time,ft,'linewidth',2.5);xlabel('Time');ylabel('CDF');
end
end
【案,大致解释一下,首先,前面给出的是韦伯分布的双参数值,分别是shap 和 scale,然后,时间段的步长为j0.1,那么就是51步,然后分别实现了PDF,CDF的公式,然后把他们画出来。】
1.2 weibull Matlab函数的使用
我们可以利用Matlab的函数,生成韦伯分布的随机数字。
1.2.1 利用wblrad函数生成韦伯分布数据
我们依旧用上面例子里面的参数,β =2,η=0.5的案例
a = wblrnd(0.5,2,[10000 1]);
time = linspace(0,1,10000);
plot(time,a,'.');
xlabel('time');
ylabel('weibull drandom number');
【如果我们变更一下】
a = wblrnd(2,2,[10000 1]);
time = linspace(0,1,10000);
plot(time,a,'.');
xlabel('time');
ylabel('weibull drandom number');
1.2.2 估算韦伯分布的参数
a = wblrnd(0.5,2,[10000 1]);
time = linspace(0,1,10000);
plot(time,a,'.');
xlabel('time');
ylabel('weibull drandom number');
parmhat = wblfit(a);
parmhat =
0.4989 2.0330
【案,验证,结果和我们之前设定结果一致】
1.2.3 绘制韦伯分布的图像
继续输入
wblplot(a)
图形接近线性,说明图像是韦伯分布。
1.2.4 获取韦伯分布的数学期望和均方差
[M V]= wblstat(0.5,2);
M =
0.4431
V =
0.0537
1.2.5 获取PDF
继续输入
fx = wblpdf(time,0.5,2);
plot(time,fx,'linewidth',2.5);xlabel('Time');ylabel('PDF');
1.2.6 获取CDF
继续输入
ft = wblcdf(time,0.5,2);
plot(time,ft,'linewidth',2.5);
xlabel('Time');
ylabel('CDF');
1.2.7 生成Histogram 数据:
histfit(a,100,'weibull');
1.3 weibull Matlab对象的使用
1.3.1 创建weibull分布对象
注册一个weibulll分布的对象,
pd = makedist('weibull','beta',2,'eta',0.5);
histfit(r,1000,'weibull');
然后,用这个对象可以做很多事情
pdf = pdf(pd,time);
plot(time,pdf,'linewidth',2.5);
例二:显示面板分析和源码
本节我们尝试用Matlab的统计工具和机器学习工具,对平面面板显示器的生产数据进行分析。
问题提出:已知平板显示器,经验数据为符合韦伯分布,且,缩放参数scale η=5000 ,形状参数shap β =0.5,求
1 面板寿命>20000hrs的概率
2 面板寿命<10000hrs的概率
3 10000<面板寿命<20000hrs的概率
4 面板报废的数学期望
问题转为求:
1 P(x > 20000)
2 P(x < 10000)
3 P(10000<x<20000)
4 E(x)
由第一章,我们知道故障的累计分布函数CDF公式为:
F
(
t
)
=
1
−
e
−
(
t
η
)
β
(
t
>
0
)
largedisplaystyle F(t) = 1 - e^{-(frac{t}{eta })^{beta }} (t>0)
F(t)=1−e−(ηt)β(t>0)
我们在Matlab中这么做,
crc;clear;
eta = 5000; beta = 0.5;
p2w_below = wblcdf(20000,eta ,beta);
p1w_below = wblcdf(10000,eta ,beta);
p2w_beyond = 1 - p2w_below ;
p1wto2w = wblcdf(20000,eta ,beta) - wblcdf(10000,eta ,beta);
[M V] = wblstat(eta ,beta);
time = 1:100:25000;
subplot(211); fx = wblpdf(time,eta,beta);
plot(time,fx,'linewidth',2.5);xlabel('Time Hour');ylabel('PDF');
subplot(212); ft = wblcdf(time,eta,beta);
plot(time,ft,'linewidth',2.5);xlabel('Time Hour');ylabel('CDF');
计算结果如下:
1 P(x > 20000) = 0.135335283236613 【p2w_beyond 】
2 P(x < 10000) = 0.756883265565786 【p1w_below 】
3 P(10000<x<20000) = 0.107781451197602 【p1wto2w 】
4 E(x)= 10000 【M】
例二:轮胎失效分析和源码
20000个轮胎,满足韦伯分布,历史数据有:缩放参数scale η=5000 ,形状参数shap β =3.5
求:
6000公里的时候,有多少轮胎报废?
crc;clear;
eta = 5000; beta = 3.5;
TotalTires = 20000;
cdf = wblcdf(6000,eta ,beta);
FailTireNumbert = TotalTires * cdf;
time = 1:10:10000;
subplot(211); fx = wblpdf(time,eta,beta);
plot(time,fx,'linewidth',2.5);xlabel('Kilometer');ylabel('PDF');
subplot(212); ft = wblcdf(time,eta,beta);
plot(time,ft,'linewidth',2.5);xlabel('Kilometer');ylabel('CDF');
FailTireNumbert = 1.698740114017983e+04
cdf = 0.849370057008992 【84.9%的轮胎可以保证到6000公里没问题】
用Matlab工具箱分析韦伯分布
Statistics and Machine Learning Toolbox
Probability Distribution Plotter
【直接上图吧,不解释了】
——————————————————————————————————————————————————————————
Matlab weibull 函数简介
wblrnd
生成韦伯分布的随机数值。Weibull random numbers
R = wblrnd(A,B)
参数 A, scale parameter
参数 B,shape parameter
[parmhat,parmci] = wblfit(data)
韦伯分布参数获取
parmhat = 缩放参数,
parmci = 形状参数
wblplot(x)
通过数据X,绘制韦伯概率分布图形.
图形如果是线性的话,说明是韦伯分布,否则不是。
[M,V] = wblstat(A,B)
输入:A, scale缩放参数 ; B,shape形状参数
返回:M,均值数学期望;V,均方差
Y = wblpdf(X,A,B)
输入:A, scale缩放参数 ; B,shape形状参数
p = wblcdf(x,a,b)
输入:X,样本周期空间,A, scale缩放参数 ; B,shape形状参数
WeibullDistribution
韦伯分布对象生成,这里他输入变量必须为,a,b
Matlab weibull 分布的函数列表:
wblcdf Weibull cumulative distribution function
wblpdf Weibull probability density function
wblinv Weibull inverse cumulative distribution function
wbllike Weibull negative log-likelihood
wblstat Weibull mean and variance
wblfit Weibull parameter estimates
wblrnd Weibull random numbers
wblplot Weibull probability plot
——————————————————————————————————————————————————————————
其他代码参考:
【这里COPY了其他开发者,用C++ 和 Python实现的源码,希望对分析有帮助】
韦伯分布python实现
import numpy as np
import matplotlib.pyplot as plt
# define the pdf of weibull distribution
def weib(x, scale, shape):
return (shape / scale) * (x / scale)**(shape - 1) * np.exp(-(x / scale) ** shape)
scale = 50
shape = 1.5
x = np.arange(1, scale*2)
y = np.zeros(len(x)) # [0 for i in range(len(x))]
for i in range(len(x)):
y[i] = weib(x[i], scale, shape)
scale = 50
shape = 2.5
y1 = np.zeros(len(x)) # [0 for i in range(len(x))]
for i in range(len(x)):
y1[i] = weib(x[i], scale, shape)
scale = 50
shape = 4
y2 = np.zeros(len(x)) # [0 for i in range(len(x))]
for i in range(len(x)):
y2[i] = weib(x[i], scale, shape)
scale = 30
shape = 2.5
y3 = np.zeros(len(x)) # [0 for i in range(len(x))]
for i in range(len(x)):
y3[i] = weib(x[i], scale, shape)
scale = 70
shape = 2.5
y4 = np.zeros(len(x)) # [0 for i in range(len(x))]
for i in range(len(x)):
y4[i] = weib(x[i], scale, shape)
plt.subplot(2, 1, 1)
plt.plot(x, y, 'r', label='scale=50, shape=1.5')
plt.plot(x, y1, 'b', label='scale=50, shape=2.5')
plt.plot(x, y2, 'g', label='scale=50, shape=4')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(x, y3, 'r', label='scale=30, shape=2.5')
plt.plot(x, y1, 'b', label='scale=50, shape=2.5')
plt.plot(x, y4, 'g', label='scale=70, shape=2.5')
plt.legend()
plt.show()
————————————————
版权声明:本文为CSDN博主「心态与做事习惯决定人生高度」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
原文链接:https://blog.csdn.net/robert_chen1988/article/details/100785494
韦伯分布C++实现:
public class WeibullDistribution : ContinuousDistribution
{
private class MMSystem
{
private double mean;
private double variance;
public MMSystem(NumericalVariable variable)
{
this.mean = variable.Mean;
this.variance = variable.Variance;
}
private Vector Evaluate(Vector input, Vector output)
{
double num;
num = input[0];
double num2;
num2 = input[1];
double num3;
num3 = GammaFunctions.Gamma(1.0 + 1.0 / num2);
output[0] = this.mean - num * num3;
output[1] = this.variance - num * num * (GammaFunctions.Gamma(1.0 + 2.0 / num2) - num3 * num3);
return output;
}
public void Solve(ref double a, ref double b)
{
DoglegSystemSolver doglegSystemSolver;
doglegSystemSolver = new DoglegSystemSolver();
doglegSystemSolver.TargetFunction = Evaluate;
doglegSystemSolver.InitialGuess = new GeneralVector(a, b);
Vector vector;
vector = doglegSystemSolver.Solve();
a = vector[0];
b = vector[1];
}
}
private class MLSystem
{
private NumericalVariable variable;
private NumericalVariable log;
public MLSystem(NumericalVariable variable)
{
this.variable = variable;
this.log = variable.Transforms.Log();
}
private Vector Evaluate(Vector input, Vector output)
{
double num;
num = input[0];
double num2;
num2 = input[1];
int length;
length = this.variable.Length;
double num3;
num3 = 0.0;
double num4;
num4 = 0.0;
double num5;
num5 = 0.0;
double num6;
num6 = Math.Log(num);
for (int i = 0; i < length; i++)
{
double num7;
num7 = Math.Pow(this.variable[i] / num, num2);
double num8;
num8 = this.log[i] - num6;
num3 += num8;
num5 += num8 * num7;
num4 += num7;
}
output[0] = (double)length + num2 * num3 - num2 * num5;
output[1] = (double)length - num4;
return output;
}
public void Solve(ref double a, ref double b)
{
DoglegSystemSolver doglegSystemSolver;
doglegSystemSolver = new DoglegSystemSolver();
doglegSystemSolver.TargetFunction = Evaluate;
doglegSystemSolver.InitialGuess = new GeneralVector(a, b);
doglegSystemSolver.Solve();
a = doglegSystemSolver.CurrentPoint[0];
b = doglegSystemSolver.CurrentPoint[1];
}
}
private double scale;
private double shape;
private double location;
public override double Mean => this.location + this.scale * GammaFunctions.Gamma(1.0 + 1.0 / this.shape);
public override double Variance => this.scale * this.scale * (GammaFunctions.Gamma(1.0 + 2.0 / this.shape) - Math.Pow(GammaFunctions.Gamma(1.0 + 1.0 / this.shape), 2.0));
public override double Skewness
{
get
{
double num;
num = GammaFunctions.Gamma(1.0 + 1.0 / this.shape);
double num2;
num2 = GammaFunctions.Gamma(1.0 + 2.0 / this.shape);
double num3;
num3 = num2 - num * num;
return (num * (2.0 * num * num - 3.0 * num2) + GammaFunctions.Gamma(1.0 + 3.0 / this.shape)) / (num3 * Math.Sqrt(num3));
}
}
public override double Kurtosis
{
get
{
double num;
num = GammaFunctions.Gamma(1.0 + 1.0 / this.shape);
double num2;
num2 = GammaFunctions.Gamma(1.0 + 2.0 / this.shape);
double num3;
num3 = num2 - num * num;
return (6.0 * num * num * (2.0 * num2 - num * num) - 3.0 * num2 * num2 - 4.0 * num * GammaFunctions.Gamma(1.0 + 3.0 / this.shape) + GammaFunctions.Gamma(1.0 + 4.0 / this.shape)) / (num3 * num3);
}
}
public double LocationParameter => this.location;
public double ScaleParameter => this.scale;
public double ShapeParameter => this.shape;
public WeibullDistribution(double shape, double scale, double location)
{
if (shape <= 0.0)
{
ThrowException.ArgumentOutOfRange("shape");
}
if (scale <= 0.0)
{
ThrowException.ArgumentOutOfRange("scale");
}
this.scale = scale;
this.shape = shape;
this.location = location;
}
public WeibullDistribution(double shape, double scale)
: this(shape, scale, 0.0)
{
}
public WeibullDistribution(double shape)
: this(shape, 1.0, 0.0)
{
}
public WeibullDistribution(NumericalVariable variable, EstimationMethod method)
{
if (variable == null)
{
ThrowException.ArgumentNull("variable");
}
double a;
a = variable.Mean;
double b;
b = 1.0;
switch (method)
{
case EstimationMethod.Default:
case EstimationMethod.MatchingMoments:
new MMSystem(variable).Solve(ref a, ref b);
break;
case EstimationMethod.MaximumLikelihood:
new MLSystem(variable).Solve(ref a, ref b);
break;
default:
ThrowException.ArgumentOutOfRange("method");
break;
}
this.scale = a;
this.shape = b;
}
public WeibullDistribution(NumericalVariable variable)
: this(variable, EstimationMethod.MatchingMoments)
{
}
public override double ProbabilityDensityFunction(double x)
{
if (x < this.location)
{
return 0.0;
}
if (x == this.location)
{
if (this.shape == 1.0)
{
return 1.0 / this.scale;
}
if (this.shape < 1.0)
{
return double.PositiveInfinity;
}
return 0.0;
}
if (this.shape == 1.0)
{
return Math.Exp((0.0 - (x - this.location)) / this.scale) / this.scale;
}
return this.shape / this.scale * Math.Exp(0.0 - Math.Pow((x - this.location) / this.scale, this.shape) + (this.shape - 1.0) * Math.Log((x - this.location) / this.scale));
}
public override double DistributionFunction(double x)
{
if (x <= this.location)
{
return 0.0;
}
return 1.0 - Math.Exp(0.0 - Math.Pow((x - this.location) / this.scale, this.shape));
}
public override double MomentFunction(int order, double x)
{
if (order < 0)
{
ThrowException.ArgumentOutOfRange("order");
}
if (order == 0)
{
return this.DistributionFunction(x) - this.DistributionFunction(0.0);
}
if (x <= this.location)
{
return 0.0;
}
if (this.shape == 1.0)
{
return Math.Exp((0.0 - this.location) / this.scale) * ElementaryFunctions.Pow(0.0 - this.scale, order) * GammaFunctions.IncompleteGamma(1 + order, (0.0 - x) / this.scale);
}
if (this.location == 0.0)
{
return ElementaryFunctions.Pow(this.scale, order) * (GammaFunctions.Gamma(1.0 + (double)order / this.shape) - GammaFunctions.IncompleteGamma(1.0 + (double)order / this.shape, Math.Pow(x / this.scale, this.shape)));
}
return base.MomentFunction(order, x);
}
public override double InverseDistributionFunction(double probability)
{
if (probability < 0.0 || probability > 1.0)
{
ThrowException.ArgumentOutOfRange("probability");
}
if (probability == 0.0)
{
return this.location;
}
if (probability == 1.0)
{
return double.PositiveInfinity;
}
return this.location + this.scale * Math.Pow(0.0 - Math.Log(1.0 - probability), 1.0 / this.shape);
}
public override double GetRandomVariate(System.Random random)
{
if (random == null)
{
ThrowException.ArgumentNull("random");
}
return this.location + this.scale * Math.Pow(0.0 - Math.Log(random.NextDouble()), 1.0 / this.shape);
}
}
上一章参考:
第1章
第2章
第3章
最后
以上就是彩色康乃馨为你收集整理的Weibull Distribution韦布尔分布的深入详述(4)源码实践和实例前言:例一:一个纺线张力的韦伯分布模型和MT应用例二:显示面板分析和源码例二:轮胎失效分析和源码用Matlab工具箱分析韦伯分布Matlab weibull 分布的函数列表:其他代码参考:上一章参考:的全部内容,希望文章能够帮你解决Weibull Distribution韦布尔分布的深入详述(4)源码实践和实例前言:例一:一个纺线张力的韦伯分布模型和MT应用例二:显示面板分析和源码例二:轮胎失效分析和源码用Matlab工具箱分析韦伯分布Matlab weibull 分布的函数列表:其他代码参考:上一章参考:所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复