Hi! I'm trying to approximate probability distributions with weighted polynomials weighted with exp(- x'x). Basically, the idea is that polynomials are "conjugate" (multiply two polynomials, and you get another polynomial), and easy to integrate -- two things that could make Bayesian inference very easy. So I'm now trying to write a few things to compare how well it does in accuracy & speed vs MCMC.
The integration is where Reduce comes in. It lets me generate functions for marginalizing out parameters or CDFs without having to write out all the different expressions for the integral of a*exp(-x^2), b*x*exp(-x^2), c*x^2*exp(-x^2), etc....
Two changes I would like to move towards eventually are splines in place of polynomials, and fitting log densities instead of densities. There's way too much math I don't understand before I can get there, so even though both of those things (especially splines) would be a massive improvement, it'll be months before I get there. Just trying a proof of concept now.
rcall(:(int( x^3 * exp(-x^2), x))). It's not well documented yet, but you can also get the same result by doing this: int(:(x^3 * exp(-x^2)), :x).
Cool, I like that syntax more.
Definitely looking forward to seeing where it goes. Also depends on how powerful it is as a CAS.
For example, I'd eventually like to write a macro that lets people specify a statistical model. I could see Reduce being helpful, potentially simplifying user-provided expressions, or even marginalizing out nuisance parameters entirely.
That'd be more accurate than doing it numerically, but it'd be complicated to implement and would have to be pretty efficient to actually save time (don't want things to take an eternity to compile).
"there is still a lot of functionality that can be built on top of it. For now, I just wanted to make sure that the basic functionality is robust, but more features will come."
What sort of things do you have in mind?
Yes, the compile-time is currently slower than I would prefer, but once the parsers are compiled the performance is quicker. This is because I had to add in features to the parser to make it more robust for handling more complicated nested expressions, functions, begin/end blocks (so the code for it is large). I did this at the expense of initial compilation time. It's possible to improve the performance of it, but it requires a lot of concentration since it grew to hundreds of lines of code. I have some ideas for improving the performance but it will not happen over night.
You can take a look at the Reduce documentation for more information on its capabilities: http://www.reduce-algebra.com/manual/manual.html
Recently I implemented the parsing of matrices and there are some REDUCE packages that manipulate matrices, so I would like to build some more of those features into the Julia package. also, I would like to try out some of the differential equations solving capability in the future. I don't have a specific road plan laid out, but as the need for more features arises and as time passes, it is becoming more clear how it can be integrated into Julia.
The REDUCE software is actually one of the fastest computer algebra systems on the planet, unfortunately the text communcation via Pipe and parsing into Julia aspect of it slows it down significantly; but it still might be faster overall than some of the other symbolic packages available for Julia. One notable difference is that it doesn't define a new "symbolic type" within Julia, instead the package works on Julia expression objects.
Julia does not have any fully native computer algebra packages, so this can help work around that. I wouldn't recommend using the package for doing things that you could otherwise do natively with Julia; but it should be a good tool for those situations where an external symbolic computation is essential for code generation.
julia> rcall(:(int( x^2 * exp(-x^2/2 ) , x, -Inf, Inf)))
:(sqrt(π) * sqrt(2))
julia> int(:( x^2 * exp(-x^2/2 ) ), :x, -Inf, Inf)
:((sqrt(π) * e ^ (inf ^ 2 / 2) * sqrt(2) * erf(inf / sqrt(2)) - 2inf) / e ^ (inf ^ 2 / 2))Seems to be a bug. Besides the second not being what I'm looking for, "inf" isn't defined, so it also throws errors.
Although, most people may actually prefer errors over Infs, so they know asap something didn't go quite right.
Also, while Reduce itself doesn't need SpecialFunctions.jl, a lot of people may want it. Not sure for how long the helpful error will stick around, so may be worth at least mentioning it somewhere in the docs.
Was hoping going back to rcall would fix this:
julia> ex2 = :(exp(-a*exp(-x)-b*x)*x^2)
:(exp(-a * exp(-x) - b * x) * x ^ 2)
julia> int(ex2,:x,-Inf,Inf)
:(((((((((((e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=inf, int(x ^ 2 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b ^ 2 + 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=inf, int(x / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b + 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=inf, int(1 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a) - e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=-inf, int(x ^ 2 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b ^ 2) - 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=-inf, int(x / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a * b) - 2 * e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * sub(x=-inf, int(1 / e ^ ((e ^ x * b * x + e ^ x * x + a) / e ^ x), x)) * a) + e ^ ((2 * e ^ inf * b * inf + a) / e ^ inf) * b ^ 2 * inf ^ 2) - 2 * e ^ ((2 * e ^ inf * b * inf + a) / e ^ inf) * b * inf) + 2 * e ^ ((2 * e ^ inf * b * inf + a) / e ^ inf)) - e ^ (e ^ inf * a) * b ^ 2 * inf ^ 2) - 2 * e ^ (e ^ inf * a) * b * inf) - 2 * e ^ (e ^ inf * a)) / (e ^ ((e ^ (2inf) * e ^ inf * b * inf + 2a) / e ^ inf) * b ^ 3))But:
julia> rcall(:(int(exp(-a*exp(-x)-b*x)*x^2),x,-Inf,Inf))
ERROR: Reduce: Nested :tuple block structure not supported
Stacktrace:
[1] show_expr(::Base.AbstractIOBuffer{Array{UInt8,1}}, ::Expr) at /home/chris/.julia/v0.6/Reduce/src/parser.jl:505
[2] unparse(::Expr) at /home/chris/.julia/v0.6/Reduce/src/parser.jl:547
[3] Type at /home/chris/.julia/v0.6/Reduce/src/rexpr.jl:38 [inlined]
[4] #rcall#14(::Array{Symbol,1}, ::Array{Symbol,1}, ::Function, ::Expr) at /home/chris/.julia/v0.6/Reduce/src/rexpr.jl:292
[5] rcall(::Expr) at /home/chris/.julia/v0.6/Reduce/src/rexpr.jl:292
julia> rcall(:(int($ex2,x,-Inf,Inf) ))
:(int(exp(-a * exp(-x) - b * x) * x ^ 2, x, -inf, inf))The answer to that involves the di- and tri-gamma functions, and requires a > 0 && b > 0.
"Julia does not have any fully native computer algebra packages, so this can help work around that. I wouldn't recommend using the package for doing things that you could otherwise do natively with Julia; but it should be a good tool for those situations where an external symbolic computation is essential for code generation."
That's probably true for generated functions in general. Your code would probably compile faster if instead you wrote out explicit versions for all of them.
julia> int(:( x^2 * exp(-x^2/2 ) ), :x, "-infinity", "infinity")
:(sqrt(π) * sqrt(2))
julia> rcall(:(int( x^2 * exp(-x^2/2 ) , x, -Inf, Inf)))
:(sqrt(π) * sqrt(2))As for the Nested :tuple block structure not supported error, it seems you made a mistake with your parenthesis, it should have been entered like this:
julia> rcall(:(int(exp(-a*exp(-x)-b*x)*x^2,x,-Inf,Inf)))
:(int(exp(-a * exp(-x) - b * x) * x ^ 2, x, -inf, inf))
julia> int(ex2,:x,"-infinity","infinity")
:(int(exp(-a * exp(-x) - b * x) * x ^ 2, x, -inf, inf))which also appears to be an unsupported integral (from the upstream REDUCE end) at first glance
when an integral cannot be evaluated by REDUCE, it just returns the expression back
However, the REDUCE documentation does describe how to set the positivity of variables:
http://www.reduce-algebra.com/manual/manualse68.html#x83-15900010.2
On the di- and tri-gamma functions, I see your Julia docs already mention support for polygamma. Polygamma(0,x) == digamma(x) and polygamma(1,x) == trigramma(x).
So the answer to the likely unsupported integral is: a^-b*Gamma(b)*( (log(a)-polygamma(0,b))^2 + polygamma(1,b) ).
Not sure how many special functions like the bessel, hypergeometrics, etc, it does support. Not sure how likely Reduce developers are to be able to increase support for things like that.
Likely unsupported, because it requires a > 0 && b > 0. Perhaps it could return an answer, given those conditions. From the Wolfram sandbox:
https://sandbox.open.wolframcloud.com/app/objects/384d34c1-ac46-46a5-ad28-b91ed2876d5e#sidebar=computeAssuming[a > 0 && b > 0, Integrate[Exp[-(a Exp[-x]) - b x] x^2, {x, -Infinity, Infinity}]]
Is where I got the answer. If I leave off the constraint, it'll return the same answer within a conditional expression.
Not sure how well Reduce handles conditionals; I'll get around to looking through the docs more thoroughly as I try to use it more. For now, http://www.reduce-algebra.com/manual/manualse16.html#x25-380005.3 appears to only cover user-provided conditionals.
Also, I'm glad I can pass reduce strings as a fall back. Trying to follow the example on the conditionals:
julia> rcall( :(abs(a*b*c)), "let sign(a) => 1, sign(b) => 1" )
ERROR: Reduce:
***** nil not defined as switch
abs(c)*a*b
***** nil not defined as switch
Stacktrace:
[1] ReduceCheck(::String) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/Reduce.jl:62
[2] read(::Reduce.PSL) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/Reduce.jl:84
[3] readsp(::Reduce.PSL) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/Reduce.jl:88
[4] #rcall#13(::Array{Symbol,1}, ::Array{Symbol,1}, ::Function, ::Reduce.RExpr) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/rexpr.jl:245
[5] (::Reduce.#kw##rcall)(::Array{Any,1}, ::Reduce.#rcall, ::Reduce.RExpr) at ./<missing>:0
[6] #rcall#14(::Array{Symbol,1}, ::Array{Symbol,1}, ::Function, ::Expr) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/rexpr.jl:278
[7] (::Reduce.#kw##rcall)(::Array{Any,1}, ::Reduce.#rcall, ::Expr) at ./<missing>:0
[8] rcall(::Expr, ::String, ::Vararg{String,N} where N) at /home/celrod/Downloads/JuliaPro-0.6.0.1/JuliaPro/pkgs-0.6.0.1/v0.6/Reduce/src/rexpr.jl:282So, I see Reduce itself returned the right answer, but the parser failed.
rcall incorrectly there. Extra arguments to rcall have to be REDUCE switch names passed as Julia symbols, you can't just send it any random string. That is why the REDUCE error tells you that you did not enter a defined switch. For int and df it's different, since those can accept extra args. You have to make 2 separate calls to rcall if you want to first evaluate a REDUCE string and then a Julia expression object like this:julia> rcall("let sign(a) => 1, sign(b) => 1")
""
julia> rcall( :(abs(a*b*c)))
:(abs(c) * a * b)julia> solve(((:m - :x) ^ 2 - :s) / :s ^ 2, :s)
1-element Array{Any,1}:
:(s = (m ^ 2 - 2 * m * x) + x ^ 2)