Matlab density function evaluation using symprod

For a homework question we are required to write a matlab function to evaluate the above density function, Im having trouble thinking how to evaluate a function like this with series sums and products in it, below is my attempt about going about solving it, however it appears you can not lookup vector values in symprod?, would it be easier to use a for loop approach? I would appreciate some help on getting started on this one as Im very stuck.

%{ 
(a) Write a Matlab function

c = Clayton c(u, theta)

which returns the value of c(u)

with parameter θ=theta.

%} 

 function c = Clayton_c(u,theta) 
    d=length(u); 
    %Power of second braceted expression(Nb: No pwr1) 
    pwr2=(-1-theta); 
    %Power of the third bracketedexpression(Nb: No pwr1) 
    pwr3=(-d-(1/theta)); 

    syms k
    part1=symprod(1+theta*(k-1),k,1,d);  
    part2=(symprod(u(k),k,1,d).^pwr2);
    part3=(1-d+(symsum((u(k)).^(-theta)),k,1,d).^pwr3); 

    c=part1.*part2.*part3; 

 end

[David Goodmanson]

Hi Tom, Is it part of the problem to use sym variables, or is a straight numerical evaluation all right? If it's the latter, you can keep u as a vector, build things up in the first part using the vector (0:d-1), do term-by-term exponentiation of u in the third part as in u.^-theta, and use the sum and prod commands on vectors where needed.

[Tom Craven]

Hi David, a straight numerical evaluation is fine, would you possibly be able to elaborate on how to do the first part ? at all evaluating from 0:d-1?.

[David Goodmanson]

Hi Tom, I wish I had looked at this thread earlier. Maybe there is a way to be notified when someone adds a comment but if so I don't know what it is. Anyway, here is the idea.

j = 1:d;                      % vector of values of j
A = 1 + theta*(j-1);          % vector of terms you need
first_part = prod(A);         % product of all those

This is probably clearer than the original suggestion which would have gone j = 0:d-1; A = 1 + theta*j;

Answer 1

[Walter Roberson]

You are correct, you cannot look up vector values using the symbolic index using symsum or symprod . Instead you need to generate the entire vector and sum() or prod() it. So instead of

symsum( u(j).^(-theta), j, 1, d)

you would

sum( u(1:d) .^ (-theta) )

and since d is length of u, that could be simplified to

sum( u .^ (-theta) )