# compose (non-)symbolic functions

I wonder how to compose with a non-symbolic function, or, more generally, how to compose several functions, if some of them are non-symbolic.

In my concrete case, I wish to plot a transformed Gaussian distribution function (concretely, the Camp-Paulson approximation), comparing with a plot of probabilities of the binomial distribution this approximates.

What works, but is not flexible enough is:

import scipy.stats
n=20
p=0.05
binom_dist = scipy.stats.binom(n,p)
Bin=bar_chart([binom_dist.pmf(x) for x in range(n+1)],width=1)
T = RealDistribution('gaussian', 1)
CP=plot(lambda k:T.distribution_function(-1/3*(((0.95*k + 0.95)/(-0.05*k + 1.0))^(1/3)*(1/(k + 1) - 9) + 1/(k - 20) + 9)/sqrt(((0.95*k + 0.95)/(-0.05*k + 1.0))^(2/3)/(k + 1) - 1/(k - 20))),(k,-1,21) , rgbcolor=(0.8,0,0))
show(Bin+CP)


This option has the argument of T.distribution_function hard coded. I rather would have this computed from n and p, using another function. I have two approaches so far, that both don't work.

The first one has a function K computing the input to T.distribution_function like so:

import scipy.stats
n=20
p=0.05
binom_dist = scipy.stats.binom(n,p)
Bin=bar_chart([binom_dist.pmf(x) for x in range(n+1)],width=1)
T = RealDistribution('gaussian', 1)
#
a = 1/(9*n-9*k)
b = 1/(9*k+9)
r = (k+1)*(1-p)/(n*p-k*p)
c = (1-b)*r^(1/3)
μ = 1 - a
σ = sqrt(b*r^(2/3) + a)
def K(k): return (c - μ)/σ
#
CP=plot(lambda k:T.distribution_function(K(k)),(k,-1,21), rgbcolor=(0.8,0,0))
show(Bin+CP)


another alternative tries to set up the composed function before entering plot like so:

import scipy.stats
n=20
p=0.05
binom_dist = scipy.stats.binom(n,p)
Bin=bar_chart([binom_dist.pmf(x) for x in range(n+1)],width=1)
T = RealDistribution('gaussian', 1)
#
a = 1/(9*n-9*k)
b = 1/(9*k+9)
r = (k+1)*(1-p)/(n*p-k*p)
c = (1-b)*r^(1/3)
μ = 1 - a
σ = sqrt(b*r^(2/3) + a)
def distCPapprox(k): return T.distribution_function((c - μ)/σ)
#
CP=plot(lambda k:distCPapprox(k),(k,-1,21), rgbcolor=(0.8,0,0))
show(Bin+CP)


Both approaches don't work. I searched extensively for advice for function composition but couldn't find anything helpful. Your advice is greatly appreciated.

edit retag close merge delete

The following idea seems to go in the correct direction:

import scipy.stats
n=20
p=0.3
binom_dist = scipy.stats.binom(n,p)
Bin=bar_chart([binom_dist.pmf(x) for x in range(n+1)],width=1)
T = RealDistribution('gaussian', 1)
#
def distCPapprox(k):
#
a = 1/(9*n-9*k)
b = 1/(9*k+9)
r = (k+1)*(1-p)/(n*p-k*p)
c = (1-b)*r^(1/3)
μ = 1 - a
σ = sqrt(b*r^(2/3) + a)
return T.distribution_function((c - μ)/σ)
#
CP=plot(lambda k:distCPapprox(k),(k,-1,n+1), rgbcolor=(0.8,0,0))
show(Bin+CP)


And, finally, what I intended seems to be the following. (However, I still would like to know how composition of functions, some of which are non-symbolic works.)

import scipy.stats
n=30
p=0.7
binom_dist = scipy.stats.binom(n,p)
Bin=bar_chart([binom_dist.cdf(x) for x in range(n+1)],width=1)
T = RealDistribution('gaussian', 1)
#
def cumCPapprox(k):
a = 1/(9*n-9*k)
b = 1/(9*k+9)
r = (k+1)*(1-p)/(n*p-k*p)
c = (1-b)*r^(1/3)
μ = 1 - a
σ = sqrt(b*r^(2/3) + a)
return T.cum_distribution_function((c - μ)/σ)
#
CP=plot(lambda k:cumCPapprox(k),(k,-1,n+1), rgbcolor=(0.8,0,0))
show(Bin+CP)


Sort by » oldest newest most voted Your first (original) attempt can be fixed by making k and K symbolic, and using substitution:

import scipy.stats
n=20
p=0.05
binom_dist = scipy.stats.binom(n,p)
Bin=bar_chart([binom_dist.pmf(x) for x in range(n+1)],width=1)
T = RealDistribution('gaussian', 1)
#
var('k')
a = 1/(9*n-9*k)
b = 1/(9*k+9)
r = (k+1)*(1-p)/(n*p-k*p)
c = (1-b)*r^(1/3)
μ = 1 - a
σ = sqrt(b*r^(2/3) + a)
K = (c - μ)/σ
#
CP=plot(lambda k: T.distribution_function(K.subs(k=k)),(k,-1,21), rgbcolor=(0.8,0,0))
show(Bin+CP)


Note that in K.subs(k=k) the first k is the name of a keyword parameter (which will be interpreted by Sage as the name of a symbolic variable), and the second k is the argument which is passed into the lambda function.

Your second (original) attempt can be fixed by making k symbolic and doing the substitution inside distCPapprox:

import scipy.stats
n=20
p=0.05
binom_dist = scipy.stats.binom(n,p)
Bin=bar_chart([binom_dist.pmf(x) for x in range(n+1)],width=1)
T = RealDistribution('gaussian', 1)
#
var('k')
a = 1/(9*n-9*k)
b = 1/(9*k+9)
r = (k+1)*(1-p)/(n*p-k*p)
c = (1-b)*r^(1/3)
μ = 1 - a
σ = sqrt(b*r^(2/3) + a)
def distCPapprox(k): return T.distribution_function(((c - μ)/σ).subs(k=k))
#
CP=plot(distCPapprox,(k,-1,21), rgbcolor=(0.8,0,0))
show(Bin+CP)


Finally, the solution in your own comments is to create a function that takes a number and returns a number (avoiding symbolics entirely); personally I very much prefer that approach if symbolics are not needed.

more