# 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)


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)

( 2022-08-07 17:44:40 +0200 )edit

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)

( 2022-08-07 18:13:24 +0200 )edit

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