Writing a function in MATLAB

I want to simply write this mathemtical equation in MATLAB, but I am getting some error. Please help.

[Sam Chak]

@Kashif Naukhez,

Would it be possible for you to furnish us with the MATLAB code in question and highlight the specific error that you are encountering? Such information would greatly facilitate our ability to assist you in rectifying the issue at hand.

Click this icon to insert the MATLAB code.

[Kashif Naukhez]

A = readmatrix('Double_q.xlsx');

Error using readmatrix
Unable to find or open 'Double_q.xlsx'. Check the path and filename or file permissions.

x = A(:,1);
y = A(:,2);
loglog(x,y,'r .','MarkerSize', 20)
hold on

q = optimvar('q',2);
beta = optimvar('m',2);
A = optimvar('A',1);
lambda = optimvar('m',1);
x0.q = [1.5,2];
x0.beta = [1,1];
x0.A = 0.5;
x0.lambda = 0.5;

diffun = A*(1 - (1 - q(1))*beta(1)*x).^(1/1-q(1)) + (1 - A)*(1 - (lambda/beta(2)) + (lambda/beta(1))*exp((q(2) - 1)*beta(2)*x)).^(1/1-q(2));

diffexpr = sum((diffun - y).^2);
ssqprob = optimproblem('Objective',diffexpr);
[sol,fval,exitflag,output] = solve(ssqprob,x0);
resp = evaluate(diffun,sol);
hold on
plot(x,resp)
hold off

Error:

Error in .^

Error in Double_q (line 16)

diffun = A*(1 - (1 - q(1))*beta(1)*x).^(1/1-q(1)) + (1 - A)*(1 - (lambda/beta(2)) + (lambda/beta(1))*exp((q(2) - 1)*beta(2)*x)).^(1/1-q(2));

[Kashif Naukhez]

I am trying to fit this equation to my plot. This is what I want.

[Dyuman Joshi]

Please attach the excel file so that we can run your code and reproduce the error you obtained.

In your above comment, did you copy and pasted the full error message you got i.e. all of the red text? If not, then copy and paste the full error message.

Additionally, it's not clear to me what the objective of the optimization is? Is it to minimize the sum you have defined or maximize it?

[Sam Chak]

Hi @Kashif Naukhez

Thank you for sharing the code. I suspect that the error is related to the element-wise power operator (.^) in 'diffun'. To properly test the code, we also require the spreadsheet 'Double_q.xlsx'. Could you please provide that as well?

[Kashif Naukhez]

yes sure!

[Kashif Naukhez]

@Dyuman Joshi I have already attached the excel file. I copied the fulll error message that I got. Objective is to find the constants q(1), q(2) ,lambda, beta(1), beta(2) and A by fitting the equation to the plot.

[dpb]

" Objective is to find the constants q(1), q(2) ,lambda, beta(1), beta(2) and A by fitting the equation to the plot."

For that use lsqnonlin instead...

[dpb]

diffun = A*(1-(1-q(1))*beta(1)*x).^(1/1-q(1)) + (1-A)*(1-(lambda/beta(2))+(lambda/beta(1))*exp((q(2)-1)*beta(2)*x)).^(1/1-q(2));

in the term

(1-(1-q(1))*beta(1)*x).^(1/1-q(1))

(1-(1-q(1))*beta(1)*x) --> (1-1+q(1))*beta(1)*x) --> q(1)*beta(1)*x

or is there a typo in the formula as written?

Then the exponent

.^(1/1-q(1)) --> .^(1-q(1))

Was that intended to be

.^(1/(1-q(1)))

instead?

[Kashif Naukhez]

@dpb Here is the equation. I just need help in writing down the equation in MATLAB.

[Dyuman Joshi]

Yes, it was intended as the last line of code you wrote @dpb, as can be seen in the pdf shared by OP.

However the error still persists after making the correction.

I suspect (probably) because optimization is not the way to go here.

[Kashif Naukhez]

I need to determine the constants given x and y (in the excel file attached). This is my problem.

[Sam Chak]

Hi @Kashif Naukhez,

The x-data in the 'Double_q.xlsx' spreadsheet is not sorted in ascending order; it appears to resemble an exponential decay when sorted in Excel.

Answer 1

[dpb]

% 1.  Associate coefficients with parameter array to estimate
A --> b(1)
q1--> b(2)
q2--> b(3)
beta1>b(4)
beta2>b(5)
lambd>b(6)
% 2. Write anonymous function in b(), x
fun=@(b,x) b(1)*(b(2)*b(4)*x).^(1/(1-b(2))) + (1-b(1))*(1-(b(6)/b(5))+(b(6)/b(5))*exp((b(3)-1)*b(5)*x)).^(1/(1-b(3)));
% 3. Pick initial guess for b0 as estimates
b0=[....];      % look at the data and make some reasonable guesses for b(1) thru b(6)
% 4. Solve for b vector...
b=lsqnonlin(fun,b0);

[Kashif Naukhez]

@Sam Chak Thanks!!