Ask Your Question
2

Use cumulative distribution functions

asked 2016-11-03 17:57:25 -0600

Massimo2013 gravatar image

I don't understand how to use cumulative distribution functions in sagemath. For example, I'd like to define a function N(x) as N(x)=1-F(x) where F is a cumulative log-normal distribution However, if I try

 W = RealDistribution('lognormal',[1.5,.6])
 N(x) = W.cum_distribution_function(x)

it does not seem to work at all

edit retag flag offensive close merge delete

3 answers

Sort by ยป oldest newest most voted
2

answered 2016-11-04 09:43:20 -0600

Using the actual function I get a fine plot via:

sage: CDFlogn(x,nu,sigma) = 1/2*(1+erf((log(x)-nu)/(sigma*sqrt(2))))
sage: plot(CDFlogn(x,1.5,.6),0,10)
Launched png viewer for Graphics object consisting of 1 graphics primitive

sage: CDFlogn(x,1.5,.6)
1/2*erf(-0.833333333333333*sqrt(2)*(-log(x) + 1.50000000000000)) + 1/2
sage: CDFlogn(4,1.5,.6)
1/2*erf(-0.833333333333333*sqrt(2)*(-log(4) + 1.50000000000000)) + 1/2
sage: CDFlogn(4,1.5,.6).n()
0.424846794957399

For example: https://en.wikipedia.org/wiki/Log-nor...

edit flag offensive delete link more
2

answered 2016-11-04 03:26:29 -0600

updated 2016-11-04 08:19:57 -0600

Use a Python function, either using def, or using lambda.

This will let you evaluate the function, and plot it.

Here, using lambda:

sage: W = RealDistribution('lognormal',[1.5, .6])
sage: def N(x): return W.cum_distribution_function(x)
sage: N(4)
0.4248467949573992
sage: plot(N, (x, 0, 10))

Here, using def:

sage: W = RealDistribution('lognormal',[1.5, .6])
sage: N = lambda x: W.cum_distribution_function(x)
sage: N(4)
0.4248467949573992
sage: plot(N, (x, 0, 10))

The down side is you need to do that for every variation based on N.

For instance, to plot 1 - N,

sage: NN = lambda x: 1 - N(x)
sage: plot (NN, (x, 0, 10))

Of course, to eg differentiate it, you would really need a symbolic function.

Maybe someone knows how to turn it into a symbolic function.

Maybe it's worth opening a Sage trac ticket to make cumulative distribution functions symbolic.

See https://trac.sagemath.org/wiki/symbolics/functions.

edit flag offensive delete link more

Comments

Thanks. However, this does not seem to allow me to use N(x) in an expression as I would do with any other (symbolic) function. For example, I'm not even able to 'plot(1-N,(x,0,10))' This is unfortunate, as I wanted to translate into sagemath a simulation I did in R, where I can use "proper" functions, such as 'plnorm' and 'dlnorm'

Massimo2013 gravatar imageMassimo2013 ( 2016-11-04 07:02:15 -0600 )edit

I've expanded my answer to show how to plot 1 - N. Maybe you can share your whole R simulation and we can see how well it can translate into SageMath.

slelievre gravatar imageslelievre ( 2016-11-04 08:21:38 -0600 )edit

Here is a shortened version of what I did in R, and would like to replicate in sagemath:

x=(0:80)/10
meanl=5/3
sdl=2/3
n <- function(z) { 1-plnorm(z,meanlog=meanl,sdlog=sdl) }
b <- function(z) dlnorm(z, meanlog=meanl, sdlog=sdl)
BB <- function(z) integrate(function(w) w * b(w), lower = z, upper = Inf)$value
B <- Vectorize(BB, vectorize.args='z')
h <- function(x) 40*exp(-x)
plot(x,h(x)+x,ylim=c(0,10),type="l")
lines(x,5-x-n(h(x))*h(x)+B(h(x)),xlim=c(0,max(x)),col=3)
Massimo2013 gravatar imageMassimo2013 ( 2016-11-04 14:23:56 -0600 )edit
2

answered 2016-11-03 18:36:39 -0600

ndomes gravatar image
W.cum_distribution_function

is not a symbolic function. What you can do is creating N as an alias:

N = W.cum_distribution_function
print N(4)
plot(N)
edit flag offensive delete link more

Comments

But what if I need to use it as/in a symbolic function?

Massimo2013 gravatar imageMassimo2013 ( 2016-11-04 01:01:49 -0600 )edit

Then at the moment you need to create your function in mathematical terms. Please see and rate my answer.

rws gravatar imagerws ( 2016-11-04 09:23:56 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2016-11-03 17:57:25 -0600

Seen: 194 times

Last updated: Nov 04 '16