construct a local function using a mixture of real and symbolic math

Would love to have help with two things. This local function runs, but …

• I’m getting partially symbolic output and I need real-valued output. (See output sample at bottom.)
• How to make this more efficient. This version runs in about 20 seconds, but I really want to set finer grained values in the loops; i.e., generating 100 times the current number of function calls and that will burn some time in the present form. Is it possible to turn symsum into a @ type function. I tried and failed.

Thanks so much if you can help on this. This is my first venture with the symbolic toolbox, so my guess is that I am really screwing this up and missing important pieces.

Dave C.

Here’s the function “symbolically”; G = gamma function

Sum(kk=0 to Inf) { (rho*sqrt(ab)/(1-rho^2) * [G((kk+1)/2]/[kk! G((kk+m)/2)] }

% Contours_ChiSq.m
% In progress
% ISum is the local function that needs work.
clear all; close all; clc;

rho = 0.5;
k = 1;
for a = 0.1:0.05:1
    for b = (1-a):0.05:1
        Xab(k,:) = [a b ISum(rho,a,b)];
        k = k+1;
    end
end
Xab;


%--------------------------------------------------------------------------
% Call the LOCAL infinite sum function below
%--------------------------------------------------------------------------
sympref('FloatingPointOutput',true);
function [s] = ISum(rho,a,b)
m = 4;
kk =sym('kk');
S = (rho*sqrt(a*b)/(1-rho))^kk;
Ng = gamma((kk+1)/2);
Dg = factorial(kk)*gamma((kk+m)/2);
ss = vpa(symsum(S,kk,0,inf),5);
s = ss*vpa(Ng/Dg,5);
end

Sample output:

Xab =

[0.1000, 0.9000, (1.4286*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]

[0.1000, 0.9500, (1.4455*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]

[0.1000, 1, (1.4625*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]

[0.1500, 0.8500, (1.5554*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]

[Dyuman Joshi]

The output isn't partially symbolic, it's entirely symbolic. If you want a numeric output, you will have to provide a value for kk.

However, you can try simplifying the equation using simplify. The reason I say try is because there is no guarantee that there exists a (further) simplification of a given expression.

Using vpasum instead of symsum appears to be much faster. And you should preallocate the output variable.

rho = 0.5;
k = 1;
sympref('FloatingPointOutput',true);
digits(5)


A = 0.1:0.05:1;
B = (1-A):0.05:1;

Xab = zeros(numel(A)*numel(B), 3, 'sym');
tic
for a = A
    for b = B
        Xab(k,:) = [a b ISum(rho,a,b)];
        k = k+1;
    end
end
toc

Elapsed time is 1.095712 seconds.

out = simplify(Xab)

out = 

function [s] = ISum(rho,a,b)
m = 4;
kk = sym('kk');
S = (rho*sqrt(a*b)/(1-rho))^kk;
ND = gamma((kk+1)/2)/factorial(kk)*gamma((kk+m)/2);
ss = vpasum(S,kk,0,inf);
s = ss*ND;
end

[Walter Roberson]

Note that it is not the symsum() that is leaving you with the kk . ND is constructed in terms of kk so when you multiply the symsum result by ND you get a result that is in terms of kk.

[David]

Dyuman,

Yes, that kk thing is baffling me. It's not clear to me how to instantiate it because on the one hand it is the incrementor in the infinite sum and at the same time it is a variable in the closed-form part of the overall expression. The original expression clearly shows that the closed-form part (with the gamma function) is not part of the infinite sum. (Below) Any further advice would be great. (Except for u,v (or a,b in my version) the first part of this expression is easy to evaluate. It's the part from the infinite sum on that is giving me the headache.

[Torsten]

S = (rho*sqrt(a*b)/(1-rho^2))^kk;

instead of

S = (rho*sqrt(a*b)/(1-rho))^kk;

and

ss = vpa(symsum(S*Ng/Dg,kk,0,inf),5);

instead of

ss = vpa(symsum(S,kk,0,inf),5);

(Ng and Dg depend on kk - they cannot be taken out of the sum).

Why do you work with symbolic variables at all ? Sum with doubles until the size of the summands get smaller than a specified tolerance.

[David]

Torsten,

First, thanks for the vigilent catch of the square on rho.

Second, your suggestion about avoiding symbolic variables altogether is great. My first instinct was to see if that series (which has you've noted must include the gamma factor ... even though the original mathematical expression does not seem to do so) has a finite expression that I could then use. But, just using a finite approximation and avoiding all the symbolic stuff ... wonderful.

Thanks very much.

Answer 1

[Walter Roberson]

You can simplify a lot.

rho = sym(0.5);
syms a b
assume(a > 0 & a <= 1);
assumeAlso(b > 0 & b <= 1);
m = 4;
kk =sym('kk');
S = (rho*sqrt(a*b)/(1-rho))^kk;
Ng = gamma((kk+1)/2);
Dg = factorial(kk)*gamma((kk+m)/2);
SS = symsum(S,kk,0,inf)

SS = 

The second condition will not occur.

So your ss will be infinity if a = 1 and b = 1, and will be 1/(1 - sqrt(a)*sqrt(b)) otherwise. For practical purposes, you can reduce that to just 1/(1 - sqrt(a)*sqrt(b)) as that will come out as infinity when a = 1 and b = 1.

[David]

Walter,

I appreciate your insights here. Let me be sure I understand what you are saying. Using the symbolic math toolbox, you've deduced that the infinite sum has the finite solution 1/(1 - sqrt(a)*sqrt(b)) when a,b are in the windows specified. Is this correct. In other words, I can replace the SS = symsum ( S,kk,),inf) with the 1/(1 - sqrt(a)*sqrt(b)).

Let me know if my interpretation is correct so far. I do appreciate your time and attention.