Error at finding minimum

Hi, I try to find the minimum of a sum. It does work for a fixed parameter for 'test', but I want to show it for various values of 'test'. I always get the error saying: "Index exceeds array bounds" in the line where I define 'c_hat = single(c_hat(1))', what is wrong in my code? (In the end I want to plot the results for each value of test(j), I deleted this in the code shared to make it more simple.) I would be very thankful if you can help me!

lambda = 0.95; %close to 1
n = 100; %number of samples
theta = 0.03;
alph = 2.5; %divergence level
%Generate Pareto distributed r.v.:
x_m = 1; %scale parameter for Pareto
alpha = 3; %shape parameter for Pareto
pd = makedist('uniform');
u = random(pd,[1,n]);
X = zeros(1,n);
for i=1:n
  X(i) = x_m / (1-u(i))^(1/alpha); %Pareto distributed r.v.
end
a_1 = x_m / (1-lambda).^(1/alpha)
V = zeros(n,1);
a=a_1;
for i=1:n
  V(i) = 1/(1-lambda) .* max(X(i)-a,0) + a;
end
k=20;
test=linspace(0,10,k);
E=zeros(k,1);
rel_entropy=zeros(k,1);
alph_div = zeros(k,1);
a_hat = zeros(k,1);
for j=1:k
  syms c
  summe_5=0;
  for i=1:n
  summe_5 =  summe_5 + (test(j)*(alph -1)*V(i) + c)^(1/(alph-1));
  end 
  cond = (1/n)* summe_5 == 1;
  c_hat = zeros(3,1);
  for q = 1:3 
  c_hat = vpasolve(cond,c, 'random', true)  ;
  end
  c_hat = single(c_hat(1))
  %syms b
  size(c_hat)
  M_b = @(b) ((alph-1)/alph) * sum((test(j).*(alph-1) .* (1/(1-lambda) .* max(X-b,0)+b) + c_hat).^(1/(alph-1)).*(1/(1-lambda) .* max(X-b,0)+b)) + (c_hat/(test(j) .* alph * (1-alph)));
  a_hat(j) = fminsearch(M_b,a_1);
  V_a = zeros(n,1);
  for i=1:n
      V_a(i) = 1/(1-lambda) .* max(X(i)-a_hat(j),0) + a_hat(j);
  end
  m_hat = zeros(n,1);
  for i=1:n
      m_hat(i) = (test(j) * (alph-1) * V_a(i)+ c_hat) .^(1/(alph-1));
  end
  E(j)= (1/n) * sum(V.*m_hat);
  %plotting results 
end

Answer 1

[Walter Roberson]

You would get that if vpasolve is returning empty.

What is the purpose of your q loop in which you vpasolve three times but do not permanently save the first two results?

[Verena Schnieders]

Oh I needed that for another file, I checked whether all results are the same. I do not necessarily need to solve it several times in this case, sorry. But it is still the same error, even without the loop and only computing c_hat once.

[Walter Roberson]

I confirm that when j = 5, that the generated cond has no solutions over real values -- or at least that was the case with whatever random number seed happened to get used when I tried.

You generate an expression that looks like

(1/100)*(c+5247771914842387/35184372088832)^(2/3)+97*(c+301594987663439/35184372088832)^(2/3)*(1/100)+(1/100)*(c+6102129260550525/140737488355328)^(2/3)+(1/100)*(c+290574697259709/549755813888)^(2/3) == 1

Take the minimum of the values that are being added to c, and take the negative of that value. If c goes lower than that value then you would be raising a negative value to 2/3, which would give a complex result. So by looking at your

(test(j)*(alph -1)*V(i) + c)^(1/(alph-1)

terms and breaking that down to

 test(j)*(alph -1)*V(i)

and finding the minimum of those, and taking the negative, you get a lower bound on where to search. If the value of

(1/n)* summe_5

is greater than 1 at that point, then there cannot be any real-valued solutions (because the other terms will all be greater value.)

[Verena Schnieders]

Thank you! I fixed the error, I only choose test to be in the interval of (0,1), then the sum is alright and everything works out fine :-)