# Computation of the expected value of a multidimensional array providing PMF on each dimension

**URL:** <https://discourse.julialang.org/t/computation-of-the-expected-value-of-a-multidimensional-array-providing-pmf-on-each-dimension/94404>\
**Category:** General Usage\
**Created:** [February 10, 2023, 10:20am UTC](https://discourse.julialang.org/t/computation-of-the-expected-value-of-a-multidimensional-array-providing-pmf-on-each-dimension/94404 "2023-02-10T10:20:26Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [February 10, 2023, 10:20am UTC](https://discourse.julialang.org/t/computation-of-the-expected-value-of-a-multidimensional-array-providing-pmf-on-each-dimension/94404/1 "2023-02-10T10:20:26Z")

</div>

Isn’t there something smarter already in the Julia standard lib to compute the expected value of a multidimensional tensor given a vector of PMFs than the following code ?

```julia
julia> function EV(v,p) # v is a multidimensional array, p is a vector of vectors (PMFs, one per dimension)
           (ndims(v) == length(p) && all(size(v) .== length.(p))) || error("Mismatch dimension or size between the value tensor and the probability vectors")
           (all([all(>=(0-eps()), v) for v in p]) && all(sum.(p) .> 1 .- eps() ) && all(sum.(p) .> 1 .- eps() )) || error("p is not a vector of probabilities")
           outsum = 0.0
           comp = 0.0 #https://en.wikipedia.org/wiki/Kahan_summation_algorithm
           nd = ndims(v)
           for idx in CartesianIndices(size(v)) 
               y = *(v[idx],[p[d][idx[d]] for d in 1:nd]...) - comp
               t = outsum + y
               comp = (t - outsum) - y
               outsum = t
           end
           return outsum
       end
EV (generic function with 1 method)

julia> v = collect(reshape(1:12,(3,2,2)));

julia> p = [[0.5,0.3,0.2],[0.6,0.4],[0.8,0.2,]];

julia> EV(v,p)
4.1000000000000005

julia> v = [2;2;;2;2;;;4;4;;4;4;;;]
2×2×2 Array{Int64, 3}:
[:, :, 1] =
 2 2
 2 2

[:, :, 2] =
 4 4
 4 4

julia> p = [[0.5,0.5],[0.5,0.5],[0.8,0.2]]
3-element Vector{Vector{Float64}}:
 [0.5, 0.5]
 [0.5, 0.5]
 [0.8, 0.2]

julia> EV(v,p)
2.4000000000000004

```