フィルターのクリア

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

2 ビュー (過去 30 日間)
David
David 2023 年 12 月 26 日
コメント済み: David 2023 年 12 月 26 日
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))]
  5 件のコメント
Torsten
Torsten 2023 年 12 月 26 日
編集済み: Torsten 2023 年 12 月 26 日
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
David 2023 年 12 月 26 日
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.

サインインしてコメントする。

回答 (1 件)

Walter Roberson
Walter Roberson 2023 年 12 月 26 日
編集済み: Walter Roberson 2023 年 12 月 26 日
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.
  1 件のコメント
David
David 2023 年 12 月 26 日
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.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeAssumptions についてさらに検索

製品


リリース

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by