I want to use Julia to get the end result like this:

the Expectation of (X - \mu)^{2} will be computed at the end to become:

E(X^{2}) - \mu^{2}

but how?

this is my try out:

using SymPy
# https://docs.sympy.org/latest/tutorials/intro-tutorial/simplification.html
@syms x, μ
v = expand((x-μ)^2)
function E(h)
if h == x^2
return x^2
elseif h==x*μ
return x*μ
elseif h==μ^2
return μ^2
end
end

what should I do next? I was thinking of creating a function with if and ifelse, maybe there is another way?

Not sure if I understand what you want to do in the end:

Compute the variance of a random variable based on a sample? If you have a sample of the random variable X and want to know/estimate the variance, there would be tools for that in Statistics · The Julia Language And probably lots more.

(I guess this is what you want) Represent random variables symbolically in a program? E.g. the operation \mathbb{E} which represents the mean of a random variable with the properties you described (e.g. \mathbb{E}(\mathbb{E}(X)) = \mathbb{E}(X) ). I’m not sure if there is a package that does this already, but it seems like sympy does have types that represent random variables – doesn’t seem like they have been ported to SymPy.jl yet though.

If you just want to implement a “rewriting rule” that takes an expression like E((x - mu)^2) and returns E(x^2) - mu^2 this might also be helpful: SymbolicUtils.jl

Perhaps you could define a symbolic function E and specify some rewrite rules.

It should work. Did you start a new REPL session when trying the code?
Both SymPy.jl and SymbolicUtils.jl export @syms, so what happened in your case is probably that z is a Sympy symbol and thus the rewriting does not work (it looks for SymbolicUtils symbols in the expression).