# Likelihood of a set of events originating from a set of sources

**URL:** <https://discourse.julialang.org/t/likelihood-of-a-set-of-events-originating-from-a-set-of-sources/101455>\
**Category:** Modelling & Simulations\
**Tags:** numerics, bayesian-inference\
**Created:** [July 10, 2023, 9:31pm UTC](https://discourse.julialang.org/t/likelihood-of-a-set-of-events-originating-from-a-set-of-sources/101455 "2023-07-10T21:31:25Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [July 10, 2023, 9:31pm UTC](https://discourse.julialang.org/t/likelihood-of-a-set-of-events-originating-from-a-set-of-sources/101455/1 "2023-07-10T21:31:25Z")

</div>

I’m hoping there’s a more performant way to do this, but I’m missing something.

I’m migrating some MATLAB code that uses reversible-jump MCMC to produce optimal sets of point sources to explain observed events (2d coordinates + precision).

The way a model is evaluated for a move to be accepted in MATLAB is to generate a matrix of all combinations of sources and coordinates, with each value being the normpdf described by the modeled source coordinate and event precision evaluated at the event coordinate. This is done separately for each axis and multiplied, and then the values for each localization are summed across sources, then the log of these are summed to produce the score for the model.

This is what the code looks like in Julia, partially optimized:

```julia
    xworking = zeros(μxs)
    yworking = zeros(μys)
    ws = 1 / nemitters
    p_alloc = 0.0
    for i ∈ eachindex(xs)
        xworking .= normpdf.(μxs .+ αxs .* frames[i], σxs[i], xs[i])
        yworking .= normpdf.(μys .+ αys .* frames[i], σys[i], ys[i])
        xworking .*= yworking
        p_emitter = sum(xworking)
        p_alloc += log(p_emitter * ws)
    end
    return p_alloc
end

```

(the alpha values allow for linear drift)

Profiling shows `normpdf`, within the broadcast `copy` call, is the bottleneck.

I feel like there must be a more efficient approach to this. Are there particular packages I should look into? Can KernelDensity.jl be used?

The original MATLAB code is here for reference:

```julia
function LogL = p_Alloc(SMD,Mu_X,Mu_Y,Alpha_X,Alpha_Y,Ws)
%This function calculated the probability of a given allocation set.
    X=SMD.X;
    Y=SMD.Y;
    T = repmat(SMD.FrameNum,[1,length(Mu_X)]);
    SigmaX=SMD.X_SE;
    SigmaY=SMD.Y_SE;
    
    Lx = length(X);
    Lmu = length(Mu_X);
    LogL = log(sum(repmat(Ws,[Lx,1]).*normpdf(repmat(X,[1,Lmu]),...
            repmat(Mu_X,[Lx,1])+repmat(Alpha_X,[Lx,1]).*T,...
            repmat(SigmaX,[1,Lmu])).*normpdf(repmat(Y,[1,Lmu]),...
            repmat(Mu_Y,[Lx,1])+repmat(Alpha_Y,[Lx,1]).*T,...
            repmat(SigmaY,[1,Lmu])),2));
   
    LogL = sum(LogL);
    
end

```