How to define a function with multiple handle in a loop?
显示 更早的评论
I have a situation where I need to use a function that can take different number of input (depnding upon the number of treatment considered).
The function has a general expression (Attached as Pic for readability)
Currently I am defining the function (for different number of treatment) seperately (as bellow)
%% Parameters
par.a=80/365;% %; %bitting frequiency (per day)
par.g=1/10; %mosquito death rate (per day)
par.n=12; %duration of sporogony in mopsquito
par.b=.5; %transmission probability: mosquito to human
par.c=.23; %transmission probability: human to mosquito
par.r=1/60; %rate of blood stage infection clearance
par.omega=1/425 ; %hypnozoites death rate
par.nu=8.5; %Number of hypnozoites per infection
par.alpha=1/332;%1/332;
p_blood=.9; %treatment effect
p_rad=.9; %treatment effect
% tau is the current time
% s1,s2,s3 are treatment time.
pA=@(tau) par.alpha*(exp(-par.r*tau)-exp(-(par.alpha+par.omega)*tau))/(par.alpha+par.omega-par.r); %Base function without treatment
pA_rad=@(tau,s1) (1-p_blood)*exp(-par.r*(tau-s1)).*pA(s1)+(1-p_rad)*(pA(tau)-exp(-par.r*(tau-s1)).*pA(s1)); %with one treatment
pA_rad2=@(tau,s1,s2) (1-p_blood)*exp(-par.r*(tau-s2)).*pA_rad(s2,s1)+(1-p_rad)^2*(pA(tau)-exp(-par.r*(tau-s2)).*pA(s2)); %with two treatment
pA_rad3=@(tau,s1,s2,s3) (1-p_blood)*exp(-par.r*(tau-s3)).*pA_rad2(s3,s1,s2)+(1-p_rad)^3*(pA(tau)-exp(-par.r*(tau-s3)).*pA(s3)); %with three
as I couldn't figure out how to define these in a single function. The main issue is that as the number of treatment increases, the function handle increases as well and the function with larger handles utilise the function with lower handles for different inputs (as is attached pic). This is making me crazy. Is there any way that I can define these in a loop or nested function? Any help is much appreciated.
采纳的回答
Jan
2022-9-11
Is there a need to use an anonymous function? Why not writing a standard function?
function y = pA(tau, s1, s2, s3)
a=80/365;% %; %bitting frequiency (per day)
g=1/10; %mosquito death rate (per day)
n=12; %duration of sporogony in mopsquito
b=.5; %transmission probability: mosquito to human
c=.23; %transmission probability: human to mosquito
r=1/60; %rate of blood stage infection clearance
omega=1/425 ; %hypnozoites death rate
nu=8.5; %Number of hypnozoites per infection
alpha=1/332;%1/332;
p_blood=.9; %treatment effect
p_rad=.9; %treatment effect
% tau is the current time
% s1,s2,s3 are treatment time.
y0 = alpha*(exp(-r*tau)-exp(-(alpha+omega)*tau))/(alpha+omega-par.r);
switch nargin
case 1
y = y0; %Base function without treatment
case 2
y = (1-p_blood)*exp(-r*(tau-s1)).*pA(s1)+(1-p_rad)*(y0-exp(-r*(tau-s1)).*pA(s1)); %with one treatment
case 3
y = (1-p_blood)*exp(-r*(tau-s2)).*pA(s2,s1)+(1-p_rad)^2*(y0-exp(-r*(tau-s2)).*pA(s2)); %with two treatment
case 4
y = (1-p_blood)*exp(-r*(tau-s3)).*pA(s3,s1,s2)+(1-p_rad)^3*(y0-exp(-r*(tau-s3)).*pA(s3)); %with three
otherwise
error('Too many inputs');
end
end
8 个评论
Md Nurul Anwar
2022-9-11
编辑：Md Nurul Anwar
2022-9-12
Jan
2022-9-12
For the call pA(100, par) the 2nd argument is copied to the variable s1 for this function header:
function y = pA(tau,s1,par)
and the 3rd input is undefined.
Better:
function y = pA(tau, varargin)
if nargin < 2
error('need parameters')
end
par = varargin{end};
y0 = par.alpha*(exp(-par.gamma*tau)-exp(-(par.alpha+par.mu)*tau))/(par.alpha+par.mu-par.gamma);
switch nargin
case 2
y = y0;
case 3
s1 = varargin{1};
y = (1-par.p_blood)*exp(-par.gamma*(tau-s1)).*pA(s1,par)+(1-par.p_rad)*(y0-exp(-par.gamma*(tau-s1)).*pA(s1,par)); %with one treatment
case 4
s1 = varargin{1};
s2 = varargin{2};
y = (1-par.p_blood)*exp(-par.gamma*(tau-s2)).*pA(s2,s1,par)+(1-par.p_rad)^2*(y0-exp(-par.gamma*(tau-s2)).*pA(s2,par)); %with two treatment
case 5
s1 = varargin{1};
s2 = varargin{2};
s3 = varargin{3};
y = (1-par.p_blood)*exp(-par.gamma*(tau-s3)).*pA(s3,s1,s2,par)+(1-par.p_rad)^3*(y0-exp(-par.gamma*(tau-s3)).*pA(s3,par)); %with three
otherwise
error('Too many inputs');
end
end
Cleaner:
function y = pA(tau, par, s1, s2, s3)
if nargin < 2
error('need parameters')
end
par = varargin{end};
y0 = par.alpha*(exp(-par.gamma*tau)-exp(-(par.alpha+par.mu)*tau))/(par.alpha+par.mu-par.gamma);
switch nargin
case 2
y = y0;
case 3
y = (1-par.p_blood)*exp(-par.gamma*(tau-s1)).*pA(s1,par)+(1-par.p_rad)*(y0-exp(-par.gamma*(tau-s1)).*pA(s1,par)); %with one treatment
... same as above without using varargin.
It looks strange, that the order of arguments is:
pA(s1, par)
pA(s2, s1, par)
pA(s3, s1, s2, par)
If you get 20 terms, using a loop would be smarter. What is the aczual formula?
Md Nurul Anwar
2022-9-13
编辑：Md Nurul Anwar
2022-9-13
If you provide par as 2nd input, there is no need to call varargin anymore. Omit line 6:
par = varargin{end}; % Not needed here
But then the other calls do not match anymore: pA(s2,s1,par) and pA(s3,s1,s2,par).
Jan
2022-9-13
This is a recursive impelemtation of your formula:
function y = pA(t, p, varargin)
N = nargin - 2;
if N < 0
error('need parameters')
end
if N == 0 % No s as inputs:
y = p.alpha / (p.alpha + p.mu - p.gamma) * ...
(exp(-p.gamma * t) - exp(-(p.alpha + p.mu) * t)) ;
else % s1..sN are provided as inputs - recursive call:
sN = varargin{N};
y = (1 - p.p_blood) * exp(-p.gamma * (t - sN)) * ...
pA(sN, p, varargin{1:N-1});
end
end
I've abbreviated "tau" to "t" and "par" to "p", because it improves the readability.
A recursive function can be converted to a loop in every case. But at least this implementation is much cleaner than a pile of nested anonymous functions.
During editing I struggled with the naming conventions: "par.alpha + par.omega - par.r". In the formula the "omega" is a mu and "r" is a gamma. The "t" is not a greek "tau".
I do no see, that your anonymous function for s1 matchs the formula:
pA_1 = @(tau,s1) (1 - p.p_blood) * exp(-p.gamma * (tau-s1)) * pA_0(s1) + ...
(1-p.p_rad)*(pA_0(tau) - exp(-p.gamma*(tau-s1)).*pA_0(s1)); %with one treatment
While the first line is equivalent to 
the 2nd line is not clear to me. p_rad does not appear in the formula.
Md Nurul Anwar
2022-9-14
移动：Jan
2022-9-14
Appology for the confusion. I used pA_rad for the subscript 'r'. The complete formula is below ( as in code)
.
.The 2nd line as you mentioned is the y0 (here pA :no treatment)
Thanks a lot for these idea so far.
Jan
2022-9-14
@Md Nurul Anwar: Please care for posting the correct formula directly in the question the next time.
function y = pA(t, p, varargin)
N = nargin - 2;
if N == 0 % No s as inputs:
y = p.alpha / (p.alpha + p.mu - p.gamma) * ...
(exp(-p.gamma * t) - exp(-(p.alpha + p.mu) * t)) ;
elseif N > 0 % s1..sN are provided as inputs - recursive call:
sN = varargin{N};
y = (1 - p.p_blood) * exp(-p.gamma * (t - sN)) * ...
pA(sN, p, varargin{1:N-1}) + ...
(1 - p.p_rad)^N * (pA(t, p) - exp(-p.gamma * (t - sN)) * pA(sN, p));
else
error('need parameters')
end
end
Please test this, if it produces the correct answers. It works with an arbitrary number of inputs s1, ...sN limited by the recursion depth of 500. If this is a problem or the function is time-critical, convert it to a loop.
An option would be to provide s as vector instead of a list of scalars.
Md Nurul Anwar
2022-9-15
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Manage Products 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
