Weibull Distribution韦布尔分布的深入详述（4）源码实践和实例

Franklin

已于 2022-06-13 20:53:26 修改

阅读量5.6k

点赞数 16

分类专栏： # 数学应用 # 概率论和数理统计数学文章标签： matlab python 数据挖掘

于 2022-06-13 09:36:43 首次发布

本文链接：https://blog.csdn.net/yellow_hill/article/details/125246950

版权

数学同时被 3 个专栏收录

53 篇文章 61 订阅

订阅专栏

数学应用

8 篇文章 4 订阅

订阅专栏

概率论和数理统计

8 篇文章 5 订阅

订阅专栏

前言：

本章我们将利用算法工具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公式为：
$\large\displaystyle F(t) = 1 - e^{-(\frac{t}{\eta })^{\beta }} (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章

Franklin

关注

16
点赞
踩
79

收藏

觉得还不错? 一键收藏
打赏
0
评论
Weibull Distribution韦布尔分布的深入详述（4）源码实践和实例

本章我们将利用算法工具Matlab，对常用的几个韦伯分布进行计算分析。本章实例包括显示面板的数据分析、轮胎数据分析和纺线的数据分析。文后附录，转分享韦伯分布的Python 和 C++的代码。这个例子，分析纺线的张力测量。【如果改变参数：下一段代码，展示了连续变换参数的结果】【案，大致解释一下，首先，前面给出的是韦伯分布的双参数值，分别是shap 和 scale，然后，时间段的步长为j0.1,那么就是51步，然后分别实现了PDF，CDF的公式，然后把他们画出来。】我们可以利用Matlab的函数，生成韦
复制链接

扫一扫