Unable to perform assignment because the size of the left side and right side are different

I want to estimate a 3-parameter weibul distribution. And I got the warning "Consider Preallocating for speed",

so I create the lines with the zeros in the middle.

But for "params = zeros(n,length(b_A),length(T_A));" I get the warning "Unable to perform assignment because the size of the left side is 1-by-5 and the size of the right side is 1-by-3."

Can someone help me plz?

Also if you have any ideas for improving my code, I would appreciate that. I am new on matlab.

n= 10;
t0= 0.5;
b_A= 1:5;
T_A= 1:5;
LowerBound= [0 0 0];
rng('shuffle') 

data= zeros(n,length(b_A),length(T_A));
params= zeros(n,length(b_A),length(T_A));

for k= 1:length(T_A)
for i= 1:length(b_A)
data(:,i,k) = wblrnd(b_A(i),T_A(k), [n,1]) + t0;

start= [b_A(i) T_A(k) t0];
custompdf = @(x,a,b,c) (x>c).*(b/a).*(((x-c)/a).^(b-1)).*exp(-((x-c)/a).^b);
opt = statset('MaxIter',1e5,'MaxFunEvals',1e5,'FunValCheck','off');

params(i,:,k) = mle(data(:,i,k),'pdf',custompdf,'start',start,'Options',opt,'LowerBound',LowerBound,'UpperBound',[Inf Inf min(data(:,i,k))])

end
end

for k= 1:length(T_A)
for i= 1:length(b_A)
    
params(i,4,k) = b_A(i);
params(i,5,k) = T_A(k);
params(i,6,k) = t0;
end
end

Answer 1

[Mara]

Hello Mustafa, when you preallocate the variable params, you have a little twist in there, since you defined it as having the same size as data. However, params has b_A as first, n as second and T_A as third dimension. Later you retrieve data (mle...), it seems like you want to save it as the first 3 elements of the second dim of params. But on the right side you tell MATLAB you want to put it in all 10 columns(i,:,k), which gives you an error (also, you are only specifying the second dimension to element 6, if you don't want the rest to be zeros, you have to replace n by 6).

For the rest I changed some things I could quickly recognize as being redundant:

1. you are having four for-loops, which can be put together. For-loops are slow and from your code I did not see why you made "a second round".
2. you are writing for i = 1:length(T_A), which in the case of a vector of integers from one is the same as writing for i = T_A. Also you do not have to index into T_A/b_A to retrieve the data (T_A(k)) because it is always k or i

clear;
clc;
n   = 10;
t0  = 0.5;
b_A = 1:5;
T_A = 1:5;

LowerBound= [0 0 0];
rng('shuffle');
data    = zeros(n,length(b_A),length(T_A));
params  = zeros(length(b_A),n,length(T_A));
for k= T_A
    for i= b_A
        data(:,i,k) = wblrnd(i,k, [n,1]) + t0;
        start       = [i k t0];
        custompdf   = @(x,a,b,c) (x>c).*(b/a).*(((x-c)/a).^(b-1)).*exp(-((x-c)/a).^b);
        opt         = statset('MaxIter',1e5,'MaxFunEvals',1e5,'FunValCheck','off');
        params(i,1:3,k) = mle(data(:,i,k),'pdf',custompdf,'start',start,'Options',opt,'LowerBound',LowerBound,'UpperBound',[Inf Inf min(data(:,i,k))])
        params(i,4,k)   = i;
        params(i,5,k)   = k;
        params(i,6,k)   = t0;
    end
end

I think it should still work as you want it, but you should check the arrays, better safe than sorry!

P.S I do not know about weibul distribution and whether your code does what it should, I just looked at the code itself.

sry for all the editing

[Mustafa Vural]

Hello Mara,

thank you very much for your answer, it helps me a lot, I really appreciate it.

If you are interested, I tell you what I am doing:

In order to determine the reliability of a product, service life tests are often carried out in practice. Statistical distribution functions can be used to describe the random failure behaviour of components or systems. The unknown parameters of the lifetime distributions are estimated from the existing failure times.

The Weibul distribution is a model with which a distribution function can be described mathematically.

The distribution function can be fully represented by the 2- or 3-parametric Weibull distribution.

To fit the Weibull distribution to data and find parameter estimates you use MLE. The Weibull distribution uses these parameters:

• T is scale parameter ( on the x-axis)
• b is the shape parameter
• t0 is a parameter where the function starts. (Normally, with the 2 parameter distribution, the function starts at the zero point).
• For generating the sample data, we use a size n. Like n=10 or 100 or 1000.

So replacing n with 6 would not work, but you give me a good reference. I just delete n from params, now it works.

I just have one more question. I want to bring the results together, so I can copie them to excel. So I did this at the end.

clear all
n   = 10;
t0  = 0.5;
b_A = 1:5;
T_A = 1:5;
LowerBound= [0 0 0];
rng('shuffle');
data    = zeros(n,length(b_A),length(T_A));
params  = zeros(length(b_A),length(T_A));
for k= T_A
    for i= b_A
        data(:,i,k) = wblrnd(i,k, [n,1]) + t0;
        start       = [i k t0];
        custompdf   = @(x,a,b,c) (x>c).*(b/a).*(((x-c)/a).^(b-1)).*exp(-((x-c)/a).^b);
        opt         = statset('MaxIter',1e5,'MaxFunEvals',1e5,'FunValCheck','off');
        params(i,1:3,k) = mle(data(:,i,k),'pdf',custompdf,'start',start,'Options',opt,'LowerBound',LowerBound,'UpperBound',[Inf Inf min(data(:,i,k))])
        params(i,4,k)   = i;
        params(i,5,k)   = k;
        params(i,6,k)   = t0;
    end
end
result(1:5,:)= params(:,:,1);
result(6:10,:)= params(:,:,2);
result(11:15,:)= params(:,:,3);
result(16:20,:)= params(:,:,4);
result(21:25,:)= params(:,:,5);

I did "result" to copie the results to excel.

When I change b or T from 1:5 to 1:3, I need to change the result manual.

For example: result(1:5,:) to result(1:3,:). How can I make a loop for result?

[Mara]

Hello,

that is good to know!

You can integrate it into your existing loop

clear;
n   = 10;
t0  = 0.5;
b_A = 3;
T_A = 3;
LowerBound= [0 0 0];
rng('shuffle');
data    = zeros(n,b_A, T_A);
params  = zeros(b_A, T_A);
result  = double.empty(b_A*T_A, 0);
for k= 1:T_A
    for i= 1:b_A
        data(:,i,k) = wblrnd(i,k, [n,1]) + t0;
        start       = [i k t0];
        custompdf   = @(x,a,b,c) (x>c).*(b/a).*(((x-c)/a).^(b-1)).*exp(-((x-c)/a).^b);
        opt         = statset('MaxIter',1e5,'MaxFunEvals',1e5,'FunValCheck','off');
        params(i,1:3,k) = mle(data(:,i,k),'pdf',custompdf,'start',start,'Options',opt,'LowerBound',LowerBound,'UpperBound',[Inf Inf min(data(:,i,k))]);
        params(i,4,k)   = i;
        params(i,5,k)   = k;
        params(i,6,k)   = t0;
    end
    result((k-1)*b_A+1:k*b_A, 1:size(params, 2)) = params(:,:,k);
end

It think what you need is clear rather than clear all (here)

Note that I changed the definition of b_A/T_A from vectors to scalars.

Did you notice the warning that pops up for the maximum likelihood estimation?

[Mustafa Vural]

Oh wow,

thank you a lot!

Yes I noticed it. I think this is for some combinations of b and T, but I am not sure. Because when you give b=1:3 and T=1:3, it doesnt appear. I need to look for that.

I just added the 2parameter version. But I think I am not finished with my work. But you helped me a lot!

clear;
n   = 100;
t0  = 0.5;
b_A = 1:5;
T_A = 1:5;
LowerBound= [0 0 0];
rng('shuffle');

data3p    = zeros(n,length(b_A),length(T_A));
params3p  = zeros(length(b_A),length(T_A));

data2P = zeros(n,length(b_A),length(T_A));
params2p = zeros(length(b_A),4,length(T_A));

result3p  = double.empty(length(b_A)*length(T_A), 0);
result2p  = double.empty(length(b_A)*length(T_A), 0);
for k= T_A
    for i= b_A
        data3p(:,i,k) = wblrnd(i,k, [n,1]) + t0;
        data2P(:,i,k) = wblrnd(i,k, [n,1]);
        start       = [i k t0];
        custompdf   = @(x,a,b,c) (x>c).*(b/a).*(((x-c)/a).^(b-1)).*exp(-((x-c)/a).^b);
        opt         = statset('MaxIter',1e5,'MaxFunEvals',1e5,'FunValCheck','off');
        params3p(i,1:3,k) = mle(data3p(:,i,k),'pdf',custompdf,'start',start,'Options',opt,'LowerBound',LowerBound,'UpperBound',[Inf Inf min(data3p(:,i,k))])
        params3p(i,4,k)   = i;
        params3p(i,5,k)   = k;
        params3p(i,6,k)   = t0;
        
        params2p(i,1:2,k) = wblfit(data2P(:,i,k));
        params2p(i,3,k)   = i;
        params2p(i,4,k)   = k;
    end
    result3p((k-1)*length(b_A)+1:k*length(b_A), 1:size(params3p, 2)) = params3p(:,:,k);
    result2p((k-1)*length(b_A)+1:k*length(b_A), 1:size(params2p, 2)) = params2p(:,:,k);
end

[Mustafa Vural]

Hi Mara, thank you again for yor help. It works really good, but I tried to figure out how you do you do it ( I mean the result part). Can you please explain me, how you do the line? I mean how do you knew its (k-1) and why "1:size(params3p), 2)". Can you please explain the line? Its very complicated for me, I dont understand the line. I am very noob in matlab ^^

[Mara]

Hi Mustafa,

sorry for the late answer. At that point I am defining where to put the new data. We need to define how many and which rows and columns of the array result3p will be used. For the new data we need 5 rows and 6 columns and we need to put it behind the data that is already there from the previous loops.

So (k-1)*b_A calculates how many rows have already been filled up in the rounds before. For example if you run the k-loop for the third time, the first new line must be Nr. 11. Why? Because Loop k=1 used 5 rows and k=2 used another 5 rows. So I for k = 3 it is calculating (3-1)*5+1=11

the size of the second dimension of params3p is 6. So the code 1:size(params3p, 2) will be the vector 1, 2, 3, 4, 5, 6. These will be the columns for the data in result3p.