solve 3 parameter weibull using method of moments - integral with parameters inside a solve

17 visualizaciones (últimos 30 días)
Hi all,
I am trying to use the method of moment to solve the 3 parameters of the weibull distribution. I have input all the equations however the intergal function embedded gives me an output error. I am not sure how to embed the integral into the solve function. I.e. the integral is complaining that is has a variable and not a value. Here is a link to the euqations for the moments.
Error using integral (line 85): A and B must be floating-point scalars.
Error in MOM>@(n)integral(@(s)fun(s,n),0,n)
Error in MOM (line 21) : eqn=[c +
(a*(g(1+1/b)))==u,a*sqrt(g(1+2/b)-(g(1+1/b))^2)==var,(g(1+3/b)-3*g(1+1/b)*g(1+2/b)+2*(g(1+1/b))^3)/((g(1+2/b)-(g(1+1/b))^2)^1.5)==skew;]
note I use 'name' cause i have multiply files which I would need to apply this too.
name=strcat('Location_',num2str(p,'%02d'),'_hs.mat');
load(name);
u=nanmean(values);
var=std(values);
skew=skewness(values);
s=values;
syms a b c
fun=@(s,n) s.^(n-1).*exp(-s); %integal of the gamma function
g=@(n) integral(@(s)fun(s,n),0,n); %gamma funtion
eqn=[c + (a*(g(1+1/b)))==u,a*sqrt(g(1+2/b)-(g(1+1/b))^2)==var,(g(1+3/b)-3*g(1+1/b)*g(1+2/b)+2*(g(1+1/b))^3)/((g(1+2/b)-(g(1+1/b))^2)^1.5)==skew;]
vars=[a b c];
[sola, solb, solc]=solve(eqn, vars)
  5 comentarios
Lauren Arendse
Lauren Arendse el 20 de Feb. de 2019
@Jeff - I have already used MLE but some the results were unrealistic/not good so wanted to use method of moments as an alternative.
Lauren Arendse
Lauren Arendse el 20 de Feb. de 2019
@walter thanks for your input. do you know if there is a way to work around this?
also assigned the data as s because I first created the equations and then added the load function. this can be cleaned up though.

Iniciar sesión para comentar.

Respuesta aceptada

Jeff Miller
Jeff Miller el 20 de Feb. de 2019
Lauren,
If you just want to get numerical estimates, Cupid will give them to you. Here is a script for that:
load('Location_01_hs.mat');
xbar = mean(values)
xvar = var(values)
cenmom3 = mean( (values-xbar).^3 ) % Cupid moment estimation uses the 3rd central moment
weib = Weibull(2,2,1); % Just guess some starting parameter values
weib.EstMom([xbar, xvar, cenmom3]) % This is the command for moment estimation
weib.Mean % Check the mean, variance, and 3rd moment with the estimated parameter values.
weib.Variance % These should match the values computed from data.
weib.CenMoment(3)
Here is the output:
xbar =
3.2147
xvar =
2.5325
cenmom3 =
5.3489
ans =
Weibull(2.2442,1.3118,1.1457)
ans =
3.2147
ans =
2.5325
ans =
5.3489
  2 comentarios
Lauren Arendse
Lauren Arendse el 21 de Feb. de 2019
thanks for this. I unfortunately do not have cupid toolbox. however your results are the same as mine for the work around below.
Thanks

Iniciar sesión para comentar.

Más respuestas (1)

Lauren Arendse
Lauren Arendse el 21 de Feb. de 2019
Unfortunately I do not hav the cupid tool box but I found a way that gives me reasonable results. To solve each parameter separately. a=scale , b=shape and c=location parameters. My results for the were the same for @jeff miller's above for the dataset I provide, i.e. a=2.2442, b=1.3118 and c=1.1457
data=zeros(12,3);
for p=1:12;
name=strcat('Location_',num2str(p,'%02d'),'_hs.mat');
load(name);
skew=skewness(values);
syms b positive
eqn=[(gamma(1+3/b)-3*gamma(1+1/b)*gamma(1+2/b)+2*(gamma(1+1/b))^3)/((gamma(1+2/b)-(gamma(1+1/b)^2))^1.5)==skew];
vari=[b];
[solb]=solve(eqn,vari);
b=double(solb);
data(p,2)=b;
clear eqn vari
var=std(values);
syms a positive
eqn=[a*sqrt(gamma(1+2/b)-(gamma(1+1/b))^2)==var];
vari=[a];
[sola]=solve(eqn,vari);
a=double(sola);
data(p,1)=a;
clear eqn vari
u=mean(values);
syms c
eqn=[a*gamma(1+1/b)+c==u];
vari=[c];
[solc]=solve(eqn,vari);
c=double(solc);
data(p,3)=c;
clear eqn vari
end
  4 comentarios
Jan
Jan el 21 de Feb. de 2019
@Lauren: Please explain this in detail. Are you sure that an admin of MathWorks has asked you to set a flag instead of using teh standard methods of voting or accepting the working answer? You have accepted another answer, which marks it as a working solution. Using a flag to post a comment would be a rather tedious option, because there is a number of users and admin who control the flags frequently, because this is the only method to get attention, if somebody is rude or posts illegal contents.
Lauren Arendse
Lauren Arendse el 21 de Feb. de 2019
@jan the textbox with the 'error' is no longer showing, otherwise I would post a screenshot. Basically I recall it something along these lines: 'your answer has been marked as spam, either flag it or email xxxx to correct it' So I flagged it.
I only accepted the previous answer after flaging my answer. Thus the 'error' regarding my answer being spam was showing before I accepted another answer.

Iniciar sesión para comentar.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by