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

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)

[Jeff Miller]

1. Are you just trying to get numerical estimates or do you need to derive some equations (I'm wondering why you are using symbolics)?
2. If you just want estimates, why do you prefer the method of moments over maximum likelihood? Since you have sample data values, you could use ML, which is generally better.

[Walter Roberson]

Your b is symbolic, so when you call g(1+2/b) then you are calling g with a symbolic parameter that cannot be converted to numeric. You defined g as an anonymous function with a single parameter that is to be used as the upper bound for integral(). Because you are passing in an expression that includes an unresolved symbolic variable, you are passing a symbolic upper bound that cannot be converted to numeric to integral(). integral() cannot handle that.

[Walter Roberson]

Why do you assign

s=values;

and then afterwards each place you use s as a variable, it is as an anonyous function parameter ? It is confusing to readers, who could easily confuse the s assigned from values with the "just a parameter" use of s.

[Lauren Arendse]

@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]

@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.

Answer 1

[Jeff Miller]

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

[Lauren Arendse]

thanks for this. I unfortunately do not have cupid toolbox. however your results are the same as mine for the work around below.

Thanks

[Jan]

You can simply download the toolbox from the given link: https://github.com/milleratotago/Cupid

Answer 2

[Lauren Arendse]

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

[Jan]

[MOVED from flags] Lauren Arendse wrote:

this is the solution that worked for me

@Lauren: Please use flags only to inform the admins about inappropriate content like spam or rudeness. Thanks. Use comments, votes or the accepting of an answer to mark it as working.

[Lauren Arendse]

@jan - no problem.

matlab has asked me the flag it because my 'answer' was marked as spam.

[Jan]

@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]

@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.