Need to solve the following equation with three knowns and 2 unknowns.

W0 = Wi + (1 - L)*Hi*gi

W0 and L is an unknowns.

Wi, Hi, and gi are the knowns with 172 values per each.

I want to find W0 and L and their standard deviations as well.

I want to get each values comes for W0 and L when the program runs.

But when this code runs i got 0 for standard deviation of both W10 and L. Is this code is right? Can anyone help me to solve this?

The following code was written.

Wi = Wi;

Unrecognized function or variable 'Wi'.

Hi = H_BM;  
gi = g_mean;  

% Define the objective function
fun = @(x) norm(Wi + (1 - x(1))*Hi.*gi - x(2));

% Set initial guess for L and W0
x0 = [0, 0];

% Solve the optimization problem
x = fminsearch(fun, x0);

% Extract the values of L and W0
L = x(1);
L
W0 = x(2);
W0

% Calculate standard deviation of W0 and L
std_W0 = std(W0);
std_W0
std_L = std(L);
std_L

[Rik]

You are calculating the standard deviation of a single value, not of the fit. fminsearch does not return the goodness of fit parameters with this syntax.

[Prabhath Manuranga]

Could you please tell me how to fix it?

[Rik]

Since you didn't post any data I can't show you what to do. You might want to try the fit function, since that will report the goodness of fit as the second output argument.

[Prabhath Manuranga]

Please find the attached excel file for Wi, H_BM, and g_mean values.

Answer 1

[Torsten]

H_BM = rand(100,1);
g_mean = rand(100,1);
Wi = rand(100,1);
Hi = H_BM;  
gi = g_mean;  
x = Hi.*gi;
y = Wi + Hi.*gi;
fitlm(x,y)

ans = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate       SE       tStat       pValue  
                   ________    ________    ______    __________

    (Intercept)    0.41233     0.043102    9.5664    1.0575e-15
    x1              1.2316      0.15371     8.013    2.3796e-12


Number of observations: 100, Error degrees of freedom: 98
Root Mean Squared Error: 0.294
R-squared: 0.396,  Adjusted R-Squared: 0.39
F-statistic vs. constant model: 64.2, p-value = 2.38e-12

[Prabhath Manuranga]

Dear Mr. Torsten,

y= Wi + Hi.*gi; or y = Wi + (1-L)*Hi.*gi;

how to get two unknowns W0 and L, and their standard deviations here?

[Torsten]

W0 = Estimate intercept

L = Estimate x1

SE = standard errors of the estimated parameters

For further details, consult the documentation of "fitlm".

[Prabhath Manuranga]

Dear Mr. Torsten,

For this equation can i use the same method you mentioned above. Above equation there is no any Pi (because Pi = 1). but here i need to use Pi weight for the calculation.

[Torsten]

You have a linear regression equation

W0 + (L - 1)*Hi*gi = Wi

because if you set

W0 = a, L - 1 = b, Hi*gi = xi, Wi = yi

, your equation reads

a + b*x = y

You want to estimate a and b to make

sum_{i=1}^{i=172} (yi-(a+b*xi))^2

minimal.

This can be done via "fitlm" as I posted.

There are easier methods to get the coefficient a and b, but you wrote you also want statistical information about reliability etc.

[Torsten]

format long
M = readmatrix("data123.xlsx")

M = 10x3

1.0e+07 * 6.263682480329910 0.000000205900000 0.097810856393885 6.263657062323940 0.000002802300000 0.097812165123015 6.263669605757160 0.000001524500000 0.097811320851740 6.263624373910401 0.000006129800000 0.097810680733666 6.263620997894920 0.000006478400000 0.097810156445693 6.263650330102171 0.000003480200000 0.097811560786981 6.263682602864220 0.000000192100000 0.097810613817760 6.263680818837160 0.000000373000000 0.097808759764806 6.263617228951360 0.000006862700000 0.097809507129465 6.263661105500140 0.000002381500000 0.097808297704551

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Wi = M(:,1);
H_BM = M(:,2);
g_mean = M(:,3);
Hi = H_BM;  
gi = g_mean;  
x = Hi.*gi;
y = Wi;
fitlm(x,y)

ans = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                         Estimate           SE    tStat    pValue
                   _____________________    __    _____    ______

    (Intercept)         62636844.9371338    0      Inf       0   
    x1             -1.00231759363032e-05    0     -Inf       0   


Number of observations: 10, Error degrees of freedom: 8
R-squared: 1,  Adjusted R-Squared: 1
F-statistic vs. constant model: 4.65e+06, p-value = 2.4e-24

n = numel(Wi);
Mat = [ones(n,1),x];
rhs = y;
sol = Mat\rhs;
W0 = sol(1)

W0 = 
     6.263684493713374e+07

L = 1 + sol(2)

L = 
   0.999989976824064

[Prabhath Manuranga]

Dear Mr. Torsten, Could you please tell me how to set this using fitlm function? is the following code is ok to calculate W0 and L?

% Create a table or dataset with the known variables

data = table(Wi, Hi, gi, Pi);

% Calculate Pi

Pi = 1./Hi;

% Calculate the numerator and denominator sums

numerator = sum(Pi .* (Wi + (1 - L) .* Hi .* gi));

denominator = sum(Pi);

% Calculate W0

W0 = numerator / denominator;

% Fit a linear model to the data

model = fitlm(data, 'Hi ~ (W0 - Wi)./gi + L*Hi');

% Get the estimated coefficients

W0 = model.Coefficients.Estimate(1);

L = model.Coefficients.Estimate(2);

[Torsten]

As far as I can see, W0 = numerator/denominator. Thus you already computed W0. What do you want to fit here ?

As far as I understood your problem, you have values for Wi, Hi and Gi and you want to compute L and W0 to minimize the sum of squared differences

sum_{i=1}^{i=172} (W0 + (L - 1)*Hi*gi - Wi)^2

That's what fitlm does if you call it as I did.

[Prabhath Manuranga]

The question is i could not understand where i should place the Pi value in the equation you mentioned here. Doesn't it need for the calculation? if i get sum_{i=1}^{i=172} Pi = K, the equation will be

K*W0 = sum_{i=1}^{i=172} (1/Hi).* (Wi - (1-L).*Hi.*gi).

[Torsten]

I don't understand how you come to the equation for W0, especially because you state that L on the right-hand side is also unknown.

Thus

sum(Pi .* (Wi + (1 - L) .* Hi .* gi))/sum(Pi)

cannot be used to determine W0.

In order to use fitlm, you must specify the independent variable x and the independent variable y in the regression model y = a*x + b. Given (in your case 172) realizations of x and y, fitlm will then try to adjust a and b.

Maybe the original regression equation was

P*W0 = P*W + (1-L)*g

?

[Prabhath Manuranga]

The problem i faced can be described as follows.

By using Eq 1 i found the valyue for Eq2 (W0 hat).

Then using Eq 4 ei was found.

then the problem start from Eq 5. because as i mentioned earlier W0 and L are unknowns. only Wi, gi, and Hi are the knowns and Eq6 Wo hat to be found.

As you mentioned earlier using your fitlm equation value W0 found. these equations were copied from a published research paper. they said W0 hat to be found using Equaton 6 and lease square adjustment was applied for Pi = 1/Hi, Pi = 1/Hi^2, and Pi = 1/Hi^0.5.

Could you please assist me to solve this problem sir?

[Torsten]

If this is really what is written in the article, I don't understand what the authors did.

As written, I would assume that L was treated as known because only a value W0^ was computed.

[Prabhath Manuranga]

No, they mentioned both W0 hat and L(lamda) were calculated using the Eq 5.

[Torsten]

If L and W0 were computed using equation (5) with a least-squares approach, then the two ways to do so are given in my answer where I used your data.

As said: I can't understand what the authors of the article did. E.g. I cannot find any reason why the authors multiplied equation (1) by Pi before summing it from i=1 to i=m and finally dividing by sum_i_Pi. Maybe they gave weights to the measurement data this way, but this should be mentionned in the article.

[Prabhath Manuranga]

Yes you are right. They did use Pi to weight the measurements (as discussed with my supervisor he also agreed with your idea). when applying different Pi values as i mentioned earlier (Pi = 1, Pi = 1/Hi, Pi = 1/Hi^2, and Pi = 1/Hi^0.5) they calculated four different W0 values. Is it possible to find four different W0 values with their L as you mentioned earlier.

here sum4 means summation of Hi values ({i=1}^{i=172} Pi = sum4).

sum4*W0 = (1/Hi)*(Wi + (1-L)*Hi*gi)

W0 = (Wi/Hi)*(1/sum4) + (1-L)*gi/sum4

i tried to get similar y =a + bx linear equation to make fit function. but did not get the correct answer for W0 and L.

x = gi;

y = (1/sum4).*(Wi./Hi) + gi;

fitlm(x,y)

[Torsten]

Note that the weights for "fitlm" must be the usual weights squared.

format long
M = readmatrix("data123.xlsx");
Wi = M(:,1);
H_BM = M(:,2);
g_mean = M(:,3);
Hi = H_BM;  
gi = g_mean;  
x = Hi.*gi;
y = Wi;
weights = 1./Hi / sum(1./Hi);
s = fitlm(x,y,'Weights',weights.^2)

s = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                         Estimate                    SE                   tStat                 pValue       
                   _____________________    ____________________    _________________    ____________________

    (Intercept)         62636844.9002729       0.535047793115774     117067756.761534    3.17482396862639e-62
    x1             -1.00215505545489e-05    1.32603813010255e-07    -75.5751311146224    1.04712464605067e-12


Number of observations: 10, Error degrees of freedom: 8
Root Mean Squared Error: 0.217
R-squared: 1,  Adjusted R-Squared: 1
F-statistic vs. constant model: 7.87e+05, p-value = 2.92e-21

Wohat = s.Coefficients.Estimate(1)

Wohat = 
     6.263684490027287e+07

Lhat = 1 + s.Coefficients.Estimate(2)

Lhat = 
   0.999989978449445

n = numel(Wi);
Mat = [weights,weights.*x];
rhs = weights.*y;
sol = Mat\rhs;
W0 = sol(1)

W0 = 
     6.263684490027288e+07

L = 1 + sol(2)

L = 
   0.999989978449443

syms W0 L
f = sum((weights.*(W0-(Wi + (1 - L)*Hi.*gi))).^2);
dfdW0 = diff(f,W0);
dfdL = diff(f,L);
sol = solve([dfdW0==0,dfdL==0],[W0,L])

sol = struct with fields:
    W0: 12040498260855305744269849528779227529092954537156824557392557339465700185329217840862105864303847686272978470720568140168141/192227087427943729367066220859610660420276...
     L: 57287562195798430748523865830601191935769536305576007891889226861407737951969429202951436538059542930491753/572881363117527155181962909876139940083851623062070447252942...

W0 = vpa(sol.W0)

W0 = 

62636844.90027287747361312966808

L = vpa(sol.L)

L = 

0.99998997844944439239438665630244

[Prabhath Manuranga]

Thanks Mr Torsten. I greatly appriciate your support reagrding my issue.

sir, when Pi = 1/Hi^2, can it be the weight like as follows,

weights = 1./Hi.^2 / sum(1./Hi.^2);

s = fitlm(x,y,'Weights',weights.^2)

when Pi = 1/Hi^0.5, can it be the weight like as follows,

weights = 1./Hi.^0.5 / sum(1./Hi.^0.5);

s = fitlm(x,y,'Weights',weights.^2)

[Torsten]

Yes.

[Prabhath Manuranga]

The following code was written for using above two equations. Here n = 2159 and n* = 78. Could you please let me know any errors in my written code sir? I could not attach EGM08 file, because the size of the file is 228MB.

clear, clc

format longg

GM = 0.3986004415E+15;

alpha = 0.63781363E+07;

% Load data

data = readmatrix('EGM08.txt');

% Extract sigma and beta values

sigma = data(:,5); % Represent sigma_Cnm of the data set

beta = data(:,6); % Represent sigma_Snm of the data set

sigma1 = sigma(end-2190:end); % Access only last 2190 data set

beta1 = beta(end-2190:end); % Access only last 2190 data set

sum = 0;

for m = 0:2159

sum(m+1) = (GM/alpha).^2*(sigma1(m+1).^2 + beta1(m+1).^2);

% disp(sum);

end

sum = sum.';

Dw = 0;

for n = 2:78

Dw = sqrt(Dw + sum(n));

end

Dw

[Torsten]

Rename your variable "sum". "sum" is an inbuilt MATLAB function to sum elements, and you overwrite this function in your code.

If you are sure that sigma_e^2(V_n) equals what you compute as sum(n) (I doubt it because you don't sum anything when computing sum(n)), your Dw should be computed as

Dw = sqrt(sum(s(2:78)))

Here, I renamed you variable "sum" to "s" and used the MATLAB function "sum" to sum elements in an array.

According to the mathematical formulae, the sigma1 and beta1 are lower triangular matrices, and your "sum" variable is computed by summing both sigma1 and beta1 over their columns, adding the result and multiply it by (GM/alpha)^2.

[Prabhath Manuranga]

Thank you sir.

Answer 2

[Rupesh]

Hi Prabhat,

I understand that you want to calculate the values of W0 and L, and their standard deviations, from your given data using an optimization approach. The initial code calculates W0 and L but results in zero standard deviations because it operates on single scalar values instead of distributions (set of values) and single scalar values don’t have standard deviation associated with them.

To solve this, we can use a bootstrapping method. By repeatedly sampling your data and performing the optimization for each sample, we can obtain distributions of W0 and L. From these distributions, we can calculate the standard deviations.

Wi = Wi;  % Your known values 
Hi = H_BM;  % Your known values 
gi = g_mean;  % Your known values 

% Number of bootstrap samples 
num_samples = 100; 

% Arrays to store the results 
L_values = zeros(1, num_samples); 
W0_values = zeros(1, num_samples); 

% Perform bootstrap sampling 
for i = 1:num_samples 
    % Randomly sample indices with replacement 
    sample_indices = randsample(length(Wi), length(Wi), true); 
    
    % Sample the data 
    Wi_sample = Wi(sample_indices); 
    Hi_sample = Hi(sample_indices); 
    gi_sample = gi(sample_indices); 
    
    % Define the objective function for this sample 
    fun = @(x) norm(Wi_sample + (1 - x(1))*Hi_sample.*gi_sample - x(2)); 
    
    % Set initial guess for L and W0 
    x0 = [0, 0]; 
    
    % Solve the optimization problem for this sample 
    x = fminsearch(fun, x0); 
    
    % Store the results 
    L_values(i) = x(1); 
    W0_values(i) = x(2); 
end 

% Calculate mean and standard deviation of L and W0 
L_mean = mean(L_values); 
L_std = std(L_values); 
W0_mean = mean(W0_values); 
W0_std = std(W0_values); 

% Display the results 
disp(['L mean: ', num2str(L_mean)]); 
disp(['L std: ', num2str(L_std)]); 
disp(['W0 mean: ', num2str(W0_mean)]); 
disp(['W0 std: ', num2str(W0_std)]); 

This script repeatedly samples your data and runs the optimization to build up distributions for W0 and L, allowing for the calculation of their standard deviations. You can refer to below documentation to get clear understanding of inbuilt sample bootstrapping.

https://in.mathworks.com/help/stats/bootstrp.html

Hope it helps!

Thanks